31
$\begingroup$

Say I want to study the deformation of a pitchfork when you have it fixed on the bottom and push one side.

<< NDSolve`FEM`
Ω = 
  RegionDifference[Cuboid[{0, -5}, {1, 0}], 
   Cuboid[{0.45, -4.5}, {.55, 0}]];
Ω // DiscretizeRegion
bcs = DirichletCondition[#, {z <= -4.99}] & /@ {u[y, z] == 0., 
    v[y, z] == 0.};
mesh = ToElementMesh[Ω, "MaxCellMeasure" -> 0.005]
planeStress = {Inactive[
       Div][{{0, -((Y*ν)/(1 - ν^2))}, {-(Y*(1 - ν))/(2*(1 \
- ν^2)), 0}}.Inactive[Grad][v[y, z], {y, z}], {y, z}] + 
     Inactive[
       Div][{{-(Y/(1 - ν^2)), 
         0}, {0, -(Y*(1 - ν))/(2*(1 - ν^2))}}.Inactive[Grad][
        u[y, z], {y, z}], {y, z}], 
    Inactive[
       Div][{{0, -(Y*(1 - ν))/(2*(1 - ν^2))}, {-((Y*ν)/(1 \
- ν^2)), 0}}.Inactive[Grad][u[y, z], {y, z}], {y, z}] + 
     Inactive[
       Div][{{-(Y*(1 - ν))/(2*(1 - ν^2)), 
         0}, {0, -(Y/(1 - ν^2))}}.Inactive[Grad][
        v[y, z], {y, z}], {y, z}]} /. {Y -> 10^3, ν -> 33/100};

That's my domain

enter image description here

and this my results:

{uif, vif} = 
  NDSolveValue[{planeStress == {NeumannValue[1, y == 0 && z > -.1], 
      0}, DirichletCondition[u[y, z] == 0, z == -5], 
    DirichletCondition[v[y, z] == 0, z == -5]}, {u, 
    v}, {y, z} ∈ mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
Show[{mesh["Wireframe"], 
  dmesh["Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]

enter image description here

I was expecting that the right arm of the pitchfork would somehow know when the other one was touching it, but this is not the case at all.

I don't have much experience with FEM so I don't know if what I'm asking is impossible due to the nonlocality of the problem or what not. Is there a way to obtain the correct behaviour in mma, ie how to make the left arm collide with the right one and that they both deform due to the force applied only on the left one?

Improve this question
$\endgroup$
5
  • 6
    $\begingroup$ The topology of the solution region changes when the two arms come into contact. You may have to terminate the calculation (use WhenEvent when the this event occurs, then use the output as initial condition for a continued calculation with the two arms joined. Even this may not be quite correct, because the contact slides. Interesting question. $\endgroup$
    – bbgodfrey
    Commented Aug 19, 2016 at 18:27
  • $\begingroup$ mmmh, one could integrate till it touches for the first time, then compute the force it applies on the right arm, move the right arm, then integrate the left arm till the next contact and proceed like that till the balance of forces is reached? That could be an option to avoid resolving the time dependent problem... $\endgroup$
    – tsuresuregusa
    Commented Aug 19, 2016 at 20:19
  • $\begingroup$ The equations you have given to the Mathematica solver to discretize and solve via F.E.M don't take into account collisions/overlapping in the mesh. Adding collisions must be done via crafting boundary conditions or penalty forces which change dynamically. This is definitely non-trivial: it's quite complicated, and is normally done iteratively. Additional note: Furthermore the Mathematica solver linearizes your equations. The linearized equation is only valid in the small deformation case: distances and angles barely change (imagine for instance 10 m steel beam bending a few millimeters) ... $\endgroup$
    – darmual
    Commented Aug 7, 2018 at 11:50
  • $\begingroup$ [second part of the comment]: ... which is definitely not your case. So your results, even without the overlapping, wouldn't be valid. You could see a clear example of this by changing the Neumann value to -50. $\endgroup$
    – Kuba
    Commented Aug 7, 2018 at 13:09
  • $\begingroup$ @darmual I moved your answer to comments section since it didn't really provide an answer $\endgroup$
    – Kuba
    Commented Aug 7, 2018 at 13:10

1 Answer 1

Reset to default
15
$\begingroup$

Here it is necessary to determine not the impact, but the additional stress arising at the point of contact of two elastic elements. The simplest solution method is to determine the force that acts on each element in contact. I will show the simplest code that can be complicated to infinity.

<< NDSolve`FEM`
Ω = RegionDifference[Cuboid[{0, -5}, {1, 0}], Cuboid[{0.45, -4.5}, {.55, 0}]];

Ω // DiscretizeRegion

bcs = DirichletCondition[#, {z <= -4.99}] & /@ {u[y, z] == 0., v[y, z] == 0.};
mesh = ToElementMesh[Ω, "MaxCellMeasure" -> 0.005];
planeStress = {Inactive[
       Div][{{0, -((Y*ν)/(1 - ν^2))}, {-(Y*(1 - ν))/(2*(1 \
- ν^2)), 0}}.Inactive[Grad][v[y, z], {y, z}], {y, z}] + 
     Inactive[
       Div][{{-(Y/(1 - ν^2)), 
         0}, {0, -(Y*(1 - ν))/(2*(1 - ν^2))}}.Inactive[Grad][
        u[y, z], {y, z}], {y, z}], 
    Inactive[
       Div][{{0, -(Y*(1 - ν))/(2*(1 - ν^2))}, {-((Y*ν)/(1 \
- ν^2)), 0}}.Inactive[Grad][u[y, z], {y, z}], {y, z}] + 
     Inactive[
       Div][{{-(Y*(1 - ν))/(2*(1 - ν^2)), 
         0}, {0, -(Y/(1 - ν^2))}}.Inactive[Grad][
        v[y, z], {y, z}], {y, z}]} /. {Y -> 10^3, ν -> 33/100};
{uif, vif} = 
  NDSolveValue[{planeStress == {NeumannValue[1, y == 0 && z > -.1] + 
       NeumannValue[-1/3, y == 0.45 && z > -.1] + 
       NeumannValue[1/3, y == 0.55 && z > -.1], 0}, 
    DirichletCondition[u[y, z] == 0, z == -5], 
    DirichletCondition[v[y, z] == 0, z == -5]}, {u, 
    v}, {y, z} ∈ mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
{mesh["Wireframe"], 
 dmesh["Wireframe"[
   "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}

fig1

Note that in this solution, the elements do not crawl together, but are deformed together. Consider the general case of an arbitrary combination of parameters of force and elasticity.

<< NDSolve`FEM`
\[CapitalOmega] = 
  RegionDifference[Cuboid[{0, -5}, {1, 0}], 
   Cuboid[{0.45, -4.5}, {.55, 0}]];
\[CapitalOmega] // DiscretizeRegion;
bcs = DirichletCondition[#, {z <= -4.99}] & /@ {u[y, z] == 0., 
    v[y, z] == 0.};
mesh = ToElementMesh[\[CapitalOmega], 
  "MaxCellMeasure" -> 0.005]; f = 2;
planeStress = {Inactive[
       Div][{{0, -((Y*\[Nu])/(1 - \[Nu]^2))}, {-(Y*(1 - \[Nu]))/(2*(1 \
- \[Nu]^2)), 0}}.Inactive[Grad][v[y, z], {y, z}], {y, z}] + 
     Inactive[
       Div][{{-(Y/(1 - \[Nu]^2)), 
         0}, {0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}}.Inactive[Grad][
        u[y, z], {y, z}], {y, z}], 
    Inactive[
       Div][{{0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}, {-((Y*\[Nu])/(1 \
- \[Nu]^2)), 0}}.Inactive[Grad][u[y, z], {y, z}], {y, z}] + 
     Inactive[
       Div][{{-(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)), 
         0}, {0, -(Y/(1 - \[Nu]^2))}}.Inactive[Grad][
        v[y, z], {y, z}], {y, z}]} /. {Y -> 10^3, \[Nu] -> 33/100};
sol = ParametricNDSolveValue[{planeStress == {NeumannValue[f, 
        y == 0 && z > -.1] + NeumannValue[-g, y == 0.45 && z > -.1] + 
       NeumannValue[g, y == 0.55 && z > -.1], 0}, 
    DirichletCondition[u[y, z] == 0, z == -5], 
    DirichletCondition[v[y, z] == 0, z == -5]}, {u, 
    v}, {y, z} \[Element] mesh, {g}];


sol1 = 
  ParametricNDSolveValue[{planeStress == {NeumannValue[f, 
        y == 0 && z > -.1] + NeumannValue[-g, y == 0.45 && z > -.1] + 
       NeumannValue[g, y == 0.55 && z > -.1], 0}, 
    DirichletCondition[u[y, z] == 0, z == -5], 
    DirichletCondition[v[y, z] == 0, z == -5]}, 
   u, {y, z} \[Element] mesh, {g}];

  g0=g/.FindRoot[
 sol1[g][.45, 0] - sol1[g][.55, 0] == .1, {g, f/2}] // Quiet

Out[]= 0.834936
dmesh = ElementMeshDeformation[mesh, sol[g0], 
   "ScalingFactor" -> 1];
{mesh["Wireframe"], 
 dmesh["Wireframe"[
   "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}

fig2 For f=1, we find g0=0.335075, which is close to 1/3 found by another method.

Improve this answer
$\endgroup$
9
  • $\begingroup$ Could you explain how you did "determine the force that acts on each element in contact"? $\endgroup$
    – user21
    Commented Nov 20, 2018 at 9:39
  • $\begingroup$ @user21 I took it as a first approximation. You can further complicate the code using contact theory. $\endgroup$
    – Alex Trounev
    Commented Nov 20, 2018 at 12:55
  • $\begingroup$ How did you take this first approximation? $\endgroup$
    – user21
    Commented Nov 20, 2018 at 13:20
  • $\begingroup$ @user21 I searched for the minimum value when there is no creep and there is contact. 1/3 is close to the exact solution in this problem. $\endgroup$
    – Alex Trounev
    Commented Nov 20, 2018 at 13:41
  • 1
    $\begingroup$ This is much better. Thank you. Now it's clear how your method works. $\endgroup$
    – user21
    Commented Nov 21, 2018 at 7:02

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.