Errors Solving Elliptic PDES with FEM

Question

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.

Answer 1

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

<< NDSolve`FEM`
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.