I am trying to solve the equation below governing transversely isotropic plane strain in cartesian coordinates with the given boundary conditions based on code found here using Mathematica 10.1 on OSX El Capitan.
$$ \begin{equation} \begin{aligned} \begin{bmatrix}\varepsilon_x\\\varepsilon_y\\\gamma_{xy}\end{bmatrix}&= \begin{bmatrix} \frac{\partial}{\partial x} & 0\\ 0 & \frac{\partial}{\partial y}\\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y}\\ \end{bmatrix} \begin{bmatrix}u_x (x,y)\\u_y (x,y)\end{bmatrix}\\ \begin{bmatrix}\sigma_x\\\sigma_y\\\tau_{xy}\end{bmatrix}&= \frac{1}{\left(1-v^2\right)} \begin{bmatrix} E&E v&0\\E v&E&0\\0&0&E\left(1-v\right) \end{bmatrix} \begin{bmatrix}\varepsilon_x\\\varepsilon_y\\\gamma_{xy}\end{bmatrix}\\ \begin{bmatrix} \frac{\partial}{\partial x}&0&\frac{\partial}{\partial y}\\ 0&\frac{\partial}{\partial y}&\frac{\partial}{\partial x} \end{bmatrix} \begin{bmatrix}\sigma_x\\\sigma_y\\\tau_{xy}\end{bmatrix}&= \begin{bmatrix}0\\0\end{bmatrix} \end{aligned} \end{equation} $$ Putting it all together to obtain the 2D Navier's equation (displacement formulation): $$ \frac{1}{\left(1-v^2\right)} \begin{bmatrix} \frac{\partial}{\partial x}&0&\frac{\partial}{\partial y}\\ 0&\frac{\partial}{\partial y}&\frac{\partial}{\partial x} \end{bmatrix} \begin{bmatrix} E&E v&0\\E v&E&0\\0&0&E\left(1-v\right) \end{bmatrix} \begin{bmatrix} \frac{\partial}{\partial x} & 0\\ 0 & \frac{\partial}{\partial y}\\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y}\\ \end{bmatrix} \begin{bmatrix}u_x (x,y)\\u_y (x,y)\end{bmatrix} = \begin{bmatrix}0\\0\end{bmatrix} $$ Subject to the boundary conditions: $$ \begin{equation} \begin{aligned} u_x(0,-4.5)&=0.0\\ u_y(0,-4.5)&=0.0\\ u_x(0,4.5)&=0.0\\ u_y(0,4.5)&=0.2 \end{aligned} \end{equation} $$ In the region $x^2+y^2=3.5^2$ subtracted from $x^2+y^2=4.5^2$.
I have written the following code to solve it but get a NDSolve::dgsvars
error.
σ = {σx[x, y], σy[x, y], τxy[x, y]};
ε = {εx[x, y], εy[x, y], γxy[x, y]};
u={ux[x,y],uy[x,y]};
displacement2StrainOp[disp_?VectorQ]:=
Block[{x,y,null},
{Inactive[Div][{disp[[1]]},{x}],
Inactive[Div][{disp[[2]]},{y}],
Inactive[Div][{disp[[2]],disp[[1]]},{x,y}] 1/2}];
strain2StressOp[strain_?VectorQ,Ep_,v_]:=Block[
{multiplier=(1-v^2)^-1,
matrix={{Ep,Ep v,0},{Ep v,Ep,0},{0,0,Ep(1-v)}}},
multiplier matrix.strain];
stress2GovOp[stress_?VectorQ]:=Block[{x,y},
{Inactive[Div][{stress[[1]],stress[[3]]},{x,y}],
Inactive[Div][{stress[[3]],stress[[2]]},{x,y}]}];
govEqns=
stress2GovOp[
strain2StressOp[
displacement2StrainOp[u],1.69 10^9,0.31]]//FullSimplify;
Needs["NDSolve`FEM`"];
tRegion=ImplicitRegion[
((x - 0)^2 + (y - 0)^2 <= 4.5^2) &&
((x - 0)^2 + (y - 0)^2 >= 3.5^2),
{{x, -5, 5}, {y, -5, 5}}];
tMesh=
ToElementMesh[tRegion, MaxCellMeasure->.25, "MeshOrder"->2];
fixedBC=
DirichletCondition[
{ux[x,y]==0., uy[x,y]==0.}, x==0 && y==-4.5];
prescribedDBC=
DirichletCondition[
{ux[x,y]==0., uy[x,y]==-0.2}, x==0 && y==4.5];
tE=1.69 10^9;
tv=0.31;
Show[
RegionPlot[tRegion, AspectRatio -> Automatic],
tMesh["Wireframe"["ElementMeshDirective"->
Directive[Thick,EdgeForm[LightRed],FaceForm[Opacity[0.3],LightBlue]]]],
ListPlot[tMesh["Coordinates"],
{Axes->True, AspectRatio->1,
PlotStyle->Directive[Red, PointSize[0.01]]}],
ImageSize->Large]
{tState}=NDSolve`ProcessEquations[
{govEqns=={0,0}, fixedBC, prescribedDBC},
{ux,uy},
{x, y} ∈ tMesh,
Method->"FiniteElement"]
NDSolve`Iterate[tState]
NDSolve`ProcessSolutions[tState]
uf = ux /. NDSolve`ProcessSolutions[tState]
vf = uy /. NDSolve`ProcessSolutions[tState]
dat = {Table[x, {x, {-4.5, 0, 4.5, 0}}],
Table[y, {y, {0, -4.5, 0, 4.5}}]} // Transpose
inf = Map[{uf[#[[1]], #[[2]]], vf[#[[1]], #[[2]]]} &, dat]
Show[Graphics[{Directive[LightRed, Thick], Circle[{0, 0}, 4.5]}],
ListPlot[{dat, dat + inf}, AspectRatio -> 1]]
tProbMesh = uf["ElementMesh"]
tDeformedProbMesh = ElementMeshDeformation[uMesh, {uf, vf}]
tDeformedProbMeshx5 =
ElementMeshDeformation[uMesh, {uf, vf}, "ScalingFactor" -> 5]
Show[
tProbMesh[
"Wireframe"[
"ElementMeshDirective" ->
Directive[Thin, EdgeForm[LightGreen], FaceForm[]]]],
tDeformedProbMeshx5[
"Wireframe"[
"ElementMeshDirective" ->
Directive[Thick, EdgeForm[LightRed],
FaceForm[Opacity[0.3], LightBlue]]]]]
But then I get the error message NDSolve::dgsvars: "The differentiation variables {x} given for Inactive[Div] should be the spatial independent variables {x,y}."
and the equation is not solved.
Modifying the code with changes to the functions below and reevaluating gives a different error message: NDSolve`ProcessEquations::femper: -- Message text not found -- (Div[{6.45038*10^8 Div$7030,5.796*10^8 ux$7029+1.86968*10^9 uy$7028}])
.
displacement2StrainOp[disp_?VectorQ]:=
Block[{x,y,null},
{Inactive[D][disp[[1]],x],
Inactive[D][disp[[2]],y],
Inactive[Div][{disp[[2]],disp[[1]]},{x,y}] 1/2}];
I do not know what else to try and will appreciate help on this issue. Am I missing something fundamental? I am open to other approaches to solve the problem as long as it uses FEM.
The reason for defining the problem this way is that I plan to later introduce a Young's modulus functionally dependent on x and y so I need a way to differentiate whatever expression passed in.
Thank you.
Edit 1:
As per the comment below (thanks), there is a zero traction condition at x^2 + y^2 == 3.5^2
. I added it with
zeroTractionBC=NeumannValue[0,x^2 + y^2 == 3.5^2];
and then in the solver:
{tState}=NDSolveProcessEquations[{
govEqns=={ zeroTractionBC, zeroTractionBC }, fixedBC, prescribedDBC},
{ux,uy},
{x,y}\[Element]tMesh, Method->"FiniteElement"]
but still have the same error:
NDSolve`ProcessEquations::temper:
-- Message text not found -- (Div[{6.45038*10^8 Div$17301,5.796*10^8 ux$17300+1.86968*10^9 uy$17299}])
Edit 2: Based on issues from the post here I have tried the same code on Mathematica 10.3 running on Windows 10 but I get the same set of errors. Please I will appreciate any help at this point.
x^2 + y^2 ==4.5^2
andx^2 + y^2 == 3.5^2
. $\endgroup$Inactive
. It also may not be necessary to use the components ofNDSolve
instead ofNDSolve
itself. If so, this could save you much work. $\endgroup$x^2 + y^2 == 3.5^2
. I added it withzeroTractionBC=NeumannValue[0,x^2 + y^2 == 3.5^2];
and then in the solver:{tState}=NDSolve
ProcessEquations[{govEqns=={zeroTractionBC,zeroTractionBC}, fixedBC, prescribedDBC},{ux,uy}, {x,y}[Element]tMesh, Method->"FiniteElement"]` but still have the same error. Sorry for the delay; timezone differences and the comment box ate the backslash beforeElement
. $\endgroup${0, 4.5}
and{0, -4.5}
, where values are specified forux
anduy
? $\endgroup$