# Errors Solving Elliptic PDES with FEM

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{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} 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{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} 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["NDSolveFEM"];

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}=NDSolveProcessEquations[
{govEqns=={0,0}, fixedBC, prescribedDBC},
{ux,uy},
{x, y} ∈ tMesh,
Method->"FiniteElement"]

NDSolveIterate[tState]

NDSolveProcessSolutions[tState]

uf = ux /. NDSolveProcessSolutions[tState]
vf = uy /. NDSolveProcessSolutions[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: NDSolveProcessEquations::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: NDSolveProcessEquations::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.

• Boundary conditions need to be applied at x^2 + y^2 ==4.5^2 and x^2 + y^2 == 3.5^2. – bbgodfrey Feb 29 '16 at 1:13
• Also, the alternative (second) representation of displacement2StrainOp seems preferable, although you may not need to use Inactive. It also may not be necessary to use the components of NDSolve instead of NDSolve itself. If so, this could save you much work. – bbgodfrey Feb 29 '16 at 1:48
• @bbgodfrey you're right. 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. Sorry for the delay; timezone differences and the comment box ate the backslash before Element. – seyisulu Feb 29 '16 at 6:37
• Are you saying, then, that a zero-traction condition is to be applied everywhere on the interior and exterior boundaries except at the points {0, 4.5} and {0, -4.5}, where values are specified for ux and uy? – bbgodfrey Feb 29 '16 at 13:58
• @bbgodfrey yes! But I do not know how to combine Neumann and Dirichlet boundary conditions. – seyisulu Feb 29 '16 at 14:27

Here is how I'd do it. First, let's write a function that generates the plane stress PDE:

<< NDSolveFEM
ClearAll[PlaneStress];
PlaneStress[{Y_, nu_}, {u_, v_}, X : {x_, y_}] :=
Module[{pStress},
pStress = -Y/(1 -
nu^2)*{{{{1, 0}, {0, (1 - nu)/2}}, {{0, nu}, {(1 - nu)/2,
0}}}, {{{0, (1 - nu)/2}, {nu, 0}}, {{(1 - nu)/2, 0}, {0, 1}}}};

{Inactive[Div][pStress[[1, 1]].Inactive[Grad][u, X], X] +
Inactive[Div][pStress[[1, 2]].Inactive[Grad][v, X], X],
Inactive[Div][pStress[[2, 1]].Inactive[Grad][u, X], X] +
Inactive[Div][pStress[[2, 2]].Inactive[Grad][v, X], X]}]

Now, you'd be able to call:

tE = 1.69 10^9;
tv = 0.31;
eqn = PlaneStress[{tE, tv}, {ux[x, y], uy[x, y]}, {x, y}] == {0, 0};

With the region and the boundary conditions:

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}}];
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];
{uxfun,uyfun}=NDSolveValue[{eqn,fixedBC,prescribedDBC},{ux,uy},{x,y} ∈ tRegion];

Visualizing:

mesh = uxfun["ElementMesh"];
Show[{
mesh["Wireframe"[ "MeshElement" -> "BoundaryElements"]],
ElementMeshDeformation[mesh, {uxfun, uyfun} ][
"Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]},
ImageSize -> Large]

This looks reasonable. Adding a material dependency is then doing something like:

eqn = PlaneStress[{tE*x, tv}, {ux[x, y], uy[x, y]}, {x, y}] == {0, 0}

I just chose tE*x but you can go wild here, as long as it's linear it's fine.

Update:

Derivation of the PlaneStress operator. We start from the stress strain relation and replace the derivatives ex, ey and gxy with the actual definitions to give sx, sy and txy:

{sx, sy, txy} = Y/(1 - \[Nu]^2) ( {{1, \[Nu], 0},{\[Nu], 1, 0},{0, 0, (1 - \[Nu])/2}}).({{ex},{ey},{gxy}}) /. {ex -> D[u[x, y], x], ey -> D[v[x, y],y],gxy -> (D[u[x, y], y] + D[v[x, y], x])}

(* {{(Y*(\[Nu]*Derivative[0, 1][v][x, y] + Derivative[1, 0][u][x, y]))/(1 - \[Nu]^2)},
{(Y*(Derivative[0, 1][v][x, y] + \[Nu]*Derivative[1, 0][u][x, y]))/(1 - \[Nu]^2)},
{(Y*(1 - \[Nu])*(Derivative[0, 1][u][x, y] + Derivative[1, 0][v][x, y]))/(2*(1 - \[Nu]^2))}} *)

Now, take the derivatives of sx, sy and txy:

{{-D[sx, x] - D[txy, y]},{-D[txy, x] - D[sy, y]}} // MatrixForm

This will give you the plane stress.

• Nice analysis, which looks quite credible (+1). However, I cannot reproduce the equations given in LaTex format at the top of the question from your Activate[eqn]. Is the question in error, or have I missed something? By the way, thanks for your earlier comment. – bbgodfrey Feb 29 '16 at 21:16
• You may wish to add << NDSolveFEM (backquotes do not reproduce well here) to the beginning of your answer. – bbgodfrey Feb 29 '16 at 21:21
• @bbgodfrey, thanks for the for [\Element] fix and yes, a <<NDSolveFEM is needed. I'll look at your other question later, I need to attend to something else first. – user21 Feb 29 '16 at 22:34
• @bbgodfrey I got similar results by simplifying both mine and user21's with the only difference being a multiplication by -1; also, I don't mind applying the Dirichlet conditions over an angle of say 5 degrees if you can show me how. user21, pardon my lack of mathematical prowess, but is that a fourth order tensor you used? Thank you guys; I wish I could buy you all beers! – seyisulu Mar 1 '16 at 5:48
• @user21 though I've established that in the end the results are similar, I would very much appreciate an explanation of how you arrived at the construct you have so I could apply the same to other equations I might have. Also, does using NDSolveValue instead of NDSolve have any effect? Thanks in advance. – seyisulu Mar 1 '16 at 7:06