Bug introduced in 13.1 and fixed in 13.2.0
I am trying to compute a deformation and a stress of a material with a custom material law. I have successfully computed the deformation using SolidMechanicsPDEComponent and NDSolveValue. Now I would like to compute and visualize the stress. I have to compute the strain using SolidMechanicsStrain first, which works fine. However the SolidMechanicsStress then sometimes returns, that the stress is zero, which is wrong.
Here is a part of my code:
strain = SolidMechanicsStrain[vars, pars, displacement];
cauchy = SolidMechanicsStress[vars, Join[pars, <|"OutputStressMeasure" -> "Cauchy"|>], strain, displacement];
firstPK = SolidMechanicsStress[vars, Join[pars, <|"OutputStressMeasure" -> "FirstPiolaKirchhoff"|>], strain, displacement];
vars are variables, pars are parameters, displacement is a displacement computed using NDSolveValue.
Here is a the visualization of one of the components of the Cauchy stress tensor. All other components are also equal to zero. The first Piola-Kirchhoff is also equal to zero.
ContourPlot[cauchy[[1, 1]], {x, y} ∈ deformedMesh, PlotRange -> All, PlotLegends -> Automatic]
Here are other parts of the code:
ClearAll["Global`*"]
Needs["NDSolve`FEM`"]
coords = {x, y};
deformation = {u[x, y], v[x, y]};
Young = 10^9;
ν = 1/3;
n = 1.5;
α = Pi/4;
Q = {{Cos[α], -Sin[α]}, {Sin[α], Cos[α]}};
GetK[Young_, ν_] := Young/(3*(1 - 2*ν));
GetM[Young_, ν_] := Young*(1 - ν)/((1 + ν)*(1 - 2* ν));
GetG[Young_, ν_] := Young/(2*(1 + ν));
K = GetK[Young, ν];
M = GetM[Young, ν];
G = GetG[Young, ν];
rectangle = Rectangle[{0, 0}, {1, 1}];
mesh = ToElementMesh[rectangle];
vars = {deformation, coords};
pars = <|"MaterialModelFunction" -> QRElasticity, "ModelForm" -> "PlaneStrain", "Thickness" -> 1, "ConstitutiveStressMeasure" -> "FirstPiolaKirchhoff", "MassDensity" -> 980, "BulkModulus" -> K, "PWaveModulus" -> M, "ShearModulus" -> G, "n" -> n, "Q" -> Q|>;
pdeQRElasticity = SolidMechanicsPDEComponent[vars, pars];
pde = {pdeQRElasticity == SolidBoundaryLoadValue[x == 1, vars, pars, <|"Pressure" -> {p, 0}|>], DirichletCondition[{u[x, y] == 0}, x == 0], DirichletCondition[{v[x, y] == 0}, x == 0]};
AbsoluteTiming[displacement = NDSolveValue[pde /. p -> 300000000, {u, v}, {x, y} \[Element] mesh];]
deformedMesh = ElementMeshDeformation[mesh, displacement, "ScalingFactor" -> 1];
The QRElasticity is the custom material law for fiber-reinforced materials. The parameter n describes, how much more stiffer are fibres than the material. The matrix Q is a matrix of rotation that describes the orientation of fibres in the material.
QRElasticity[vars_, pars_, data_] := Module[{u, x, dim, idm, n, Q, K, M, G, F, FAni, R, sinθ, cosθ, a, b, γ, U, RU, δ, ε, pi, σ, τ, STilde11, STilde22, STilde12, STilde, UInvTrans, stressMatrix},
u = vars[[1]];
x = vars[[-1]];
K = pars["BulkModulus"];
M = pars["PWaveModulus"];
G = pars["ShearModulus"];
n = pars["n"];
Q = pars["Q"];
dim = Length[u];
idm = IdentityMatrix[dim];
(*Print["K = ",K];
Print["M = ",M];
Print["G = ",G];*)
F = ConstantArray[0, {dim, dim}];
F[[1 ;; dim, 1 ;; dim]] = idm + Grad[u, x];
FAni = Q . F . Inverse[Q];
(*Print["F =",F//MatrixForm];
Print["FAni =",FAni//MatrixForm];*)
R = ConstantArray[0, {dim, dim}];
sinθ = -FAni[[2, 1]]/Sqrt[FAni[[1, 1]]^2 + FAni[[2, 1]]^2];
cosθ = FAni[[1, 1]]/Sqrt[FAni[[1, 1]]^2 + FAni[[2, 1]]^2];
R = {{cosθ, sinθ}, {-sinθ, cosθ}};
(*Print["R = ",R//MatrixForm];*)
a = Sqrt[FAni[[1, 1]]^2 + FAni[[2, 1]]^2];
b = (FAni[[1, 1]] FAni[[2, 2]] - FAni[[1, 2]] FAni[[2, 1]])/Sqrt[FAni[[1, 1]]^2 + FAni[[2, 1]]^2];
γ = (FAni[[1, 1]] FAni[[1, 2]] + FAni[[2, 1]] FAni[[2, 2]])/(FAni[[1, 1]]^2 + FAni[[2, 1]]^2);
U = {{a, a*γ}, {0, b}};
(*Print["U = ",U//MatrixForm];*)
(*Print["RU = ",Simplify[R.U]//MatrixForm];*)
δ = Log[Sqrt[a^n*b^(1/n)]];
ε = Log[Sqrt[(a^n)/(b^(1/n))]];
(*Print["delta = ",δ];
Print["epsilon = ",ε];*)
pi = 4*K*δ;
σ = 2*M*ε;
τ = G*γ;
STilde11 = 1/2*(n*pi + n*σ);
STilde22 = 1/2*(pi/n - σ/n);
STilde12 = (b/a)*τ;
STilde = {{STilde11, STilde12}, {STilde12, STilde22}};
(*Print["STilde = ",STilde//MatrixForm];*)
UInvTrans = Inverse[Transpose[U]];
(*Print["UInvTrans = ", UInvTrans//MatrixForm];*)
stressMatrix = Inverse[Q] . R . STilde . UInvTrans . Q;
stressMatrix = Simplify[stressMatrix[[1 ;; dim, 1 ;; dim]]];
(*Print["StressMatrix = ", stressMatrix//MatrixForm];*)
stressMatrix
]
QRElasticity? People will need the full code to experiment with this. $\endgroup$