1
$\begingroup$

I'm trying to solve the following linear elastic problem:

The pairs of numbers beside each node is the coordinates in meters. The loads are in Newtons.

The part of codes that probably have no mistakes:

ClearAll["Global`*"]

(* Load the FEM package *)
<< NDSolve`FEM`

(* Define the mesh *)
node = {
    {0., 0.}, {10., 2.}, {20., 0.},
    {2., 6.}, {5., 8.}, {17., 6.}
  };
element = {
    {1, 2, 5, 4},
    {2, 3, 6, 5}
  };
mesh = ToElementMesh[
    "Coordinates" -> node,
    "MeshElements" -> {QuadElement[element]}
  ];

(* Constants *)
e = 2.*^11 (* Pa, Young's modulus *);
ν = 0.3    (* dimensionless, Poisson's ratio *);
ρ = 7.8*^3 (* kg m^-3, density *);
h = 0.1    (* m, thickness *);
g = 9.8    (* m s^-2, gravitational acceleration *);

(* Loads *)
f = {
    {0.,     0.   },
    {0.,     0.   },
    {1.*^4,  0.   },
    {0.,    -1.*^4},
    {0.,    -1.*^4},
    {0.,     0.   }
  } (* N, nodal load *);
b = -ρ g (* N m^-3, body force *);

(* Strain (2x2) from displacement (1x2) *)
ε[u_] := Grad[u[x, y], {x, y}];

(* Strain from stress (2x2) *)
σ2ε[σ_] := 1/e ((1 + ν) σ - ν Total@Diagonal[σ] IdentityMatrix[2]);

(* Stress from displacement *)
σ[u_] := e/(1 - ν^2) ((1 - ν) ε[u] + ν Total@Diagonal[ε[u]] IdentityMatrix[2]);

The last three definitions are:

$\boldsymbol{\varepsilon} = \nabla \boldsymbol{u}$

$\boldsymbol{\varepsilon} = \displaystyle\frac{1}{E}\big[(1+\nu)\boldsymbol{\sigma} - \nu(\sigma_{xx} + \sigma_{yy})\mathbf{I}\big]$

$\boldsymbol{\sigma} = \displaystyle\frac{E}{1-\nu^2}\big[(1-\nu)\boldsymbol{\varepsilon} + \nu(\varepsilon_{xx} + \varepsilon_{yy})\mathbf{I}\big]$

The mesh, for future reference:

The part where I'm not so sure about:

The PDE operator, which is basically $\nabla\cdot\boldsymbol{\sigma}$ :

(* PDE operator *)
op = Div[σ[{ux[x, y], uy[x, y]}], {x, y}];
opx = Flatten[op][[1]] (* x component *);
opy = Flatten[op][[2]] (* y component *);

Neumann boundary condition, converting the nodal loads to distributed loads:

nbc := Module[
      {
        line = Line[node[[{#1, #2}]]],
        l, sum1, sum2, n, σ
      },

    l = ArcLength[line];
    sum1 = Total[node[[#1]]];
    sum2 = Total[node[[#2]]];
    n = RotationMatrix[90 Degree].((node[[#2]] - node[[#1]])/l);
    σ = Total[
        f[[{#1, #2}]] {x + y - sum1, -x - y + sum2} / (l h (sum2 - sum1))
      ] 
      IdentityMatrix[2];

      {
        NeumannValue[(n.σ2ε[σ])[[1]], {x, y} ∈ line] (* x component *),
        NeumannValue[(n.σ2ε[σ])[[2]], {x, y} ∈ line] (* y component *)
      }
  ] &;
ΓN = Total@{
    nbc @@ {4, 5} (* Edge 4-5 *),
    nbc @@ {6, 3} (* Edge 6-3 *)
  };

This converts the nodal load to pressure times the shape function of that node, and then finds the normal component.

Dirichlet boundary condition:

ΓD =
  {
    DirichletCondition[{ux[x, y] == 0, uy[x, y] == 0}, {x, y} == node[[1]]],
    DirichletCondition[uy[x, y] == 0, {x, y} == node[[2]]],
    DirichletCondition[uy[x, y] == 0, {x, y} == node[[3]]]
  };

Solving the PDEs

uHat = NDSolveValue[
    {opx == ΓN[[1]], opy + b == ΓN[[2]], ΓD},
    {ux[x, y], uy[x, y]},
    {x, y} ∈ mesh
  ];

which is where the error occurs:

Compile::argcompten: The comparison, LessEqual, is invalid for tensor arguments. Compile::argcompten: The comparison, LessEqual, is invalid for tensor arguments. Compile::argcompten: The comparison, LessEqual, is invalid for tensor arguments. General::stop: Further output of Compile::argcompten will be suppressed during this calculation.

$\endgroup$
  • $\begingroup$ Why are you using a two-element mesh? $\endgroup$ – Alex Trounev Nov 27 '19 at 0:47
  • $\begingroup$ @AlexTrounev It's a course assignment that specifies this mesh. I wanted to try using the FEM package which lead to this problem here. I ended up hard coding the algorithms instead. $\endgroup$ – Yukai Qian Nov 28 '19 at 0:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.