7
$\begingroup$

Suppose you have an elastic cylinder of radius $R$ and height $H$ and you want to solve the mechanical 3D equilibrium with Mathematica's FEM.

How do you impose a pressure $p$ only on the cylinder sides at $r=R$ (only on the sides, not on the top)?

If possible, I would prefer to stay in the description with Cartesian coordinates.

Cylinder mesh

I am well aware, that a simple analytical solution is easily computable ($u_r(r) = a r + b/r$) for the isotropic case, but I want to extend this problem to some other cases with anisotropic materials and more stuff.

--

Working example with pressure on top

Parameters: geometry, loading and material

(*Geometry [m]*)
R = 1;
H = 5;

(*Pressure [N]*)
p = 10000;

(*Material*)
Em = 2.1*10^9;(*Young's modulus [N/m^2]*)
nu = 0.3;(*Poisson's ratio [-]*)

FEM setup:

Geometry

region = Cylinder[{{0, 0, 0}, {0, 0, H}}, R];
mesh = ToElementMesh[region, "MaxCellMeasure" -> 1.0, 
   "MeshOrder" -> 1];
Show[mesh["Wireframe"], Axes -> True, AxesLabel -> {x1, x2, x3}]

Material (full fourth-order stiffness tensor $\mathbb{C}$ for isotropic material with given Young's modulus Em and Poisson's ratio nu, see parameters)

(*Isotropic material stiffness-fourth-order tensor*)
(*Identities*)
I2 = IdentityMatrix@3;
IdI = TensorProduct[I2, I2];
I4 = TensorTranspose[IdI, {1, 3, 2, 4}];
IS = (I4 + TensorTranspose[I4, {1, 2, 4, 3}])/2;
(*Isotropic projectors*)
P1 = 1/3*IdI;
P2 = IS - P1;
(*Isotropic stiffness*)
C4 = l1*P1 + l2*P2;
l1 = 3*Km;
l2 = 2*Gm;
Km = 1/3*Em/(1 - 2*nu);
Gm = 1/2*Em/(1 + nu);

PDE: linear momentum without body forces ($\mathrm{div}(\boldsymbol{\sigma}) = (0,0,0)_{\{x_1,x_2,x_3\}}$ with $\sigma_{ij} = C_{ijkl} u_{k,l}$)

pde = Table[
   Inactive[
      Div][-C4[[i, ;; , 1, ;;]].Inactive[Grad][
        u1[x1, x2, x3], {x1, x2, x3}], {x1, x2, x3}] + 
    Inactive[
      Div][-C4[[i, ;; , 2, ;;]].Inactive[Grad][
        u2[x1, x2, x3], {x1, x2, x3}], {x1, x2, x3}] + 
    Inactive[
      Div][-C4[[i, ;; , 3, ;;]].Inactive[Grad][
        u3[x1, x2, x3], {x1, x2, x3}], {x1, x2, x3}]
   , {i, 3}
   ];

BCs: Dirichlet to fix bottom ($x_3=0$) and Neumann to impose pressure $p$ on top ($x_3 = H$)

bcD = {
   DirichletCondition[{u3[x1, x2, x3] == 0}, x3 == 0],
   DirichletCondition[{u1[x1, x2, x3] == 0, u2[x1, x2, x3] == 0, 
     u3[x1, x2, x3] == 0}, x1 == 0 && x2 == 0 && x3 == 0]
   };
bcN = {0, 0, NeumannValue[-p, x3 == H]};

Solution for top pressure

{u1sol, u2sol, u3sol} = 
   NDSolveValue[{pde == bcN, bcD}, {u1, u2, u3}, 
    Element[{x1, x2, x3}, mesh]]; // AbsoluteTiming

{0.374435, Null}

Plot deformation

mesh = u1sol["ElementMesh"];
Show[{mesh["Wireframe"], 
  ElementMeshDeformation[mesh, {u1sol, u2sol, u3sol}, 
    "ScalingFactor" -> 10^5][
   "Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]},
  Axes -> True, AxesLabel -> {x1, x2, x3}]

Deformed

$\endgroup$
8
$\begingroup$

You can use something like this:

pred = (0 <= x3 <= H) && (x1^2 + x2^2 >= (R^2 - 10^-2));
bcN = {NeumannValue[-p*Cos[ArcTan[x1, x2]], pred], 
   NeumannValue[-p*Sin[ArcTan[x1, x2]], pred], 0};

Note, that you have to restrict the bcN at the cabs such that there is not force on them. The -10^-2 is for the quadratic case.

A much simpler approach is to use boundary element ElementMarker.

I generated a mesh with:

region = Cylinder[{{0, 0, 0}, {0, 0, H}}, R];
mesh = ToElementMesh[region, "MaxCellMeasure" -> 1.0, 
   "MeshOrder" -> 1, "BoundaryMeshGenerator" -> "OpenCascde"];
groups = mesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
Show[
 mesh["Edgeframe"],
 mesh["Wireframe"["MeshElement" -> "BoundaryElements", 
   "MeshElementStyle" -> (Directive[EdgeForm[], FaceForm[#]] & /@ 
      colors)]]
 , Epilog -> Inset[LineLegend[colors, groups], Scaled[{0.85, 0.8}]]
 ]

enter image description here

This now has an ElementMarker 3 at the boundary in question.

bcN = {NeumannValue[-p*Cos[ArcTan[x1, x2]], ElementMarker == 3], 
   NeumannValue[-p*Sin[ArcTan[x1, x2]], ElementMarker == 3], 0};

Also consider using the "OpenCascade" boundary mesh generator. That will improve the quality of the meshed geometry.

ToElementMesh[region, "BoundaryMeshGenerator" -> "OpenCascde"]

This will be the default from V12.3 onward.

This will give you:

enter image description here

$\endgroup$
10
  • $\begingroup$ Cool, thank you. This is numerically interesting, if you use 0 <= x3 <= H instead, the amplified deformation is quite different at the top and bottom. I suppose that is due to internal numerical reasons, right? $\endgroup$ – Mauricio Fernández Apr 22 at 13:58
  • $\begingroup$ I think it's more of a modeling issue in relation to the DirichletCondition to constraint the object. For this bcD = { DirichletCondition[{u1[x1, x2, x3] == 0, u2[x1, x2, x3] == 0, u3[x1, x2, x3] == 0}, x3 == 0], DirichletCondition[{u1[x1, x2, x3] == 0, u2[x1, x2, x3] == 0}, x3 == H] }; you get a closer to similar deformation for the different pressure NeumannValues. $\endgroup$ – user21 Apr 22 at 14:19
  • $\begingroup$ Hmmm ... your alternative Dirichlet conditions are not what I want in my problem. Your alternative bcD[[1]] is far too restrictive at the bottom (total fix, e.g., you would glue the cylinder to a table) and bcD[[2]] forces u1 = 0 and u2 = 0 at the top. I was thinking of a cylinder under side pressure e.g. minding his own business standing freely on a table (only u3 = 0 at the bottom). From a mechanical point of view, I expected the top to stay flat naturally. From an implementation point of view I have to be careful with 0 <= x3 <= H in the Neumann values from now on. $\endgroup$ – Mauricio Fernández Apr 22 at 16:48
  • $\begingroup$ @MauricioFernández, oh did not want to suggest that you change the BCs I just wanted to point out that there is an interplay between them. To be honest I was also surprised by the effect < vs v<= had. $\endgroup$ – user21 Apr 23 at 6:45
  • $\begingroup$ @MauricioFernández, see update. The problem is that with <= you are also applying the NBC at the cabs. And with < it's not applied everywhere. Sorry it took me so long. I was thinking in the wrong direction. $\endgroup$ – user21 Apr 23 at 10:53

Your Answer

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

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