Pure cross-anisotropy for geotechnical elastic potentials

The pure cross-anisotropy is understood as a special scaling of strain (or stress). The scaled tensor is used as an argument in the elastic stiffness (or compliance). Such anisotropy can be overlaid on the top of any elastic stiffness, in particular on one obtained from an elastic potential with its own stress-induced anisotropy. This superposition does not violate the Second Law. The method can be also applied to other functions like plastic potentials or yield surfaces, wherever some cross-anisotropy is desired. The pure cross-anisotropy is described by the sedimentation vector and at most two constants. Scaling with more than two purely anisotropic constants is shown impossible. The formulation was compared with experiments and alternative approaches. Static and dynamic calibration of the pure anisotropy is also discussed. Graphic representation of stiffness with the popular response envelopes requires some enhancement for anisotropy. Several examples are presented. All derivations and examples were accomplished using the algebra program Mathematica.


Introduction
Elastic response is an essential part of most constitutive models for soils. It is particularly important for soil dynamics, for stability analysis [2] and for material response in the range of small strains. This range corresponds roughly to strain amplitudes of 10 À5 for sand and 10 À4 for clays. Under such loading soil can be much stiffer than at amplitudes of say 10 À3 . This paper deals with small-strain elastic (incrementally linear) stiffness only. For larger amplitudes, hysteretic [23] or cumulative models [24] are necessary. Stiffness may be a function EðrÞ of stress (or strain), but it interrelates rates (or tiny increments) of stress and strain rather than stress and strain themselves.
In the elastic regime, stress should be a continuous 1-1 function rðeÞ of strain. Otherwise, some stress could be accumulated within a closed strain loop, see Sect. 2. A thermodynamically sound elastic material model should not allow for the accumulation of stress or energy upon any closed strain loop. The energetic requirement is not trivial for soils with a barotropic (pressure-dependent) stiffness. It is well known that the barotropic elastic modulus, E $ p or E $ ffiffi ffi p p , with a constant Poisson number m violates the Second Law [13,31]. In order to avoid this problem, several elastic potentials have been proposed in the literature, see Sect. 2.1. A tangential stiffness obtained from such potential is a function of stress (not only of stress invariants), and one may speak of the stress-induced anisotropy 1 (rA). It should be distinguished from the inherent cross-anisotropy 2 (ÂA), which is caused by sedimentation process and/or geological petrification (cementation) of the geostatic K 0 state. The ÂA is independent of the current stress or strain.
Any constant cross-anisotropic stiffness E ÂA ijkl can be described by five material constants, usually denoted as E v ; E h ; m h ; m vh and G v , see Sect. 4. The main objective of this paper is to represent this stiffness in the form 3 wherein the elastic properties 4 are given in the isotropic stiffness E iso and all pure anisotropic properties are moved to the anisotropy tensor Q. The advantage of such separated description follows from the fact that the same Q can be applied to any hyperelastic (and barotropic) stiffness without violating the Second Law. This is proven in Sect. 3. In other words, any basic tangential stiffness (or compliance), possibly with its own induced anisotropy, can be superposed by the pure inherent anisotropy. Here, this pure cross-anisotropy is denoted as ÂA M wherein M is the number of constants required for the anisotropy tensor 5 Q. Two anisotropy tensors Q, for ÂA 1 and ÂA 2 , are analytically derived in Sects. 5 and 6. Unfortunately, the derivation of Q for the general case ÂA 3 is not feasible as demonstrated in Sect. 7.
Calibration of the parameters of Q from static (cyclic) triaxial tests on samples cut in different directions or from wave velocities in different directions [8,27] is commented in Sect. 8. A few remarks on experimental data for ÂA are given in Sect. 9, and the advantage of ÂA 2 over ÂA 1 is demonstrated.
The graphic representation of stiffness in the form of polar response envelopes [11] is well known in the geotechnical literature. In the case of ÂA, some complications may arise from the fact that the stress rate, _ rðr 0 ; _ e; MÞ, may not be axisymmetric for the axisymmetric initial stress, r 0 , and co-axisymmetric 6 strain rate, _ e. The problem is caused by the dependence on the direction of sedimentation, m, appearing here in the form of the sedimentation dyad, M ¼ m m. This may also cause a loss of coaxiality. Therefore, an enhanced graphic representation is proposed in Sect. 10. Some examples of extended response envelopes with ÂA 2 and polar diagrams of wave velocities are shown.
Finally, ÂA 2 is applied to stress and substituted to the Matsuoka-Nakai yield surface. The modified surface is shown graphically in Sect. 11. All relevant packages and notebooks for the algebra program Mathematica are available from the authors.

Notation
Bold-face letters like r are vectors or second rank tensors. Sans serif letters, e.g., E, are the fourth-order tensors.
Gibbs notation like _ r ¼ E : _ e or index notation _ r ij ¼ E ijkl _ e kl in the Cartesian coordinate system with usual summation over repeated (dummy) indices is used. The geotechnical sign convention is applied to r and e with compression positive. A fourth-order tensor E can appear in a form of a 9 Â 9 matrix (no Voigt 6 Â 6 notation) denoted as [E]. The 9 Â 9 form facilitates some transformations in the algebra program Mathematica. Similarly, ½r is the 3 Â 3 matrix obtained from the tensor r. The essential variables are: between the stress rate _ r ij and the strain rate _ e kl . The tangential stiffness E ijkl needs not be constant. It may be a function of stress or strain, but it cannot be a function of their rates. Such incrementally linear model is called hypoelastic.
Let the strain evolve along the path 7 e ij ðsÞ, Fig. 1a. After a 180 reversal, identical negative strain increments can be applied in the opposite sequence and the strain evolves back along exactly the same path. The relation _ r ij ðÀ _ e kl Þ ¼ À _ r ij ð _ e kl Þ holds due to the incremental linearity. Hence, the same stress path is followed and, eventually, the original state r ij ðt 0 Þ is reached. The energy density, dW ¼ r ij _ e ij dt, is also recovered. However, if one departs from e ij ðt 0 Þ upon one path and returns to e ij ðt 0 Þ upon another path, Fig. 1b, then neither the initial stress nor the energy is in general recovered. At least, one cannot conclude such recovery from incremental linearity (2) alone.
In hyperelastic models, apart from linear relation (2), some additional conditions must be imposed on E ijkl . In isothermal elastic materials, strain is the only independent state variable, i.e., e ij alone dictates the internal elastic energy W. This dependence must be a function WðeÞ, i.e., the elastic energy cannot depend on the strain path e ij ðsÞ. The change in W upon the path from e 0 and this DW is identical upon any strain path e ij ðsÞ. If the choice of a path e ij ðsÞ between e 0 ij and e 1 ij could influence the integral DW, then one could input less energy upon one path, 0 ! 1, than could be recovered on the way back, 1 ! 0. Such gain of energy without any change of state (strain returns to e 0 ij ) violates the Second Law. Even if this gain occurred at the cost of thermal energy, it would be a violence of the Second Law (a perpetuum mobile of the second kind). Hence, the integral in (3) should indeed be path-independent, which implies the existence of a function WðeÞ. Being a function, WðeÞ has the total differential From the comparison of (4) with (3) for any de ij , it follows that As a derivative of a function of strain, stress also must be a function rðeÞ. Stress rate can be calculated using the chain rule, _ r ij ¼ ðor ij =oe kl Þ _ e kl . From the comparison with (2) It is evident from (6) 2 that E ijkl must be symmetric. Note, however, that the symmetry, E klij ¼ E ijkl , is only a necessary (but not sufficient) condition for the existence of an elastic potential. Let a symmetric stiffness E klij ðeÞ be a primary function. For the existence of WðeÞ, also a function r ij ðeÞ must exist. For the integrability all mixed second derivatives of r ij ðe kl Þ must be identical which is not guaranteed by the symmetry E klij ¼ E ijkl . For example, E ijkl ðeÞ ¼ e nn 3Kmd ij d kl =ð1 þ mÞ þ 2GI ijkl Â Ã is symmetric, but it is not hyperelastic because it does not satisfy condition (8).
Functions WðeÞ cannot be directly measured. They are usually formulated by trial and error. An educated guess can be based on the measurements of second derivatives E ijkl (6) 2 at different strains. Alternatively, the complementary energy WðrÞ may be used, In granular materials, the main difficulty in the formulation of WðeÞ or WðrÞ arises from the pressure dependence (the so-called barotropy) of the stiffness.

Geotechnical hyperelastic models
Several hyperelastic models have been proposed in the literature. A critical review can be found in [20] and more recently in [9]. It is helpful to assume the hyperelastic stiffness as a homogeneous function of stress, i.e., 8k [ 0 : EðkrÞ ¼ k m EðrÞ. The order m of homogeneity is usually m % 0:6 for sand and m % 1 for clays. The compliance, C ¼ E À1 , is homogeneous of order Àm, of course. It can be proven 8 that the corresponding elastic potentials, WðrÞ and WðeÞ, are homogeneous functions of order 2 À m and ð2 À mÞ=ð1 À mÞ, respectively.
A simple hyperelasticity was proposed by Vermeer [28]. The hyperelastic potential is given explicitly, with a material constant c 1 . The order of homogeneity of EðrÞ must be m 6 ¼ 1.
Borja et. al [4] proposed a hyperelastic model based on elastic potential formulated in terms of the strain invariants, wherein e Ã is the deviatoric part of e. In this case, the stiffness appears to be inhomogeneous in stress. Niemunis and Cudny [20] introduced a potential for clays that yields stiffness EðrÞ with a homogeneity of order m ¼ 1.
The following expression for the complementary energy was proposed for sand by Niemunis et al. [21] WðrÞ ¼ c 11 P c 12 R 2ÀmÀc 12 ; wherein m 6 ¼ 1 is the order of homogeneity of EðrÞ.
Response envelopes [11] are polar representations of stiffness at different stresses, see Sect. 10. They can be measured (here for medium dense sand [14,15]) and calculated analytically, e.g., using (13). A comparison like in Fig. 2 may be used for the calibration.
Selected terms from (12) and (13) have been recently combined for kaolin by Gehring [9] into WðrÞ ¼ c 11 P c 12 R 2ÀmÀc 12 þ c 13 P lnðPÞ: This potential is suitable for cohesive materials because the second summand removes the singularity of C at m ¼ 1.
Experimental (for kaolin [9]) response envelopes are compared with the theoretical ones obtained with (14), Fig. 3. A strong inherent anisotropy was caused by K 0 consolidation of kaolin. The required anisotropy tensor Q given in (27) is described in Sect. 5. The proposed superposition of rA and ÂA is a convenient alternative to a direct postulation of Wðr; MÞ with the sedimentation dyad M ¼ m m as an additional argument. For example such function was proposed by Cudny and Staszewska [7] for m 6 ¼ 1.
Similar approach related to the microscopic description has been recently proposed by Amorosi, Houlsby and Rollo [1,12]. Instead of using an explicit potential WðrÞ, Boyce [5] postulated a 1-1 homogeneous function eðrÞ of order 1 À m. In this case, existence of the complementary elastic potential WðrÞ should be proven. For such formulation, the superposition described in the next sections can also be applied using identical tensor Q.

Anisotropy tensor Q
Stiffness E ijmn and a family of transformations E 0 ijmn ¼ a ik a jl a mr a ns E klrs with directional cosines a ij build a symmetry group, if the components of stiffness are preserved, that is, if E 0 ijmn ¼ E ijmn . For an isotropic stiffness E iso ijmn , it is true for any a ij . For an inherent cross-anisotropic stiffness E ÂA ijmn with sedimentation direction m ¼ 0; 0; 1 f g, a ij corresponds to an arbitrary rotation 9 around m by angle w, In this paper, the pure inherent cross-anisotropy ÂA in a form of tensor Q is proposed. It is a function of m and some constants. This ÂA can be ''added'' to any stiffness, e.g., to one obtained from a potential WðeÞ or WðrÞ with its own rA, see Sect. 2.1. The constants in Q can be determined from the transformation 8 For this purpose, one may use ð2 À mÞð1 À mÞ WðrÞ ¼ r : o 2 W oror : r ¼ r : C : r; which is analogous to the well-known Euler formula for homogeneous functions, here applied twice to WðrÞ. The homogeneity of WðrÞ of order 2 À m is sufficient (but not necessary) for the homogeneity of order m in EðrÞ. After adding a constant to WðrÞ, the homogeneity of WðrÞ is lost, but homogeneity of EðrÞ is preserved.
of the isotropic stiffness E iso abcd to the desired E ÂA ijkl . Tensor Q should scale any stiffness in a similar manner. All components of Q are independent of e ij ; E and m, and hence, Q stores the pure anisotropy.
Let us apply Q to the strain, e ij ¼ Q ijkl e kl , and then substitute e ij into an elastic potential Wð eÞ. Differentiating Wð eÞ with respect to e ij and using the chain rule, one obtains the stiffness with the combined effect of rA and ÂA, wherein E rA abcd is the stiffness with rA only. Note that deviations from isotropy are superposed and hence, the symmetry group is restricted rather than extended. Tensors Q have relatively simple forms for ÂA 1 and ÂA 2 with the major symmetry, Q ijab ¼ Q abij , see Sects.5 and 6.
Inverting both sides of (17), one may use Q À1 ijkl for the compliance 10 , The same Q À1 ijkl can be applied to stress, r ij ¼ Q À1 ijkl r kl , and the modified stress r ij can be substituted into the given complementary potential Wð rÞ. Differentiating with the chain rule, one obtains the compliance with superposed effects of rA and ÂA, wherein C rA abcd is the compliance with rA only. Summing up, the most important advantage of the pure anisotropy is the fact that it can be ''added'' a posteriori to any hyperelastic stiffness E rA or compliance 11 C rA without violating the Second Law. Moreover, a fairly easy implementation of Q to existing constitutive models can be expected. Tensor Q can be interpreted as a modifier of the strain tensor 12 In the case of ÂA 1 , a special form of Q derived in Sect. 5 allows to interpret this strain transformation as scaling of the displacements u i and the coordinates x i . This has already been observed by Lodge [17] and used for scaling of boundary value problems. Contrarily to the current approach, Lodge started by scaling of displacements u and coordinates x, which imposes an unnecessary constraint on the scaling of strains e. For example, the anisotropy ÂA 2 cannot be squeezed into the class of anisotropic elastic solids discussed in [17], see Sect. 6.
A different cross-anisotropic scaling was proposed by Osinov and Wu [25]. They applied a diagonal fourth rank tensor P to the resulting hypoplastic stress rate _ r as follows Our tensor Q could be applied to r, i.e., to the argument in EðrÞ in (21). The thermodynamic aspects of P : E were ignored in [25].

Cross-anisotropic constant stiffness
It is well known that constant (stress-independent) crossanisotropic elastic stiffness (22) requires five material constants, E v ; E h ; m h ; m vh and G v . The vertical coordinate is x v (=direction of sedimentation) and the horizontal coordinate is x h , Fig. 4. These material constants will be separated into two elastic parameters and three purely anisotropic ones. This pure anisotropy is denoted as ÂA 3 . For The elastic Young's moduli along x h and x v are E h and E v , respectively. Shear modulus in horizontal plane is G h ¼ E h =ð2ð1 þ m h ÞÞ and from symmetry follows Stability of the material behavior requires elastic stiffness matrix to be positive definite. This implies the following conditions on the material constants The pure anisotropy tensor Q corresponding to ÂA 3 is discussed in Sect. 7 after the presentation of ÂA 1 and ÂA 2 in Sects.5 and 6.

Anisotropy tensor for · A 1
A three-constant elastic cross-anisotropic stiffness has been proposed by Graham and Houlsby [10] using the anisotropy parameter a in the following relations The single parameter a relates the material constants in the horizontal, t h , and in the vertical (parallel to sedimentation), t v , direction. The representation of stiffness for m ¼ The total number of independent material constants is reduced from five to three: Two constants describe the isotropic elasticity and just one pertains to the pure anisotropy, and hence the notation ÂA 1 . Separation of the material constants is essential. Conversion of the isotropic stiffness E iso into ÂA 1 has been only mentioned in [10] without giving an explicit form. Anisotropy tensor Q has been recently derived in [21], viz.
Tensor Q for ÂA 1 depends on m and a only. In the special case of a ¼ 1, the anisotropy tensor is reduced to identity tensor d ik d jl . Due to the symmetry l ij ¼ l ji , the major symmetry holds. Note that l ij transforms e kl into e ij analogously as the directional cosines a ij do, i.e., e ij ¼ l ik l il e kl , see Sect. 3. Hence, l ij could be used to scale the displacements u i or the coordinate axes x i .

Anisotropy tensor for · A 2
It is argued [8,19] that ÂA 1 is overly restrictive. Therefore, an ÂA 2 with two anisotropy constants, a and b, is proposed. These constants provide more flexibility for modelling of pure anisotropy. For b ¼ 1, the ÂA 1 is recovered and for a ¼ b ¼ 1, the tensor Q ijkl is reduced to the identity. The new parameter b is added to (25) as an exponent, Two isotropic elastic parameters, E ¼ E v and m ¼ m h , are supplemented by two anisotropy constants, a and b. For such ÂA 2 , an anisotropy tensor Q ijkl must be found. If applied to constant isotropic elasticity, the resulting stiffness E ÂA2 ijkl ¼ Q abij E iso abcd Q cdkl should be with the same A, B as defined in (26) and X ¼ a 1=b , h ¼ a 2bÀ1 .
By trial and error, the following anisotropy tensor has been found and a, b, c are functions of the constants a and b, namely The major symmetry Q ijkl ¼ Q klij is preserved due to symmetry l ik ¼ l ki given in (32). For m ¼ 0; 0; 1 f g, tensor Q ijkl can be represented as a diagonal matrix and easily 13 inverted to Q À1 ijkl . Otherwise, the analytical inversion requires diagonalization 14 . The new exponent b does not affect stability condition (29). Assuming b ¼ 1 in (30), the ÂA 1 given in (25) is recovered.
The improved flexibility of ÂA 2 goes at the expense of more complex calibration. One possibility is to assume the value of b from the literature, see Sect. 9.
The class of anisotropic elastic solids proposed by Lodge [17] was based on individual scaling of displacements and coordinates. This led to e ij ¼ a ir b js e rs . Our relation e ij ¼ Q ijrs e rs with Q ijrs from (32) cannot be brought to the same form. This fact can be demonstrated using the transposition U ikjl ¼ Q ijkl . There are two nonzero eigenvalues of U, which precludes U from being a dyad.

No pure anisotropy tensor for · A 3
Boehler and Sawczuk [3] formulated the following general representation of isotropic tensorial function of two arguments We need e ¼ Fðe; MÞ to be linear with respect to e because Q ¼ o e=oe should be independent of e. Hence, (35) can be reduced to the following bilinear function wherein only f 0 and f 1 may depend on invariants tre and trðM Á eÞ, i.e., with six material constants C i . The derivative of the stress rate function _ r ¼ Fð _ e; MÞ in representation (38) leads to the linear stiffness wherein C 2 ¼ C 3 follows from the symmetry E ijkl ¼ E klij .
In our case, function e ¼ Fðe; MÞ in representation (38) is differentiated to Q ¼ o e=oe keeping C 2 6 ¼ C 3 , i.e., the tensor Q has the matrix form Of course, (40) holds for m ¼ 0; 0; 1 f g only. With (40) in hand, one may attempt to find the constants C i , for which the postulated separation of elasticity and pure anisotropy is valid. Although the matrices E ÂA3 and E iso are congruent, it can be shown that the separation of elastic constants, E ¼ E v ; m ¼ m h , and purely anisotropic constants, a; b; c, from is not possible using Q given in (40). In order to demonstrate this fact, it is convenient to investigate the compliances, C iso and C ÂA3 , rather than the stiffnesses, E iso and E ÂA3 . For the special case of E ¼ 1, the constant isotropic compliance matrix is  43) and (44), should be coupled analogously to (41). Such coupling is possible, if a set of components of the inverse anisotropy matrix ½Q À1 can be found that satisfies The inverse matrix ½Q À1 has identical formal representation (40) as ½Q. The uniqueness of the solution is not necessary. The following guess ð46Þ nearly satisfies (45). Using ½Q À1 given in (46), the product ½Q À1 T Á ½C iso Á ½Q À1 (45) is almost identical as ½C ÂA3 given in (44). Only one component of ½Q À1 T Á ½C iso Á ½Q À1 differs from the respective component of ½C ÂA3 . These components may be set equal, x 2 ¼ a 1=c , which leads to c ¼ b=2, but this corresponds to the constraint imposed on the cross-anisotropy by ÂA 2 , as described in Sect. 6. The formal structure of ½Q À1 given in (40) with only a few independent C i poses a strong limitation on the congruence relation. The congruence requires ½Q À1 to be a nonsingular matrix only. However, identical zero blocks in ½Q À1 from (40) and in ½C iso provide a major advantage for the determination of C i , namely, the search for the 9 Â 9 coupling matrix ½Q À1 can be split into two independent and smaller tasks: (1) Coupling of the upper left 3 Â 3 blocks (2) Coupling of the lower right 6 Â 6 blocks.
The solution of the second task can be taken as the lower right 6 Â 6 block of ½Q À1 from (46). Unfortunately, the first task is less trivial. The upper left 3 Â 3 block of ½C iso from (43) should be coupled with the upper left 3 Â 3 block of ½C ÂA3 from (44) using just the upper left 3 Â 3 block of ½Q À1 independently of the remaining components. Obeying the structure of ½Q À1 from (40), the first task takes the form Using the powerful command Reduce[] from Mathematica, one can algebraically reduce the system. This reduction leads to the constraint, x 2 ¼ a 1=c , imposed on a; b and c, identical as in Â A 2 described in Sect. 6. Hence, the construction of the inverse anisotropy tensor Q À1 for ÂA 3 without constraints, i.e., preserving all pure anisotropy parameters, a; b and c, is not possible. If the elastic constant m was allowed 15 to enter Q, then E ÂA3 given in (22) could be decomposed could be interpreted 16 . Tensor E v I describes the isotropic elastic stiffness for the special case with m ¼ 0 and E ¼ E v .

Calibration of pure cross-anisotropy
Two methods of calibration of the ÂA constants will be presented: static triaxial tests with small stress cycles applied in different directions and dynamic tests with different wave types propagated in different directions. In both cases, the average stress should be isotropic. Otherwise, the ÂA must be calibrated jointly with the rA, which is much more difficult.
A combined partly dynamic and partly static, cyclic calibration should be avoided because the anisotropy of the small-strain stiffness may change with the size of the amplitude. Strain amplitudes due to wave propagation are usually much smaller than the ones from static cycles.

Static calibration of · A 1
In this section, two methods to determine a; E v and m h for the ÂA 1 are presented. The first one is based on two saturated, undrained triaxial tests, and the second one needs two drained triaxial tests with measurement of the volume change. In isotropic elasticity, the volumetric and deviatoric behavior can be described separately. Isochoric (at constant volume % undrained [22]) stress paths are perpendicular to the hydrostatic axis. In anisotropic elasticity, the inclination may be measured, see Fig. 5. The inclination g is different for the v-sample cut parallel and for the h-sample cut perpendicular to the direction of sedimentation from the same material. This can be illustrated with the results from cyclic stress tests on kaolin [29], see Fig. 6. The inclinations are interrelated by and (51) holds for any ÂA. Hence g v and g h provide equivalent information for the calibration of a and m, for which two conditions are required. In the coordinate 15 No true separation of elasticity and pure anisotropy anymore. 16 The root of a symmetric matrix A can be found from spectral where D is the diagonal matrix with eigenvalues of A and G contains the corresponding orthonormalized eigenvectors in rows. system from Fig. 5 the first condition can be formulated for the v-sample Assuming E v ¼ 1, the right-hand side of (52) 2 is a function of a and m h only and g v is known. The second condition is based on the observation that identical stress amplitudes q ampl cause different strain amplitudes in the v-and hsample. The ratio r ¼ e ampl av =e ampl ah 6 ¼ 1 can be measured in the undrained test. Again, in the coordinate system from Fig. 5, the second condition can be expressed by three equations wherein _ e v and _ e h are strain rates in v-sample and h-sample caused by the same stress rate _ In the conventional undrained triaxial tests with _ r tot r ¼ 0, one may express these strain rates as wherein the effective stress rates and the rates of pore pressures _ u v 6 ¼ _ u h may be different in v-and h-samples (in spite of the same _ q). Using the conditions, one may express a and m h by analytical formulas, see Appendix A.
With a and m h in hand, one may determine the module E v ¼ s _ r tot av = _ e av . The rates _ r tot av and _ e av should be measured from the undrained v-sample. The scaling factor sðm h ; aÞ can be determined substituting into _ The system _ The complete solution is given in Appendix A.
Alternatively, the ÂA 1 parameter along with the elastic constants can be determined from the conventional drained triaxial tests (at _ r r ¼ 0). From a compression of a v-sample and a h-sample, one obtains E v ¼ _ r av = _ e av and E h ¼ _ r ah = _ e ah , respectively. The measurement of volumetric and axial deformations leads to the following system which can be solved for a; m vh ; m hv ; m hh , see Appendix A.

Dynamic calibration of · A 2
In this section only the dynamic calibration of ÂA 2 is discussed. A static calibration of b via G v is possible, but it needs a hollow-cylinder torsion test on a v-sample. Anisotropic elastic parameters can be determined from the measurements of wave velocities (dynamic tests) in different direction of propagation n. Using this direction, the acoustic tensor can be built wherein E is the stiffness and n is unit vector. The eigenvalues of C ik are related to the velocities of different waves wherein q is the mass density. Three eigenvalues qv 2 may be obtained from det C jk À qv 2 d jk À Á ¼ 0. They may correspond, in general, to three different waves with different velocities, all propagating along n. The corresponding eigenvectors A describe the polarizations of displacement amplitudes. In the case of isotropic elasticity, it is one Pwave with Akn and two S-waves with A ? n, Fig. 7. The velocities v S and v P are independent of n.
In a cross-anisotropic medium with E ÂA2 , the velocities of propagation and the polarization directions depend on the anisotropy parameters, a and b, and on the angle between n and m. The explicit expressions for C ik in the case of any n and m ¼ f0; 0; 1g are given in Appendix B. We examine two directions of propagation, nkm (index v) and n ? m (index h) with m ¼ f0; 0; 1g, Fig. 8.
For such n, the polarization A can be either perpendicular or parallel to n. The respective eigenvalues are denoted as qv 2 Sij and qv 2 Pij , wherein i is the direction of propagation and j is the direction of polarization, both taking the values h or v. The velocities for ÂA 2 can be easily found as the eigenvalues of tensors given in (72) in Appendix B Both parameters, a and b, can be calibrated from vertical and horizontal waves 17 alone, using (60), see Fig. 9.
Four independent wave velocities, v Pvv ; v Phh ; v Shh and v Svh ¼ v Shv , can be measured and (60) can be solved for two pure anisotropic parameters and two elastic parameters, E ¼ E v and m ¼ m h , Determination of all five parameters for stiffness (22) requires additionally a wave velocity in an inclined direction n, say for n Á m ¼ 1= ffiffi ffi 2 p [8,27].

Tests of · A
Recently, Mašín and Rott [19] have reviewed numerous experiments on sedimentary clays. They concluded that, using the nomenclature of (42), most clays need c [ 1=2, which can be covered by ÂA 2 or ÂA 3 but not by ÂA 1 . It is claimed [19] that the average value should be c % 4=5. This observation was based on tests which could be blurred by the rA. However, for practical purposes, such results are sufficient because ÂA has been shown to be v v A S1 S1 v Fig. 7 Direction of propagation n with two shear waves, v S1 and v S2 , and one pressure wave, v P , for isotropic elasticity Fig. 9 Setup of bender elements for the determination of ÂA 2 parameters: a waves with vertical propagation, b waves with horizontal propagation dominant over rA in highly overconsolidated clays [19] as well as in kaolin [9]. Unfortunately, only a few tests from [19] were carried out under hydrostatic stress. In consequence, not much usable data can be found. However, some results from London Clay and Gault Clay referred to in [30] confirmed the discrepancies from c ¼ 1=2 and speak for ÂA 2 rather than for ÂA 1 . The exponent c ¼ 1=2 was estimated for Bangkok Clay under isotropic stress [26]. Measured values of c are presented for different as in Fig. 10.
Some dynamic test data for Kenya Sand [8] and Hostun Sand [27] at different isotropic stress levels, p, revealed an influence of p on the parameter b. This strange effect can be attributed to errors in measurements or to partial destruction of ÂA by isotropic loading. Tests with temporary overloading (up to a high p and back) could help to confirm such a degradation. The dynamic tests prove c ! 1=2 for sands.
Parameter b and the ratio b=c are plotted as functions of a in Figs. 11 and 12, respectively. The ratio b=c ¼ 2 was assumed in ÂA 2 because of the mathematical convenience. Due to the scatter of experimental data, one can neither confirm nor reject this assumption.

Graphic representation of anisotropy
For constitutive rate-type models in the form of an isotropic function _ rðr 0 ; _ eÞ, the well-known concept [11] of response envelopes can be used for the graphic representation of stiffness. The 2D plots of response envelopes to strain disturbances require that the initial stress, r 0 , and all strain rates, _ e, are co-axisymmetric, i.e., axisymmetric with respect to the same symmetry axis.
In the case of ÂA, the sedimentation dyad, M ¼ m m, appears as an additional argument in _ rðr 0 ; _ e; MÞ. This dyad needs not be co-axisymmetric with r 0 and _ e. In such case, the usual 2D response envelopes cannot be plotted, if ÂA spoils the co-axisymmetry of r 0 and _ r. For a general graphic representation of stiffness with any ÂA, the original concept [11] can be extended. In this extension, the stress increments 18 , Dr, need not be co-axisymmetric with r 0 .

2D response envelopes
A response envelope is a polar representation of a tangential stiffness at a given stress r 0 . Starting from a diagonal and axisymmetric initial stress, r 0 ¼ diagðr 0 1 ; r 0 2 ; r 0 3 Þ with r 0 2 ¼ r 0 3 , different axisymmetric strain increments of constant length, are applied, Fig. 13a. The envelope of the corresponding stress increments, Dr ¼ Drð/Þ, is termed the response envelope. Linear elasticity maps a circle (63) in the strain space to an ellipse in the stress space, Fig. 13c. Increments Dr are co-  [30], Gault Clay (GC) [30], Bangkok Clay (BC), [26] Hostun Sand (HS) [27] and Kenya Sand (KS1, KS2) [8]  axisymmetric with r 0 , if r 0 is co-axisymmetric with De and ÂA is absent or its m is parallel to the symmetry axis. In such cases, the end stresses, r 0 þ Dr, can be plotted. These plots are quite common in the geotechnical literature. Usually, they are shown on the Rendulić plane, ffiffi ffi 2 p r r À r a , or on the plane of isometric Roscoe invariants, P À Q.
Generally, r 0 þ Dr cannot be plotted because the ÂA may spoil the co-axisymmetry between Dr and r 0 . However, all Dr are coplanar, if all De are and because the constitutive relation, _ rðr 0 ; _ e; MÞ ¼ Eðr 0 ; MÞ : _ e, is incrementally linear. Let the following orthogonal strain increments: produce stress increments, Dr P and Dr Q , respectively. These two increments span a plane in 6D stress space. All other stress responses lie in this plane due to the linearity of E. In other words, any response is a linear combination of Dr P and Dr Q . After orthonormalization of Dr P and Dr Q , they constitute the orthogonal basis fe H P ; e H Q g on the response plane and we may introduce the coordinates, DP H and DQ H , on this plane. Any stress response can be represented as for example Drð/ P Þ ¼ DP H e H P .

An example of 2D response
Experiments on kaolin [9] show that the effects from ÂA dominate over the ones from rA, Fig. 3. It turns out that, for kaolin, the ÂA 1 with a single anisotropy parameter a simulates the experiments sufficiently well and b is not necessary. In sedimentary clays, however, ÂA 1 can be inaccurate, see Sect. 9. As an example, 2D response envelopes from the superposition of rA from (14) and ÂA 2 are plotted in the DP H À DQ H plane in Fig. 14.

3D response envelopes
To plot 3D response envelopes, solely the coaxiality of r 0 and _ e in _ rðr 0 ; _ eÞ is required. If the ÂA is present, all arguments in _ rðr 0 ; _ e; MÞ must be coaxial. Starting from a given initial stress, r 0 ¼ diagðr 0 1 ; r 0 2 ; r 0 3 Þ, diagonal , axisymmetric strain increments of constant length, Deð/; wÞ ¼ rdiag sin /; cos / cos w; cos / sin w ð Þ with r ¼ const % 0:0001 and 0 /; w\2p; ð65Þ are applied, Fig. 13b. They can be encompassed by a sphere in the 3D space of principal strains. In the case of a linear elastic constitutive relation, _ rðr 0 ; _ eÞ ¼ Eðr 0 Þ : _ e, the end stresses, r 0 þ Dr, form an ellipsoidal response envelope in the 3D space of principal stresses, Fig. 13d. The respective stress increments, Dr ¼ Drð/; wÞ, are coaxial with r 0 , if r 0 and De are. Generally, the coaxiality of r 0 and Dr may be violated by the presence of the ÂA, when M is not coaxial with r 0 .
Similarly as in the 2D case, we define three orthogonal strain increments: They correspond to the following angles: The respective stress increments, Dr P ; Dr Q and Dr R , are not necessarily orthogonal, but they span a 3D subspace of the 6D stress space. Analogously as in the 2D case, these stress increments can be orthonormalized to define the basis fe H P ; e H Q ; e H R g and the coordinate system DP H À DQ H À DR H of this subspace. Due to the incremental linearity, all stress increments can be expressed as linear combinations of the basis tensors, ffiffi ffi 3 p Þ and w P ¼ p=4.

An example of 3D response
The 3D stress response envelopes were obtained with the identical constitutive model and the same material constants as for the 2D ones from Fig. 14. The 3D strain increments De were applied to plot Dr in DP H À DQ H À DR H system, Fig. 15.  (14) and with ÂA 2 : 2D isometric stress plots (b,c,d) were calculated at different diagonal initial stresses r 0 and for the same sedimentation m ¼ 1; 2; 3 f g !

Polar diagrams of wave velocity
Using the acoustic tensor C from (58), the velocities v of different waves can be plotted as functions of the direction of propagation n. The directional dependence of wave velocities can be then visualized in the form of polar diagrams for each wave type. An example of polar diagrams obtained with the superposition of ÂA 2 and rA from (13) is shown in Fig. 16.

Scaling of yield functions
The anisotropy tensor Q from ÂA 1 and ÂA 2 may have a variety of applications beyond elasticity. A yield stress criterion describes the boundary of all accessible stress states, FðrÞ 0, where FðrÞ is an isotropic function of stress. For example, Matsuoka and Nakai [18] proposed the following yield function FðrÞ trrtr r À1 À Á À 8 tan 2 u À 9; wherein u is the friction angle.
The ÂA can be imposed to stress using the anisotropy tensor from (32) and substituted into FðrÞ, i.e., F ÂA2 ðr ab Þ ¼ FðQ abcd r cd Þ. As an example, FðrÞ from (67) with the ÂA 2 was plotted in the deviatoric plane, Fig. 17.
The transformed yield function F ÂA2 ðrÞ requires calibration of the corresponding friction angle u ÂA2 .
In the literature, one may find some attempts to make a yield surface FðrÞ cross-anisotropic, e.g., [16]. In comparison, scaling with the anisotropy tensor, Q, is an elegant and easy method.

Summary
Inherent cross-anisotropy and stress-induced anisotropy can be easily superposed within the elastic range, in particular dealing with geotechnical (barotropic) elastic potentials. The pure anisotropy tensor, Q, depends on the sedimentation direction, m, and some material constants. The simplified versions, ÂA 1 and ÂA 2 , of cross-anisotropy could be used to build such Q but not the general form,  (14) and ÂA 2 : 3D isometric stress plots (b,c,d) were calculated at different diagonal initial stresses r 0 and for the same sedimentation m ¼ 1; 2; 3 f g ! Second Law, if superposed with hyperelasticity. The pure anisotropy can be applied also to any isotropic potential function, for example to a yield surface. The proposed calibration procedure for Q can be based on static, cyclic or dynamic tests. The popular concept of response envelopes [11] has been extended to provide the graphic representation of polar stiffness at presence of ÂA. For this purpose, a new isometric representation system has been proposed. The concept of pure anisotropy has been compared to some recent approaches from the literature. Visualization of the superposed ÂA 2 and rA conducted with the algebra program Mathematica has been given in examples. All notebooks and packages involved in this paper are available from the authors. and set of equations (60) can be determined from the eigenvalues of C.
Let us define three of polarization cosines P i ¼ n Á Ã i . In the case of isotropic elasticity, P ¼ f1; 0; 0g means one P-and two S-waves. At presence of ÂA 2 , one can speak of only one S-wave 19 . Its polarization is perpendicular to both n and m. Two other waves lie in the plane spanned by n and m. All three wave velocities are different. For example, a ¼ 1:8 and b ¼ 1:2 in ÂA 2 with n ¼ f1; 2; 3g ! yield P ¼ f0:94; 0; 0:33g, wherein the second polarization corresponds to the S-wave. The other two polarizations depend on a; b and on the angle between n and m.

Author Contributions Not applicable
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest Not applicable
Availability of data and material All data were digitalized from the figures published in the literature.
Code availability The relevant packages and notebooks for the algebra program Mathematica are available from the authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.