1
$\begingroup$

This question is partly based off Gonza's question and Mauricio's answer. I wanted to adapt the FEM method to simulate a point pressure being applied to distort an encapsulated material that is resting on a surface. I'll include an image of what I've achieved so far, as I'm probably not describing it very well: enter image description here I'm trying to apply different material properties to both regions of the generated mesh, following roughly the example listed in ElementMeshCreation here. The code uses a replacement rule so that different properties are called depending on the value of ElementMarker with:

(*If statement to replace material properties depending on ElementMarker*)
props = If[ElementMarker == 10, Propiedades10, Propiedades20];

(*NDSolve with boundary conditions included*)
{u1, v1, \[Sigma]x1, \[Sigma]y1, \[Tau]xy1} = NDSolveValue[{PS == bcN, hl, bccorrect} /. props, {u, v, \[Sigma]x, \[Sigma]y, \[Sigma]xy}, Element[{x, y}, Mesh4]];

but I am getting this error:

If[ElementMarker==10,Propiedades10,Propiedades20]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing

Is it the wrong place to be inserting a replacement rule? I feel that the if statement may not be referencing the correct RegionMarker value? The full code I've used/modified is:

(*Initiation*)
Needs["NDSolve`FEM`"];

(*Forces*)
q = 60;

(*Material properties*)
Propiedades10 = {Y -> 25000000, \[Nu] -> 47/100};
Propiedades20 = {Y -> 105000000000, \[Nu] -> 30/100};

(*Geometry*)
(*https://reference.wolfram.com/language/FEMDocumentation/tutorial/\
ElementMeshCreation.html*)
L = 3;h1 = 0.4;sxMin = 0.1 L; sxMax = 0.9 L; syMin = 2 h1/5; syMax = 3 h1/5;

(*Mesh*)
Reg4 = ToBoundaryMesh["Coordinates" -> {{0., 0.}, {L, 0.}, {L, h1}, {0, h1}, {sxMin,syMin}, {sxMax, syMin}, {sxMax, syMax}, {sxMin, syMax}},"BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4,1}, {5, 6}, {6, 7}, {7, 8}, {8, 5}}]}];

Mesh4 = ToElementMesh[Reg4,"RegionMarker" -> {{{L/2, h1/2}, 10},{{sxMin/2, syMin/2}, 20}},MeshQualityGoal -> 0];

(*FEM Solution*)
(*2D Hooke's law*)
hl = {\[Sigma]x[x, y] ==Y/(1 - \[Nu]^2) (D[u[x, y], x] + \[Nu] D[v[x, y], y]), \[Sigma]y[x, y] == Y/(1 - \[Nu]^2) (D[v[x, y], y] + \[Nu] D[u[x, y], x]), \[Sigma]xy[x, y] == (Y*\[Nu])/(1 - \[Nu]^2) (D[u[x, y], x] + D[v[x, y], y])};

(*Equations*)
PS = {Inactive[Div][{{0, -((Y \[Nu])/(1 - \[Nu]^2))}, {-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))), 0}}.Inactive[Grad][v[x, y], {x, y}], {x, y}] + Inactive[Div][{{-(Y/(1 - \[Nu]^2)), 0}, {0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}}.Inactive[Grad][u[x, y], {x, y}], {x, y}], Inactive[Div][{{0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}, {-((Y \[Nu])/(1 - \[Nu]^2)),0}}.Inactive[Grad][u[x, y], {x, y}], {x, y}] + Inactive[Div][{{-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))), 0}, {0, -(Y/(1 - \[Nu]^2))}}.Inactive[Grad][v[x, y], {x, y}], {x, y}]};

(*BCs*)
(*Neumann - Boundary condition of force exerted*)
bcN = {0, NeumannValue[-q, {1.25 <= x <= 1.75, y == h1}]};

(*Dirichlet - Boundary condition of held positions*)
bccorrect = {DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x <= 3 && y <= 0.01],DirichletCondition[v[x, y] == 0, x == L && y <= 0.01]};

(*If statement to replace material properties depending on ElementMarker*)
props = If[ElementMarker == 10, Propiedades10, Propiedades20];

(*NDSolve with boundary conditions included*)
{u1, v1, \[Sigma]x1, \[Sigma]y1, \[Tau]xy1} = NDSolveValue[{PS == bcN, hl, bccorrect} /. props, {u, v, \[Sigma]x, \[Sigma]y, \[Sigma]xy}, Element[{x, y}, Mesh4]];

(*Deformation*)
DMesh1 = ElementMeshDeformation[Mesh4, {u1, v1}, "ScalingFactor" -> 6*10^4];

(*Show deformed mesh in red and green*)
Show[{DMesh1["Wireframe"["MeshElementStyle" ->{Directive[FaceForm[Green]],Directive[FaceForm[Red]]}]]}, ImageSize -> Large]

Thanks!

*Edited code for clarity

*Edit 2: Line breaks in copy-pasted code were breaking it. Fixed this issue.

*Edit 3: Question wasn't clear. Edited.

$\endgroup$
3
$\begingroup$

I made a few changes to the boundary conditions and made the parameters rules, I also removed hooks law as it was not necessary to show the issue and it did not have boundary conditions:

(*Dirichlet-Boundary condition of held positions*)
bccorrect = {DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x <= 3 && y == 0], 
DirichletCondition[v[x, y] == 0, x == L || x == 0]};

(*If statement to replace material properties depending on ElementMarker*)
pY = Y -> If[ElementMarker == 10, 25000000, 105000000000];
pNU = \[Nu] -> If[ElementMarker == 10, 46/100, 3/10];

(*NDSolve with boundary conditions included*)
{u1, v1} = NDSolveValue[{PS == bcN, bccorrect} /. {pY, pNU}, {u, v}, 
   Element[{x, y}, Mesh4]];

(*Deformation*)
DMesh1 = ElementMeshDeformation[Mesh4, {u1, v1},"ScalingFactor" -> 2.5*10^5];

(*Show deformed mesh in red and green*)
Show[{DMesh1[
   "Wireframe"[
    "MeshElementStyle" -> {Directive[FaceForm[Green]], 
      Directive[FaceForm[Red]]}]]}, ImageSize -> Large]

enter image description here

| improve this answer | |
$\endgroup$

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.