anisotropicBiotElastic
This law calculates small-strain orthotropic effective stress. Pore pressure can be included by nesting it inside poroMechanicalLaw. The runtime type is:
anisotropicBiotElastic
User Guide
What it computes
The cell-centred correct obtains the total strain from the displacement gradient:
epsilon = symm(grad(D)) // total formulation
epsilon = epsilon.oldTime() + symm(grad(DD)) // incremental formulation
In the branch labelled as 2-D, the implemented stress relation is:
sigma_xx = A11*epsilon_xx + A12*epsilon_yy
sigma_yy = A21*epsilon_xx + A22*epsilon_yy
sigma_xy = A44*epsilon_xy
The other stress components are not assigned in this branch. In the 3-D branch, the relation is:
sigma_xx = A11*epsilon_xx + A12*epsilon_yy + A31*epsilon_zz
sigma_yy = A12*epsilon_xx + A22*epsilon_yy + A23*epsilon_zz
sigma_zz = A31*epsilon_xx + A23*epsilon_yy + A33*epsilon_zz
sigma_xy = A44*epsilon_xy
sigma_yz = A55*epsilon_yz
sigma_xz = A66*epsilon_xz
The same relations are applied to the boundary values. The law implements only correct(volSymmTensorField&). Its surface overload aborts with notImplemented. The inherited point overload also aborts. The inherited CompactListList overload is unavailable on foam-extend and aborts with notImplemented on the other forks.
Model options
The elastic properties are read as plain scalars, so the dictionary does not attach dimension sets to them. Use one consistent pressure unit for all Young's and shear moduli.
| Entry | Required | Description |
|---|---|---|
Ex | yes | Young's modulus in the x direction |
Ey | yes | Young's modulus in the y direction |
Ez | in the 3-D branch | Young's modulus in the z direction |
nuxy | yes | Poisson's ratio nuxy |
nuyz | in the 3-D branch | Poisson's ratio nuyz |
nuzx | in the 3-D branch | Poisson's ratio nuzx |
Gxy | yes | Shear modulus in the xy plane |
Gyz | in the 3-D branch | Shear modulus in the yz plane |
Gzx | in the 3-D branch | Shear modulus in the zx plane |
rho | for direct use | Density, [1 -3 0 0 0 0 0] |
regionName | no | Base mesh region; otherwise solid or region0 |
solvePressureEqn | no | Accepted switch; default no, unused here |
pressureSmoothingScaleFactor | no | Default 100, unused here |
Which entries are required depends on a branch flag that is currently inverted, so read the following carefully before setting up a case.
The law selects its 2-D branch when mesh.solutionD()[vector::Z] > 0, whichis true when z is a solved direction — that is, for a 3-D mesh. Thesense of the test is the opposite of the one intended, so today a 3-D casetakes the 2-D branch and reads only Ex, Ey, nuxy and Gxy, silentlyignoring Ez, nuyz, nuzx, Gyz and Gzx and leaving the out-of-planestress components at zero. A genuinely 2-D case, with z empty, takes the 3-Dbranch and aborts with keyword Ez is undefined unless all nine constants aresupplied. This is tracked asissue #334; supply allnine constants until it is fixed.
rho is read through the base-class density functions. When the law is nested inside poroMechanicalLaw, density belongs to the outer law, as in the setup below. The pressure-equation entries are stored by mechanicalLaw, but this law never calls updateSigmaHyd(), so they do not change its stress update.
Recommended dictionary setup
The following seabed-soil values follow the available poroelastic tutorial:
planeStress no;
mechanical
(
seabedSoil
{
type poroMechanicalLaw;
rho rho [1 -3 0 0 0 0 0] 2650;
biotCoeff biotCoeff [0 0 0 0 0 0 0] 1.0;
effectiveStressMechanicalLaw
{
type anisotropicBiotElastic;
Ex 1.2e7;
Ey 1.2e7;
Ez 2e7;
nuxy 0.2;
nuyz 0.24;
nuzx 0.4;
Gxy 0.5e7;
Gyz 1.2e7;
Gzx 1.2e7;
// Optional for a non-standard base mesh name
// regionName solid;
}
}
);
Field glossary
epsilon: cell-centred total small-strain tensor. It is created withNO_READandNO_WRITEand accumulated from its old-time value for an incremental solid model.grad(D),grad(DD): displacement-gradient fields looked up for total and incremental solid models respectively; the law does not own them.sigma: caller-owned cell-centred stress updated bycorrect.impK: temporary uniform pressure-dimensioned field returned to the solid model.
Developer Notes
Class role
anisotropicBiotElastic derives directly from mechanicalLaw. It stores the branch flag model2d_, ten scalar stiffness coefficients and the cell-centred epsilon_ field. It is registered in the linGeomMechLaw selection table under anisotropicBiotElastic.
The source uses #ifdef OPENFOAM_NOT_EXTEND to select writable internal and boundary field accessors. The class itself has no FOAMEXTEND-only interface guards. Its source is listed in both Make/files.openfoam and Make/files.foamextend.
Construction
The base constructor first reads solvePressureEqn with default false and pressureSmoothingScaleFactor with default 100.0. It then reads regionName, or selects an existing solid or region0 mesh. If none is available, it raises a fatal error asking for regionName.
The derived constructor sets model2d_ from whether the z solution direction is greater than zero, zeroes every stiffness coefficient and creates epsilon_. In the 2-D branch it raises a fatal error if x or y is an empty direction. It then reads Ex, Ey, nuxy and Gxy, and derives:
nuyx = nuxy*Ey/Ex
J = 1/(1 - nuxy*nuyx)
A11 = J*Ex
A22 = J*Ey
A12 = J*nuyx*Ex
A21 = J*nuxy*Ey
A44 = 2*Gxy
In the 3-D branch it reads all nine elastic entries and derives:
nuyx = nuxy*Ey/Ex
nuxz = nuzx*Ex/Ez
nuzy = nuyz*Ez/Ey
J = (1 - nuxy*nuyx - nuyz*nuzy - nuzx*nuxz
- 2*nuyx*nuzy*nuxz)/(Ex*Ey*Ez)
A11 = (1 - nuyz*nuzy)/(J*Ey*Ez)
A22 = (1 - nuxz*nuzx)/(J*Ex*Ez)
A33 = (1 - nuyx*nuxy)/(J*Ey*Ex)
A12 = (nuxy + nuzy*nuxz)/(J*Ex*Ez)
A31 = (nuzx + nuyx*nuzy)/(J*Ey*Ez)
A23 = (nuyz + nuyx*nuxz)/(J*Ex*Ey)
A44 = 2*Gxy
A55 = 2*Gyz
A66 = 2*Gzx
The constructor prints the supplied properties and derived coefficients. It does not validate moduli, Poisson's ratios or the denominators, and emits no law-specific warnings.
Key methods
impK()returns a uniform field containingmax(A11, max(A22, A33)). In the 2-D branchA33remains zero. This is the diffusivity of the solid model's implicit Laplacian term and affects the outer-iteration convergence rate rather than the converged answer.correct(volSymmTensorField&)updatesepsilon_fromgrad(D)orgrad(DD), then applies the component-wise constitutive relation to cell and boundary values.correct(surfaceSymmTensorField&)always aborts withnotImplemented.materialTangent()is not overridden. The inherited implementation aborts withnotImplemented, so this law does not supply a Newton material tangent.
Extension points
A related orthotropic law can copy this class, replace the property reads and stiffness derivation, and update the component equations in correct. It must be registered in the appropriate mechanical-law selection table and added to both build lists if all forks are supported. Implement the surface, point and quadrature overloads, and materialTangent(), when the target solid models require those interfaces.
The source is at anisotropicBiotElastic.C.
Tutorials
solids/poroelasticity/rodAndSeabed