Solid Mechanics FEM Simulation with Different Material Properties

Question

How would I assign different material properties to the "bar" and "support"? Meaning, for example, the bar would be assigned the properties of single-crystal Copper and the support single-crystal Nickel.

Update: More specifically, "MaterialModel" -> "LinearElasticAnisotropic" is needed. The If statement works for properties like "MassDensity", but I've been unsuccessful with the compliance matrix.

Update 2: Thanks to User21 comment about IF statement usage, the compliance matrix now has IF statements at every entry, and it runs. The issue now is that the reported frequencies are as if the solid part is not constrained in space. It's constrained the same way as a monolithic material which gives appropriate values. Any thoughts?

Update 3: It was a units (1*^-12) error that caused very low frequencies that made me think it was a constraint issue. Fixed in the code below (answer) and it seems correct.

Update 4: Runs significantly slower compared to the monolithic material, even with a coarse mesh. Any suggestions to speed it up would be appreciated. Please comment if a bounty would motivate anyone to take the code from my "answer" and speed it up for the first 60 Eigenmodes.

Answer 1

Update:

I recommend using Piecewise and added an example and an explanation below.

This is probably something that should be shown in the documentation and I'll add this. Here is how to do it. The compliance matrix needs to be given as a matrix. So using paradigms like If[predicate, matA, matB] do not work. The trick is to use the If (or even better Piecewise) statement inside the matrix at each entry. Here, however, we have to be careful that the entries are evaluated to numbers such that the evaluation can work efficiently. What we do not want is an entry like If[predicate, matA[[1,1]], matB[[1,1]]] In this case the evaluator will will have to look up the matrix values for each integration point. This is very inefficient. What we want is something like If[predicate, valueA11, valueB11] where the values are actual numbers. We use the compliance matrices for two materials:

matSNickel = {{5.00, -1.55, -1.55, 0.00, 0.00, 0.00}, {-1.55, 
    5.00, -1.55, 0.00, 0.00, 0.00}, {-1.55, -1.55, 5.00, 0.00, 0.00, 
    0.00}, {0.00, 0.00, 0.00, 13.10, 0.00, 0.00}, {0.00, 0.00, 0.00, 
    0.00, 13.10, 0.00}, {0.00, 0.00, 0.00, 0.00, 0.00, 13.10}};
matSCopper = {{7.69, -2.62, -2.62, 0.00, 0.00, 0.00}, {-2.62, 
    7.69, -2.62, 0.00, 0.00, 0.00}, {-2.62, -2.62, 7.69, 0.00, 0.00, 
    0.00}, {0.00, 0.00, 0.00, 20.62, 0.00, 0.00}, {0.00, 0.00, 0.00, 
    0.00, 20.62, 0.00}, {0.00, 0.00, 0.00, 0.00, 0.00, 20.62}};

Here is how we can then set up a single compliance matrix based on an If statement that uses a predicate where each of the components are active:

cmat = MapThread[
   If[Abs[y] <= Abs[barY/2], #1, #2] &, {matSNickel, 
     matSCopper}*1*^-12, 2];

pars = <|"MaterialModel" -> "LinearElasticAnisotropic", 
   "MassDensity" -> If[Abs[y] <= Abs[barY/2], 8960, 8908], 
   "ComplianceMatrix" -> cmat|>;

Note that for example

cmat[[1, 1]]

(* If[Abs[y] <= 0.0015, 5.*10^-12, 7.69*10^-12] *)

This is very efficient. An even better approach is to use Piecewise, because Piecewise does not have the HoldRest attribute the If statement has. For this reason the expressions in the branch will automatically be evaluated. So the recommended expression to use is \

cmat = MapThread[
   Piecewise[{{#1, Abs[y] <= Abs[barY/2]}}, #2] &, {matSNickel, 
     matSCopper}*1*^-12, 2];

Hope this helps.

Answer 2

This is answer for the first question about titanium and nickel usage in one installation, and for the last question about how to speed up evaluation. First step is a mesh generation and visualization with FEM

Needs["NDSolve`FEM`"]

barX = 0.0100;
barY = 0.0015;
barZ = 0.1000;

supportX = 0.0750;
supportY = 0.0050;
supportZ = 0.0010;

barSolid = 
  Cuboid[{-(barX/2), -(barY/2), -(barZ/2)}, {(barX/2), (barY/
      2), (barZ/2)}];
supportASolid = 
  Cuboid[{-(supportX/2), -supportY - (barY/2), -(supportZ/
       2)}, {(supportX/2), -(barY/2), (supportZ/2)}];
supportBSolid = 
  Cuboid[{-(supportX/2), (barY/2), -(supportZ/2)}, {(supportX/
      2), (barY/2) + supportY, (supportZ/2)}];
reg = RegionUnion[barSolid, supportASolid, supportBSolid];


mesh = ToElementMesh[reg, "MeshOrder" -> 1, MaxCellMeasure -> 10^-6];
mesh["Wireframe"]

Second step is material property definition

en1 = Entity["Element", "Titanium"];

en2 = Entity["Element", "Nickel"];

param1 = 
 QuantityMagnitude[
  en1[{"YoungModulus", "PoissonRatio", "MassDensity"}]]

param2 = 
 QuantityMagnitude[
  en2[{"YoungModulus", "PoissonRatio", "MassDensity"}]]

Please pay attention that en1, en2 have mixed bases so that if we evaluate for example en1 then we have

en1[{"YoungModulus", "PoissonRatio", "MassDensity"}]

Out[]= {Quantity[116., "Gigapascals"], 0.32, 
 Quantity[4.507, ("Grams")/("Centimeters")^3]}

For frequencies simulation we need Young modulus in Pa, while mass density in $kg/m^3$.Therefore we define parameters in different regions as follows

 cons = (-barX/2 <= x <= barX/2 && -barY/2 <= y <= 
    barY/2 && -barZ/2 <= z <= barZ/2); Y = 
 With[{Y1 = 10^9 param1[[1]], Y2 = 10^9 param2[[1]]}, 
  If[cons, Y1, Y2]]; \[Nu] = 
 With[{n1 = param1[[2]], n2 = param2[[2]]}, 
  If[cons, n1, n2]]; \[Rho] = 
 With[{r1 = 10^3 param1[[3]], r2 = 10^3 param2[[3]]}, 
  If[cons, r1, r2]];

Here we use With since we can't use If directly with FEM due to compilation problem. Last step

tvars = {{u[t, x, y, z], v[t, x, y, z], w[t, x, y, z]}, 
   t, {x, y, z}};
constraintA = 
  DirichletCondition[{u[t, x, y, z] == 0, v[t, x, y, z] == 0, 
    w[t, x, y, z] == 
     0}, (x^2 <= (barX/2)^2 && y^2 == (barY/2)^2 && 
      z^2 <= (supportZ/2)^2) || (x^2 >= (supportX/2)^2 && 
      y^2 >= (barY/2)^2 && z^2 <= (supportZ/2)^2)];
op = SolidMechanicsPDEComponent[{{u[t, x, y, z], v[t, x, y, z], 
     w[t, x, y, z]}, t, {x, y, z}}, 
   Join[<|"YoungModulus" -> Y, "PoissonRatio" -> \[Nu], 
     "MassDensity" -> \[Rho]|>]];

{evals, evecs} = 
  NDEigensystem[{op == {0, 0, 0}, constraintA}, tvars[[1]], 
   t, {x, y, z} \[Element] mesh, 6];
 Re[Sqrt[evals]/(2 Pi)]

(*Out[]= {9.90091, 9.90091, 10.1875, 10.1875, 16.6058, 16.6058}*)

Visualization of first and fifth modes

{DensityPlot3D[evecs[[1, 1]], {x, y, z} \[Element] mesh, 
  ColorFunction -> "Rainbow", OpacityFunction -> None, Boxed -> False,
   Axes -> False, BoxRatios -> Automatic, 
  PlotLabel -> "Titanium-Nickel"], 
 DensityPlot3D[evecs[[5, 1]], {x, y, z} \[Element] mesh, 
  ColorFunction -> "Rainbow", OpacityFunction -> None, Boxed -> False,
   Axes -> False, BoxRatios -> Automatic, 
  PlotLabel -> "Titanium-Nickel"]}

Update 1. The answer for the second question about copper-nickel installation:

Needs["NDSolve`FEM`"]

barX = 0.0100;
barY = 0.0030;
barZ = 0.1000;

supportX = 0.06000;
supportY = 0.00150;
supportZ = 0.00100;

barSolid = 
  Cuboid[{-(barX/2), -(barY/2), -(barZ/2)}, {(barX/2), (barY/
      2), (barZ/2)}];
supportASolid = 
  Cuboid[{-(supportX/2), -supportY - (barY/2), -(supportZ/
       2)}, {(supportX/2), -(barY/2), (supportZ/2)}];
supportBSolid = 
  Cuboid[{-(supportX/2), (barY/2), -(supportZ/2)}, {(supportX/
      2), (barY/2) + supportY, (supportZ/2)}];

solid = RegionUnion[barSolid, supportASolid, supportBSolid];

mesh = ToElementMesh[solid, "MeshOrder" -> 1, 
  MaxCellMeasure -> {"Volume" -> 5*^-10}, MeshQualityGoal -> 1, 
  "MeshElementType" -> "HexahedronElement"]

en1 = Entity["Element", "Copper"];

en2 = Entity["Element", "Nickel"]; 
 
param1 = 
 QuantityMagnitude[
  en1[{"YoungModulus", "PoissonRatio", "MassDensity", "ShearModulus"}]]

param2 = 
 QuantityMagnitude[
  en2[{"YoungModulus", "PoissonRatio", "MassDensity", "ShearModulus"}]]

cons = (-barX/2 <= x <= barX/2 && -barY/2 <= y <= barY/2 && -barZ/2 <=
     z <= barZ/2); Y1 = 
 With[{y1 = 10^9 param1[[1]], y2 = 10^9 param2[[1]]}, 
  If[cons, y1, y2]]; G1 = 
 With[{g1 = 10^9 param1[[4]], g2 = 10^9 param2[[4]]}, 
  If[cons, g1, g2]]; \[Nu]1 = 
 With[{n1 = param1[[2]], n2 = param2[[2]]}, 
  If[cons, n1, n2]]; \[Rho]1 = 
 With[{r1 = 10^3 param1[[3]], r2 = 10^3 param2[[3]]}, 
  If[cons, r1, r2]];
param1 = 
 QuantityMagnitude[
  en1[{"YoungModulus", "PoissonRatio", "MassDensity", "ShearModulus"}]]

param2 = 
 QuantityMagnitude[
  en2[{"YoungModulus", "PoissonRatio", "MassDensity", "ShearModulus"}]]

cons = (-barX/2 <= x <= barX/2 && -barY/2 <= y <= barY/2 && -barZ/2 <=
     z <= barZ/2); Y1 = 
 With[{y1 = 10^9 param1[[1]], y2 = 10^9 param2[[1]]}, 
  If[cons, y1, y2]]; G1 = 
 With[{g1 = 10^9 param1[[4]], g2 = 10^9 param2[[4]]}, 
  If[cons, g1, g2]]; \[Nu]1 = 
 With[{n1 = param1[[2]], n2 = param2[[2]]}, 
  If[cons, n1, n2]]; \[Rho]1 = 
 With[{r1 = 10^3 param1[[3]], r2 = 10^3 param2[[3]]}, 
  If[cons, r1, r2]];
(*compliance matrix in general case*)
D0 = {{1/Subscript[Y, x*x], -(Subscript[\[Nu], y*x]/Subscript[Y, y*y]), -(Subscript[\[Nu], z*x]/Subscript[Y, z*z]), 0, 0, 0}, 
    {-(Subscript[\[Nu], x*y]/Subscript[Y, x*x]), 1/Subscript[Y, y*y], -(Subscript[\[Nu], z*y]/Subscript[Y, z*z]), 0, 0, 0}, 
    {-(Subscript[\[Nu], x*z]/Subscript[Y, x*x]), -(Subscript[\[Nu], y*z]/Subscript[Y, y*y]), 1/Subscript[Y, z*z], 0, 0, 0}, 
    {0, 0, 0, 1/Subscript[G, y*z], 0, 0}, {0, 0, 0, 0, 1/Subscript[G, x*z], 0}, {0, 0, 0, 0, 0, 1/Subscript[G, x*y]}}; Compliance matrix in this case
mat = D0 /. {Subscript[Y, x*x] -> Y1, Subscript[Y, y*y] -> Y1, 
    Subscript[Y, z*z] -> Y1, Subscript[\[Nu], y*x] -> \[Nu]1, 
    Subscript[\[Nu], z*x] -> \[Nu]1, Subscript[\[Nu], x*y] -> \[Nu]1, 
    Subscript[\[Nu], z*y] -> \[Nu]1, Subscript[\[Nu], x*z] -> \[Nu]1, 
    Subscript[\[Nu], y*z] -> \[Nu]1, Subscript[G, y*z] -> G1, 
    Subscript[G, x*z] -> G1, Subscript[G, x*y] -> G1};

Compliance matrix test

mat /. {x -> 0, y -> 0, z -> 0}

(*Out[]= {{7.69*10^-12, -2.6*10^-12, -2.6*10^-12, 0, 0, 
  0}, {-2.6*10^-12, 7.69*10^-12, -2.6*10^-12, 0, 0, 
  0}, {-2.6*10^-12, -2.6*10^-12, 7.69*10^-12, 0, 0, 0}, {0, 0, 0, 
  2.1*10^-11, 0, 0}, {0, 0, 0, 0, 2.1*10^-11, 0}, {0, 0, 0, 0, 0, 
  2.1*10^-11}}*)

mat /. {x -> 0, y -> (barY/2) + supportY/2, z -> 0}

(*Out[]= {{5.0*10^-12, -1.6*10^-12, -1.6*10^-12, 0, 0, 
  0}, {-1.6*10^-12, 5.0*10^-12, -1.6*10^-12, 0, 0, 
  0}, {-1.6*10^-12, -1.6*10^-12, 5.0*10^-12, 0, 0, 0}, {0, 0, 0, 
  1.3*10^-11, 0, 0}, {0, 0, 0, 0, 1.3*10^-11, 0}, {0, 0, 0, 0, 0, 
  1.3*10^-11}}*)

Eigen system computation

tvars = {{u[t, x, y, z], v[t, x, y, z], w[t, x, y, z]}, t, {x, y, z}};

pars = <|"MaterialModel" -> "LinearElasticAnisotropic", 
   "MassDensity" -> \[Rho]1, "ComplianceMatrix" -> mat|>;

eigenmodeOperator = 
  SolidMechanicsPDEComponent[tvars, 
   Join[pars, <|"AnalysisType" -> "Eigenmode"|>]];

constraintB = 
  SolidFixedCondition[
   x^2 >= (supportX/2)^2 && y^2 >= (barY/2)^2 && 
    z^2 <= (supportZ/2)^2, tvars, pars];

{evals, evecs} = 
  NDEigensystem[{eigenmodeOperator == {0, 0, 0}, constraintB}, 
   tvars[[1]], t, {x, y, z} \[Element] mesh, 6];

eigenfrequencies = Re[Sqrt[evals]/(2 Pi)]

(*Out[]= {155.855, 319.74, 513.075, 551.638, 1642.23, 2346.6}*)

Visualization

 {DensityPlot3D[
  evecs[[1, 1]] // Re // Evaluate, {x, y, z} \[Element] mesh, 
  ColorFunction -> "Rainbow", OpacityFunction -> None, Boxed -> False,
   Axes -> False, BoxRatios -> Automatic, 
  PlotLabel -> "Copper-Nickel"], 
 DensityPlot3D[
  evecs[[5, 1]] // Re // Evaluate, {x, y, z} \[Element] mesh, 
  ColorFunction -> "Rainbow", OpacityFunction -> None, Boxed -> False,
   Axes -> False, BoxRatios -> Automatic, 
  PlotLabel -> "Copper-Nickel"]}

Answer 3

With the help of user21, the following code simulates the Eigenmodes of a solid component with two different densities and compliance matrices based on position within the mesh.

Clear["Global`*"]

barX = 0.0100;
barY = 0.0030;
barZ = 0.1000;

supportX = 0.06000;
supportY = 0.00150;
supportZ = 0.00100;

barSolid = 
  Cuboid[{-(barX/2), -(barY/2), -(barZ/2)}, {(barX/2), (barY/
      2), (barZ/2)}];
supportASolid = 
  Cuboid[{-(supportX/2), -(barY/2), -(supportZ/2)}, {(supportX/
      2), -supportY - (barY/2), (supportZ/2)}];
supportBSolid = 
  Cuboid[{-(supportX/2), (barY/2), -(supportZ/2)}, {(supportX/
      2), (barY/2) + supportY, (supportZ/2)}];

solid = RegionUnion[barSolid, supportASolid, supportBSolid];

mesh = ToElementMesh[solid, "MeshOrder" -> 1, 
  MaxCellMeasure -> {"Volume" -> 5*^-10}, MeshQualityGoal -> 1, 
  "MeshElementType" -> "HexahedronElement"]

matSNickel = {
   {5.00, -1.55, -1.55, 0.00, 0.00, 0.00},
   {-1.55, 5.00, -1.55, 0.00, 0.00, 0.00},
   {-1.55, -1.55, 5.00, 0.00, 0.00, 0.00},
   {0.00, 0.00, 0.00, 13.10, 0.00, 0.00},
   {0.00, 0.00, 0.00, 0.00, 13.10, 0.00},
   {0.00, 0.00, 0.00, 0.00, 0.00, 13.10}};
matSCopper = {
   {7.69, -2.62, -2.62, 0.00, 0.00, 0.00},
   {-2.62, 7.69, -2.62, 0.00, 0.00, 0.00},
   {-2.62, -2.62, 7.69, 0.00, 0.00, 0.00},
   {0.00, 0.00, 0.00, 20.62, 0.00, 0.00},
   {0.00, 0.00, 0.00, 0.00, 20.62, 0.00},
   {0.00, 0.00, 0.00, 0.00, 0.00, 20.62}};

matSMulti[loc_] := {
  {If[Abs[loc] <= Abs[barY/2], matSCopper[[1, 1]], 
    matSNickel[[1, 1]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[1, 2]], 
    matSNickel[[1, 2]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[1, 3]], 
    matSNickel[[1, 3]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[1, 4]], 
    matSNickel[[1, 4]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[1, 5]], 
    matSNickel[[1, 5]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[1, 6]], 
    matSNickel[[1, 6]]]},
  {If[Abs[loc] <= Abs[barY/2], matSCopper[[2, 1]], 
    matSNickel[[2, 1]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[2, 2]], 
    matSNickel[[2, 2]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[2, 3]], 
    matSNickel[[2, 3]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[2, 4]], 
    matSNickel[[2, 4]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[2, 5]], 
    matSNickel[[2, 5]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[2, 6]], 
    matSNickel[[2, 6]]]},
  {If[Abs[loc] <= Abs[barY/2], matSCopper[[3, 1]], 
    matSNickel[[3, 1]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[3, 2]], 
    matSNickel[[3, 2]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[3, 3]], 
    matSNickel[[3, 3]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[3, 4]], 
    matSNickel[[3, 4]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[3, 5]], 
    matSNickel[[3, 5]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[3, 6]], 
    matSNickel[[3, 6]]]},
  {If[Abs[loc] <= Abs[barY/2], matSCopper[[4, 1]], 
    matSNickel[[4, 1]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[4, 2]], 
    matSNickel[[4, 2]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[4, 3]], 
    matSNickel[[4, 3]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[4, 4]], 
    matSNickel[[4, 4]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[4, 5]], 
    matSNickel[[4, 5]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[4, 6]], 
    matSNickel[[4, 6]]]},
  {If[Abs[loc] <= Abs[barY/2], matSCopper[[5, 1]], 
    matSNickel[[5, 1]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[5, 2]], 
    matSNickel[[5, 2]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[5, 3]], 
    matSNickel[[5, 3]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[5, 4]], 
    matSNickel[[5, 4]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[5, 5]], 
    matSNickel[[5, 5]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[5, 6]], 
    matSNickel[[5, 6]]]},
  {If[Abs[loc] <= Abs[barY/2], matSCopper[[6, 1]], 
    matSNickel[[6, 1]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[6, 2]], 
    matSNickel[[6, 2]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[6, 3]], 
    matSNickel[[6, 3]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[6, 4]], 
    matSNickel[[6, 4]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[6, 5]], 
    matSNickel[[6, 5]]], 
   If[Abs[loc] <= Abs[barY/2], matSCopper[[6, 6]], matSNickel[[6, 6]]]}
  }

tvars = {{u[t, x, y, z], v[t, x, y, z], w[t, x, y, z]}, t, {x, y, z}};

pars = <|
   "MaterialModel" -> "LinearElasticAnisotropic",
   "MassDensity" -> If[Abs[y] <= Abs[barY/2], 8960, 8908],
   "ComplianceMatrix" -> matSMulti[y]*1*^-12
   |>;

eigenmodeOperator = 
  SolidMechanicsPDEComponent[tvars, 
   Join[pars, <|"AnalysisType" -> "Eigenmode"|>]];

constraintB = 
  SolidFixedCondition[
   x^2 >= (supportX/2)^2 && y^2 >= (barY/2)^2 && 
    z^2 <= (supportZ/2)^2, tvars, pars];

{evals, evecs} = 
  NDEigensystem[{eigenmodeOperator == {0, 0, 0}, constraintB}, 
   tvars[[1]], t, {x, y, z} \[Element] mesh, 6];

eigenfrequencies = Re[Sqrt[evals]/(2 Pi)]