Finite element method for 3D built-up structures using two materials

Question

I wish to do finite element calculations of a steel structure supported on rubber blocks. As a warm up problem I make an all steel structure as follows.

Needs["NDSolve`FEM`"];
Len = 1.0; (*Length of cylinder*)
r = 0.1;  (* Radius of cylinder*)
w2 = 0.12;(* Half width of supports*)
h = 0.1;(* Full height of supports*)
h2 = 0.05;(* Support height under cylinder *)
tk = 0.05; (* Thickness of suports *)

cyl = Cylinder[{{0, 0, 0}, {Len, 0, 0}}, r];
sup1 = Cuboid[{Len/5, -w2, -(r + h2)}, {Len/5 + tk, w2, h - (r + h2)}];
sup2 = Cuboid[{4 Len/5, -w2, -(r + h2)}, {4 Len/5 - tk, w2, 
h - (r + h2)}];
Graphics3D[{cyl, sup1, sup2},
Axes -> True]

The blocks under the cylinder are meant to be rubber but for now I am going to use steel. Make the mesh:

geom = RegionUnion[cyl, sup1, sup2];
mesh = ToElementMesh[geom];
mesh["Wireframe"]

Next define the parameters, variables and that this is a solid mechanics problem. Then I calculate the vibration properties (eigenvalues and vectors). The latter takes about 10 seconds.

pars = <|"YoungModulus" -> 210  10^9, "PoissonRatio" -> 3/10, 
"MassDensity" -> 7850, "AnalysisType" -> "Eigenmode"|>;
vars = {{u[t, x, y, z], v[t, x, y, z], w[t, x, y, z]}, t, {x, y, z}};
ps = SolidMechanicsPDEComponent[vars, pars];
Subscript[\[CapitalGamma], wall] = 
SolidFixedCondition[z == -r - h2, vars, pars];
Timing[{vals, vecs} = 
NDEigensystem[{ps == {0, 0, 0}, Subscript[\[CapitalGamma], 
 wall]}, {u[t, x, y, z], v[t, x, y, z], w[t, x, y, z]}, 
t, {x, y, z} \[Element] mesh, 4];]

Now look at the eigen shapes

These look about right. So far so good. How to include the supports as rubber?

Start of full problem

I start as before but work with OpenCascade

Needs["OpenCascadeLink`"]
Needs["NDSolve`FEM`"];
Len = 1.0; (*Length of cylinder*)
r = 0.1;  (* Radius of cylinder*)
w2 = 0.12;(* Half width of supports*)
h = 0.1;(* Full height of supports*)
h2 = 0.05;(* Support height under cylinder *)
tk = 0.05; (* Thickness of suports *)

cyl = Cylinder[{{0, 0, 0}, {Len, 0, 0}}, r];
sup1 = Cuboid[{Len/5, -w2, -(r + h2)}, {Len/5 + tk, w2, h - (r + h2)}];
sup2 = Cuboid[{4 Len/5, -w2, -(r + h2)}, {4 Len/5 - tk, w2, 
h - (r + h2)}];

I make the supports as separate components and join everything up

supA = RegionDifference[sup1, cyl];
supB = RegionDifference[sup2, cyl];

shape1 = OpenCascadeShape[cyl];
shape2 = OpenCascadeShape[supA];
shape3 = OpenCascadeShape[supB];

faces1 = OpenCascadeShapeFaces[shape1];
faces2 = OpenCascadeShapeFaces[shape2];
faces3 = OpenCascadeShapeFaces[shape3];
union = OpenCascadeShapeUnion[Flatten[{faces1, faces2, faces3}]];

bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[union];
mesh = ToElementMesh[bmesh];
mesh["Wireframe"]

Careful examination of the mesh shows that there is a boundary between the supports and the cylinder. This should enable the material geometry to be modelled properly. The Youngs Modulus of rubber should be about 1.0*10^6 with a Poisson Ratio of 0.4 and a mass density of 1500 but I don't know how to put these in. Looking at other posts it is suggested that a Piecewise construction is preferred but that does not seem possible here. How do I proceed? I could not find a 3D example in Help.

Thanks

Answer 1

welcome back to FEM...

The multi material section seems relevant for this. Here is how I'd do it.

We start in exactly the same way you did:

Needs["OpenCascadeLink`"]
Needs["NDSolve`FEM`"];
Len = 1.0; (*Length of cylinder*)
r = 0.1;  (*Radius of cylinder*)
w2 = 0.12;(*Half width of supports*)
h = 0.1;(*Full height of supports*)
h2 = 0.05;(*Support height under cylinder*)
tk = 0.05; (*Thickness of suports*)

cyl = Cylinder[{{0, 0, 0}, {Len, 0, 0}}, r];
sup1 = Cuboid[{Len/5, -w2, -(r + h2)}, {Len/5 + tk, w2, h - (r + h2)}];
sup2 = Cuboid[{4  Len/5, -w2, -(r + h2)}, {4  Len/5 - tk, w2, 
    h - (r + h2)}];

supA = RegionDifference[sup1, cyl];
supB = RegionDifference[sup2, cyl];

shape1 = OpenCascadeShape[cyl];
shape2 = OpenCascadeShape[supA];
shape3 = OpenCascadeShape[supB];

faces1 = OpenCascadeShapeFaces[shape1];
faces2 = OpenCascadeShapeFaces[shape2];
faces3 = OpenCascadeShapeFaces[shape3];
union = OpenCascadeShapeUnion[Flatten[{faces1, faces2, faces3}]];

bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[union];

Now, I set up a helper association. This will make the code more readable:

material = <|"Rubber" -> 1, "Steel" -> 2|>
(* <|"Rubber" -> 1, "Steel" -> 2|> *)

We insert markers into mesh:

mesh = ToElementMesh[bmesh, 
   "RegionMarker" -> {{RegionCentroid[cyl], 
      material["Steel"]}, {RegionCentroid[sup1], 
      material["Rubber"]}, {RegionCentroid[sup2], 
      material["Rubber"]}}];
mesh["Wireframe"]

And check that the markers are present:

mesh["MeshElementMarkerUnion"]
(* {1, 2} *)

We visualize the individual material groups:

parts = Map[
  mesh["Wireframe"[ElementMarker == #[[1]], 
     "MeshElement" -> "MeshElements", 
     "ElementMeshDirective" -> 
      Directive[EdgeForm[], FaceForm[#[[2]]]]]] &, {{material[
     "Rubber"], Black}, {material["Steel"], Gray}}]

And visualize the combination:

Show[parts]

Now, to the actual question:

vars = {{u[t, x, y, z], v[t, x, y, z], w[t, x, y, z]}, t, {x, y, z}};
pars = <|
    "YoungModulus" -> Piecewise[{
      {210*10^9, ElementMarker == material["Steel"]},
      {1*10^6, ElementMarker == material["Rubber"]}}],
   "PoissonRatio" -> Piecewise[{
      {3/10, ElementMarker == material["Steel"]},
      {4/10, ElementMarker == material["Rubber"]}}],
   "MassDensity" -> Piecewise[{
      {7850, ElementMarker == material["Steel"]},
      {1500, ElementMarker == material["Rubber"]}}], 
   "AnalysisType" -> "Eigenmode"|>;

Set up the equation just like you did:

ps = SolidMechanicsPDEComponent[vars, pars];
Subscript[\[CapitalGamma], wall] = 
  SolidFixedCondition[z == -r - h2, vars, pars];

Get the eigenvalues and functions:

numEigen = 4;
AbsoluteTiming[{vals, vecs} = 
   NDEigensystem[{ps == {0, 0, 0}, 
     Subscript[\[CapitalGamma], wall]}, {u[t, x, y, z], v[t, x, y, z],
      w[t, x, y, z]}, t, {x, y, z} \[Element] mesh, numEigen];]

Compute the frequencies:

eigenfrequencies = Sqrt[vals]/(2 \[Pi])
(* {2.41738, 2.92912, 3.40086, 6.16735} *)

Visualize the result:

GraphicsGrid[Partition[MapThread[Show[
     mesh[
      "Edgeframe"[
       PlotLabel -> Style["freq=" <> ToString[#1] <> " Hz", 8]]],
     deformedMesh2 = ElementMeshDeformation[mesh, #2];
     deformedMesh2["Edgeframe"],
     deformedMesh2["Wireframe"[
       "ElementMeshDirective" -> 
        Directive[EdgeForm[], FaceForm[LightGray]]]]
     ] &, {eigenfrequencies, vecs}], numEigen/2]]

What you be helpful for me to know, is what you tried to figure this out and why you could not find the solution. Where did you look in the help, what would you have expected to find where? That way I can try to improve the documentation.