10
$\begingroup$

For stress analysis in the finite element method you need a stress operator. For two dimensional plane stress it can be found here. For two dimensional plain strain it can be found here. I need the two dimensional case where there is axial symmetry so that in cylindrical coordinates this is displacements and stress in terms of {r,z}. Getting these stress operators is not so straight forward since one needs a careful application of Inactive. I also need the calculation of stress from the displacement fields. Can you help please?

For stress in 3D the stress operator may be found here.

Here are the stress operators and the method for calculating stress for plane stress and plane strain. I also give vonMises stress for plane stress

ClearAll[planeStressOperator, Y, ν];
planeStressOperator[{u_, v_}, {x_, 
   y_}, {Y_, ν_}] := {Inactive[
     Div][({{0, -((Y ν)/(1 - ν^2))}, {-((Y (1 - ν))/(2 (1 \
- ν^2))), 0}}.Inactive[Grad][v[x, y], {x, y}]), {x, y}] + 
   Inactive[
     Div][({{-(Y/(1 - ν^2)), 
        0}, {0, -((Y (1 - ν))/(2 (1 - ν^2)))}}.Inactive[Grad][
       u[x, y], {x, y}]), {x, y}], 
  Inactive[Div][({{0, -((Y (1 - ν))/(2 (1 - ν^2)))}, {-((Y \
ν)/(1 - ν^2)), 0}}.Inactive[Grad][u[x, y], {x, y}]), {x, y}] +
    Inactive[
     Div][({{-((Y (1 - ν))/(2 (1 - ν^2))), 
        0}, {0, -(Y/(1 - ν^2))}}.Inactive[Grad][
       v[x, y], {x, y}]), {x, y}]};


 ClearAll[planeStrainOperator];
planeStrainOperator[{u_, v_}, {x_, 
    y_}, {Y_, ν_}] := {Inactive[
      Div][({{0, -((Y ν)/((1 - 
                2 ν) (1 + ν)))}, {-(Y/(2 (1 + ν))), 
         0}}.Inactive[Grad][v[x, y], {x, y}]), {x, y}] + 
    Inactive[
      Div][({{-((Y (1 - ν))/((1 - 2 ν) (1 + ν))), 
         0}, {0, -(Y/(2 (1 + ν)))}}.Inactive[Grad][
        u[x, y], {x, y}]), {x, y}], 
   Inactive[
      Div][({{0, -(Y/(2 (1 + ν)))}, {-((Y ν)/((1 - 
                2 ν) (1 + ν))), 0}}.Inactive[Grad][
        u[x, y], {x, y}]), {x, y}] + 
    Inactive[
      Div][({{-(Y/(2 (1 + ν))), 
         0}, {0, -((Y (1 - ν))/((1 - 
                2 ν) (1 + ν)))}}.Inactive[Grad][
        v[x, y], {x, y}]), {x, y}]};


ClearAll[stress2D]
stress2D[{uif_InterpolatingFunction, 
   vif_InterpolatingFunction}, {Y_, ν_}, ps_] :=
 Block[{fac, dd, df, mesh, coords, dv, ux, uy, vx, vy, ex, ey, gxy, 
   sxx, syy, sxy},
  dd = Outer[(D[#1[x, y], #2]) &, {uif, vif}, {x, y}];
  Which[
   ps == "PlaneStress", 
   fac = Y/(1 - ν^2)*{{1, ν, 0}, {ν, 1, 0}, {0, 
       0, (1 - ν)/2}},
   ps == "PlaneStrain", 
   fac = Y/((1 + ν) (1 - 2 ν)) *{{1 - ν, ν, 
       0}, {ν, 1 - ν, 0}, {0, 0, (1 - 2 ν)/2}}
   ];

  df = Table[Function[{x, y}, Evaluate[dd[[i, j]]]], {i, 2}, {j, 2}];
  (*the coordinates from the ElementMesh*)
  mesh = uif["Coordinates"][[1]];
  coords = mesh["Coordinates"];
  dv = Table[df[[i, j]] @@@ coords, {i, 2}, {j, 2}];
  ux = dv[[1, 1]];
  uy = dv[[1, 2]];
  vx = dv[[2, 1]];
  vy = dv[[2, 2]];
  ex = ux;
  ey = vy;
  gxy = (uy + vx);
  sxx = fac[[1, 1]]*ex + fac[[1, 2]]*ey;
  syy = fac[[2, 1]]*ex + fac[[2, 2]]*ey;
  sxy = fac[[3, 3]]*gxy;
  {ElementMeshInterpolation[{mesh}, sxx],
   ElementMeshInterpolation[{mesh}, syy],
   ElementMeshInterpolation[{mesh}, sxy]
   }]

ClearAll[vonMisesPlaneStress];
vonMisesPlaneStress[{sxx_, syy_, sxy_}] := 
 Function[{x, y}, Sqrt[
  sxx[x, y]^2 - sxx[x, y] syy[x, y] + syy[x, y]^2 + 3 sxy[x, y]^2]]

Here is an example of their use. I am looking at stress concentrations near a cut-out. I include a MeshRefinementFunction at the cut-out which I am using to show that improving the mesh does not help with the stress concentration in this case.

Needs["NDSolve`FEM`"];
L = 2; (* Length *)
h = 1; (* Height *)
Y = 10^3; (* Modulus of elasticity *)
ν = 33/100;  (* Poisson ratio *)
stress = 10;  (* stress *)
gls = 0.1; (* Grid length scale *)
nr = 10;  (* Factor for refinement *)
r0 = 2 gls;  (* radius for refinement *)


ClearAll[cf];
cf = Compile[{{c, _Real, 2}, {a, _Real, 0}},
   Block[{com, c1, c2, r1, r2, r},
    com = Total[c]/3;
    c1 = {0.2 L, 0.6 h};
    c2 = {0.4 L, 0.6 h};
    r1 = Norm[com - c1];
    r2 = Norm[com - c2];
    r = Min[{r1, r2}];

    If[r < r0 && a > (Sqrt[3]/4) (gls/nr + r/r0 gls)^2, True, False]

    ]
   ];

mesh = ToElementMesh[ 
   RegionDifference[Rectangle[{0, 0}, {L, h}], 
    RegionUnion[Rectangle[{0.2 L, 0.6 h}, {0.4 L, h}], 
     Disk[{0.75 L, 0.4 h}, 0.2 h]]],
   MaxCellMeasure -> {"Length" -> .1}, MeshRefinementFunction -> cf];

Show[mesh["Wireframe"],
  Graphics[{
   {Thickness[0.01], Line[{{0, -0.2 h}, {0, 1.2 h}}]},
   Table[Arrow[{{L, y}, {L + L/2 (h - y), y}}], {y, 0, h, h/12}],
   Table[Arrow[{{x, 3 h/2}, {x, h}}], {x, 3 L/4, L, L/20}]
   }]
 ]

Mathematica graphics

I can solve for plane stress and work out the displacements as follows.

{uif, vif} = NDSolveValue[{
    planeStressOperator[{u, v}, {x, y}, {Y, ν}] == {
      NeumannValue[stress (h - y), x == L && 0 <= y <= h],
      NeumannValue[-stress, 3 L/4 <= x <= L && y == h]},
    DirichletCondition[u[x, y] == 0, x == 0], 
    DirichletCondition[v[x, y] == 0, x == 0]
    }, {u, v}, {x, y} ∈ mesh];



  dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
ContourPlot[uif[x, y], {x, y} ∈ mesh, 
 AspectRatio -> Automatic, ColorFunction -> "Temperature", 
 PlotRange -> All]
Show[{
  mesh["Wireframe"],
  dmesh["Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]

Mathematica graphics

Mathematica graphics

Finally I calculate the stresses and the von Mises stress.

{σxx, σyy, σxy} = 
  stress2D[{uif, vif}, {Y, ν}, "PlaneStress"];
vm = vonMisesPlaneStress[{σxx, σyy, σxy}];

Plot3D[vm[x, y], {x, y} ∈ mesh, BoxRatios -> {L, h, h}, 
 PlotRange -> All]

Mathematica graphics

All this shows that the finite element method is looking to be in good shape. What I need to complete my stress analysis is the axisymmetric stress operator for {r, z} and a method of calculating the stresses. Has anyone done this?

I have found a good reference for this problem. To help get started here are the strain-displacement and strain-stress equations for three dimensions { r, θ, z}.

$\epsilon _{\text{rr}}=\frac{\partial u_r}{\partial r}$

$\epsilon _{\theta \theta }=\frac{\partial u_{\theta }}{r\partial \theta }+\frac{u_r}{r}$

$\epsilon _{\text{r$\theta $}}=\frac{\partial u_{\theta }}{\partial r}-\frac{u_{\theta }}{r}+\frac{\partial u_r}{r\partial \theta }$

$\epsilon _{\text{zz}}=\frac{\partial u_z}{\partial z}$

$\epsilon _{\text{$\theta $z}}=\frac{\partial u_{\theta }}{\partial z}+\frac{\partial u_z}{r\partial \theta }$

$\epsilon _{\text{rz}}=\frac{\partial u_r}{\partial z}+\frac{\partial u_z}{\partial r}$

$Y \epsilon _{\text{rr}}=\sigma _{\text{rr}}-\nu \left(\sigma _{\theta \theta }+\sigma _{\text{zz}}\right)$

$Y \epsilon _{\theta \theta }=\sigma _{\theta \theta }-\nu \left(\sigma _{\theta \theta }+\sigma _{\text{zz}}\right)$

$Y \epsilon _{\text{zz}}=\sigma _{\text{zz}}-\nu \left(\sigma _{\theta \theta }+\sigma _{\text{rr}}\right))$

$\tau _{\text{r$\theta $}}=G \epsilon _{\text{r$\theta $}}$

$\tau _{\text{rz}}=G \epsilon _{\text{rz}}$

$\tau _{\text{$\theta $z}}=G \epsilon _{\text{$\theta $z}}$

Here Y, G and ν are the modulus of elasticiy, the shear modulus and Poisson ratio respectively.

Putting these into Mathematica format gives

strainDisplacement = {
   err[r, θ, z] == D[ur[r, θ, z], r],
   eθθ[r, θ, z] == 
    D[uθ[r, θ, z], θ]/r + ur[r, θ, z]/r,
   erθ[r, θ, z] == 
    D[uθ[r, θ, z], r]/r - uθ[r, θ, z]/r + 
     D[ur[r, θ, z], θ]/r,
   ezz[r, θ, z] == D[uz[r, θ, z], z],
   eθz[r, θ, z] == 
    D[uθ[r, θ, z], z] + D[uz[r, θ, z], θ]/
     r,
   erz[r, θ, z] == 
    D[ur[r, θ, z], z] + D[uz[r, θ, z], z] + 
     D[uz[r, θ, z], r]
   };
strainStress = {
   Y err[r, θ, z] == σrr[r, θ, 
      z] - ν (σθθ[r, θ, 
         z] + σzz[r, θ, z]),
   Y eθθ[r, θ, z] == σθθ[
      r, θ, 
      z] - ν (σzz[r, θ, 
         z] + σθθ[r, θ, z]),
   Y ezz[r, θ, z] == σzz[r, θ, 
      z] - ν (σrr[r, θ, 
         z] + σθθ[r, θ, z]),
   τrθ[r, θ, z] == G erθ[r, θ, z], 
   τrz[r, θ, z] == G erz[r, θ, z],
   τθz[r, θ, z] == G eθz[r, θ, z]
   };

I guess that as I want a two dimensional axisymmetric {r,z} formulation all the derivatives with respect to θ will be zero. Not sure where to go then... Somehow we need to get the equations in the same form as above involving just derivatives of the two displacements ur and uz. Also we need the Inactive in the correct place.

$\endgroup$
5
  • 1
    $\begingroup$ Hugh, I am confused whether you are asking about the fundamental physics, or how to implement them in Mathematica. Do you have an analytical expression / a starting point for the stress operator you are interested in? $\endgroup$
    – MarcoB
    Apr 30, 2018 at 19:11
  • 1
    $\begingroup$ I need the mathematica implementation. I will look to see if I can find the analytical expression. $\endgroup$
    – Hugh
    Apr 30, 2018 at 19:12
  • $\begingroup$ @MarcoB I have added the equations which I hope will give you a start. $\endgroup$
    – Hugh
    Apr 30, 2018 at 20:09
  • $\begingroup$ I have a question how do you deal with the boundary condition of shear stress? $\endgroup$ Mar 10, 2020 at 1:52
  • $\begingroup$ @PleaseCorrectGrammarMistakes You can put in tangential as well as normal stresses. However, watch the directions very carefully. $\endgroup$
    – Hugh
    Mar 10, 2020 at 12:51

1 Answer 1

3
$\begingroup$

This question has been around unanswered for quite a while. In the mean time much has happened and the SolidMechanicsPDEComponent can be used to set up axisymmetric solid mechanics models. I'll show the model first, and then I'll add a few comments at the end.

Create the geometry:

Needs["NDSolve`FEM`"];
geometricParameters = <|L -> 2, h -> 1|>;
region = 
  RegionDifference[Rectangle[{0, 0}, {L, h}], 
    RegionUnion[Rectangle[{0.2  L, 0.6  h}, {0.4  L, h}], 
     Disk[{0.75  L, 0.4  h}, 0.2  h]]] /. geometricParameters;

Create the mesh:

(meshAxi = ToElementMesh[region])["Wireframe"]

enter image description here

Add the variables and parameters:

varsAxi = {{u[r, z], 0, w[r, z]}, {r, \[Theta], z}};
parsAxi = 
  Join[<|"RegionSymmetry" -> "Axisymmetric", "YoungModulus" -> 10^4, 
    "PoissonRatio" -> 33/100|>, geometricParameters];

Note the specific way I set up the axisymmetric case in that I used 0 as the displacement in the azimuthal direction and theta as a independent variable. There are other ways to set this up. But this allows us to get the 0 displacement in the v-direction.

Solve the PDE:

displacementAxi = 
 NDSolveValue[{SolidMechanicsPDEComponent[varsAxi, parsAxi] == 
    SolidBoundaryLoadValue[r == L, varsAxi, 
      parsAxi, <|"Pressure" -> {10  (h - z), 0, 0}|>] + 
     SolidBoundaryLoadValue[z == h && 3 L/4 <= r <= L, varsAxi, 
      parsAxi, <|"Pressure" -> {0, 0, -10}|>]
   , SolidDisplacementCondition[z == 0 && r < L/10, varsAxi, 
    parsAxi, <|"Displacement" -> {0, None, 0}|>]
   , SolidDisplacementCondition[z == h && r < L/10, varsAxi, 
    parsAxi, <|"Displacement" -> {0, None, 0}|>]
   }, varsAxi[[1]], {r, z} \[Element] meshAxi
  ]

Visualize the solution:

Show[
 meshAxi["Wireframe"["MeshElement" -> "BoundaryElements"]],
 ElementMeshDeformation[meshAxi, displacementAxi[[{1, 3}]], 
   "ScalingFactor" -> 5]["Wireframe"]]

enter image description here

Note that I changed the boundary conditions a bit (and the Young modulus value). For an axisymmetric case you need to be careful with the boundary conditions at r==0 since that is actually within the rotational region and you can not really clamp the object there. So I choose to hold it fixed in a small part at the top and bottom. Hope that makes sense.

Compute and visualize the vonMises stress:

strain = SolidMechanicsStrain[varsAxi, parsAxi, displacementAxi];
stress = 
  SolidMechanicsStress[varsAxi, parsAxi, strain, displacementAxi];
vms = Head[PDEModels`VonMisesStress[varsAxi, parsAxi, stress]];
ContourPlot[vms[r, z], {r, z} \[Element] meshAxi, PlotRange -> All, 

AspectRatio -> Automatic, Contours -> 25]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.