Prescribe the rotation as the boundary condition and the plane section constraint in AceFEM

Question

I am going to analyze a rectangular beam under bending moment in AceFEM. To prescribe the bending I would like to prescribe the rotation angle of the section at the point $z=L_z$ (see the AceFEM file below). In addition, I would like to constrain this section to remain plane after bending. Theoretically, I know how to do it but in AceFEM I need help.

$\textbf{In theory}$:

Supposing that the prescribed rotation angle of the section at $z=L_z$ is $\theta_s$, the following multi-point constraint imposes the plane section remain plane hypothesis on this section, i.e.,

$\tan(\theta_s)=\frac{z_0-z_i}{x_0-x_i}$, $\qquad$ for each node on this section.

where $x_0$ and $z_0$ are the bottom nodes of this section.

The above constraint can be easily imposed by writing a quasi-potential that makes $\tan(\theta_s)-\frac{z_0-z_i}{x_0-x_i}$ equal to zero. In other words,

$\Pi = \int_s \lambda (\tan(\theta_s)-\frac{z_0-z_i}{x_0-x_i}) \rm{d}s$

where $\lambda$ is a Lagrangian multiplier that makes the above equation zero at each node of the section.

$\textbf{In AceFEM}$:

I have made a simple example in AceFEM of a rectangular beam. The beam is clamped at the section located at $z=0$. While the rotation should be prescribed the section located at $z=L_z$.

<< AceFEM`

{Lx, Ly, Lz} = {5, 5, 50};
{Nx, Ny, Nz} = {10, 10, 100};
points = {{0, 0, 0}, {Lx, 0, 0}, {Lx, Ly, 0}, {0, Ly, 0}, {0, 0, 
    Lz}, {Lx, 0, Lz}, {Lx, Ly, Lz}, {0, Ly, Lz}};

setup1[] := (
  SMTInputData[];
  SMTAddDomain[{"A","OL:SED3H2DFLEH2Hooke", {"E *" -> 35000, "ν *" -> 0.1}}];
  SMTAddMesh[Hexahedron[points], "A", "H2", {Nx, Ny, Nz}];
  SMTAddEssentialBoundary[
   Polygon[ {{0, 0, 0}, {Lx, 0, 0}, {Lx, Ly, 0}, {0, Ly, 0}}, "D"], 
   1 -> 0, 2 -> 0, 3 -> 0]; SMTAnalysis["Output" -> "Example.out"];
  )

setup1[];

Now I am stuck at this point. Because I don't know how to impose the above constraints discussed. Of course, it's not mandatory for me to use the method discussed and any other idea is welcomed.

Thanks in advance

Answer 1

After a fun bicycle ride today I had an Idea how to simply enforce the constraint I initially had in mind.

Say you prescribe the rotation of the cross section in terms of the normal on this cross section. The rotation of such a normal can easily described in terms of an Euler angle, at least on a 2D rotation.

Then the constraint may state that the in plane vector of each node with respect to a reference node in the plane is perpendicular to the prescribed surface normal and thus their scalar product vanishes.

My idea was to implement this concept within line elements, which avoids an over constraint situation. The finite element is very simple and reads:

<< AceGen`;
SMSInitialize["C1Boundary", "Environment" -> "AceFEM", 
  "Mode" -> "Prototype"];
AdditionalGraphics = Function[{e, m, b, x}, {
    If[b,
     {{Red, Arrow[{
         SMTNodeData[e[[3, 1]], "X"] + SMTNodeData[e[[3, 1]], "at"],
         SMTNodeData[e[[3, 1]], "X"] + SMTNodeData[e[[3, 1]], "at"] + 
          SMTDomainData[e[[2]], "Data"]
         }]
       }}
     ]
    }];
SMSTemplate[
  "SMSTopology" -> "C1",
  "SMSNoNodes" -> 3,
  "SMSDOFGlobal" -> {3, 3, 1},
  "SMSNodeID" -> {"D", "D", "Lagrange -LP -L"},
  "SMSAdditionalNodes" -> Hold[{Null} &],
  "SMSDefaultIntegrationCode" -> 0,
  "SMSAdditionalGraphics" -> AdditionalGraphics,
  "SMSDomainDataNames" -> {"nx", "ny", "nz"},
  "SMSDefaultData" -> {1, 0, 0}
  ];
SMSStandardModule["Tangent and residual"];
XI ⊢ 
  Table[SMSReal[nd$$[i, "X", j]], {i, SMSNoNodes}, {j, 
    SMSNoDimensions}];
uI ⊢ 
  SMSReal[Table[
    nd$$[i, "at", j], {i, SMSNoNodes}, {j, SMSDOFGlobal[[i]]}]];
DOFVector ⊨ Flatten[uI];
\[DoubleStruckN] ⊢ 
  SMSReal[{es$$["Data", 1], es$$["Data", 2], es$$["Data", 3]}];
λLagrange ⊨ uI[[3, 1]];
xI ⊨ XI[[;; 2]] + uI[[;; 2]];
\[DoubleStruckT] ⊨ xI[[1]] - xI[[2]];
constraint ⊨ \[DoubleStruckT].\[DoubleStruckN];
Πconstraint ⊨ λLagrange constraint;

δΠ ⊨ 
  SMSD[Πconstraint, DOFVector, "Constant" -> {}];
 ΔδΠ ⊨ 
  SMSD[δΠ, DOFVector];
SMSExport[δΠ , p$$, "AddIn" -> True];
SMSExport[ΔδΠ, s$$, "AddIn" -> True];

SMSWrite[];

The only issue here is that you need to know all the nodes on your cross section at the mesh input step. However this can be addresses for example as shown below:

<< AceFEM`
{Lx, Ly, Lz} = {5, 5, 80};
{Nx, Ny, Nz} = {4, 4, 10};
points = {{0, 0, 0}, {Lx, 0, 0}, {Lx, Ly, 0}, {0, Ly, 0}, {0, 0, 
    Lz}, {Lx, 0, Lz}, {Lx, Ly, Lz}, {0, Ly, Lz}};
(*Generate Solid Mesh to get End-Section Nodes*)
SMTInputData[];
SMTAddDomain[{"A","OL:SED3H2DFHYH2NeoHooke", {"E *" -> 35000, "ν *" -> 0.}}];
SMTAddMesh[Hexahedron[points], "A", "H2", {Nx, Ny, Nz}];
SMTAnalysis[];
(*Find Nodes on End-Section*)

EndSectionCoors = SMTNodeData[SMTFindNodes["Z" == Lz &], "X"];
SectionMesh = 
  Table[{1, i + 1}, {i, (EndSectionCoors // Length) - 1}];
(*Start Acual Problem*)
SMTInputData[];
SMTAddDomain[{"A", 
   "OL:SED3H2DFHYH2NeoHooke", {"E *" -> 35000, "ν *" -> 0.}}];
SMTAddDomain[{"B", "C1Boundary", {}}];
SMTAddMesh[Hexahedron[points], "A", "H2", {Nx, Ny, Nz}];
SMTAddMesh[EndSectionCoors, {"B" -> SectionMesh}];
SMTAddEssentialBoundary[Point[{2.5, 2.5, 0}, "D"], 1 -> 0, 2 -> 0, 
  3 -> 0];
SMTAddEssentialBoundary["Z" == 0 &, 3 -> 0];
SMTAnalysis[];
SMTDomainData["B", "Data", {0, 0, 1}]; θ = 0;

The load steering in this case is done by modification of the section normal, which is a DomainData.

Do[
 SMTNextStep[1, 1];
 Δθ = π/100 + 10^-8;
 θ = θ + Δθ;
 n = SMTDomainData["B", "Data"];
 n = n.RotationMatrix[Δθ, {0, 1, 0}];
 SMTDomainData["B", "Data", n];
 While[SMTConvergence[10^-6, 20], SMTNewtonIteration[]];
 θ (180/π) // N
 , {step, 100}]
SMTDomainData["B", "Data", 10 n];
Row[{
  SMTShowMesh["DeformedMesh" -> True, "BoundaryConditions" -> True],
  SMTShowMesh["DeformedMesh" -> True, "BoundaryConditions" -> True, 
   ViewPoint -> Bottom]
  }]

Let me know what you think about this.

Answer 2

First of all, let me clarify the boundary condition that I mentioned. In the answer that you posted you didn't impose the rotation of the section. While I was searching for a boundary condition in which one can impose the rotation of the section while it remains plane. The idea of using a rigid element to impose plane section is not good because the section should just remain plane in one plane while in the perpendicular view it can expand or shrink. This is illustrated in the image below. Using a rigid element will make it plane (rigid) in all directions.

Regarding the rotation of the section, I searched the documentation of AceFEM and I found an example (the torsion test example) that shows how to impose the rotation on a section. Thus I used the same idea to impose bending moment. In this case, since you manually prescribe the displacement of the nodes in the yz plane to follow a rotational path, the section will remain plane in this plane, while in the other direction (x) it may deform.

<< AceFEM`
SMTInputData[];
SMTAddDomain[{"A", 
   "OL:SED3H1DFHYH1NeoHooke" , {"E *" -> 2500, "ν *" -> 0.35}}];
L = 20; a = 2; b = 1; ϕM = π/3 // N;
SMTAddMesh[
  Raster3D[{{{{-b, -a, 0}, {b, -a, 0}}, {{-b, a, 0}, {b, a, 
       0}}}, {{{-b, -a, L}, {b, -a, L}}, {{-b, a, L}, {b, a, L}}}}], 
  "A", "H1", {8, 8, 20}];
SMTAddEssentialBoundary[
  Polygon[{{-b, -a, 0}, {b, -a, 0}, {b, a, 0}, {-b, a, 0}}], 1 -> 0, 
  2 -> 0, 3 -> 0];
SMTAddEssentialBoundary[
  Polygon[{{-b, -a, L}, {b, -a, L}, {b, a, L}, {-b, a, L}}], 2 -> 0, 
  3 -> 0];
SMTAnalysis[];
SMTNextStep["λ" -> π/300.];
Mϕ = {{0, 0}};
ExtremeSurface = SMTFindNodes["Z" == L &];
Bp = SMTNodeData[ExtremeSurface, "Bp"];
While[
 ϕ = SMTRData["Multiplier"]; Δϕ = 
  SMTRData["MultiplierIncrement"];
 Bt = Map[{0, #[[2]] (Cos[ϕ] - 1) - #[[3]] Sin[ϕ], #[[
        3]] (Cos[ϕ] - 1) + #[[2]] Sin[ϕ]} &, 
   SMTNodeData[ExtremeSurface, "X"]];
 SMTNodeData[ExtremeSurface, "dB", (Bt - Bp)/Δϕ];
 While[step = 
   SMTConvergence[10^-8, 
    16, {"Adaptive BC", 12, 0.0001, ϕM/100, ϕM}],
  SMTNewtonIteration[];];
 If[step[[4]] == "MinBound" , SMTStatusReport["Analyze"]; 
  SMTStepBack[];];
 If[Not[step[[1]]], Bp = Bt;
  Mi = SMTResidual[ExtremeSurface][[;; , 3]];
  Zi = SMTNodeData[ExtremeSurface , "X"][[;; , 2]] // Abs;
  M = 2*Mi.Zi;
  AppendTo[Mϕ, {ϕ, M}];
  SMTShowMesh["DeformedMesh" -> True, "Show" -> "Window" , 
   "Field" -> "Mises stress", ViewVertical -> {1, 0, 0}];
  ];
 step[[3]]
 , If[step[[1]], SMTStepBack[];];
 SMTNextStep["Δλ" -> step[[2]]];
 ]
GraphicsRow[{SMTShowMesh["BoundaryConditions" -> True , 
   ViewPoint -> Left, "DeformedMesh" -> True , ImageSize -> 500] , 
  SMTShowMesh["BoundaryConditions" -> True , ViewPoint -> Top , 
   "DeformedMesh" -> True, ImageSize -> 2000]}]

Please comment your idea and points about this code. In my opinion this code is doing the job fine.

Answer 3

i can tell the task your describing is not trivial. Just recently I was facing the same problem. Your description of the constraint is not completely clear to me, although one way to implement such would be a boundary element. However, it turned out that for me it was sufficient to model the effect of a bending moment on a beam, by application of a follower load distribution on the one end. Additionally I used an approximately rigid solid mesh to enforce the plane cross section at the end. Maybe such is also suitable in your case, take a look at the code below. First the boundary element:

<< AceGen`;
SMSInitialize["S2SurfaceLoad", "Environment" -> "AceFEM"];
SMSTemplate[
  "SMSTopology" -> "S2"
  , "SMSSymmetricTangent" -> False
  , "SMSDomainDataNames" -> {
    "tx -traction in global X direction"
    , "ty -traction load in global Y direction"
    , "tz -traction load in global Z direction"}
  , "SMSDefaultData" -> {0, 0, 0}
  , "SMSDefaultIntegrationCode" -> 3];
SMSStandardModule["Tangent and residual"];
SMSDo[Ig, 1, SMSInteger[es$$["id", "NoIntPoints"]]];
{ξ, η, ζ, ω} ⊢ 
  Array[SMSReal[es$$["IntPoints", #1, Ig]] &, 4];
XI ⊢ 
  Table[SMSReal[nd$$[i, "X", j]], {i, SMSNoNodes}, {j, 
    SMSNoDimensions}];
uI ⊢ 
  SMSReal[Table[
    nd$$[i, "at", j], {i, SMSNoNodes}, {j, SMSDOFGlobal[[i]]}]];
DOFVector ⊨ Flatten[uI];
\[DoubleStruckT] ⊨ 
  SMSReal[{es$$["Data", 1], es$$["Data", 2], es$$["Data", 3]}];
λload ⊨ SMSReal[rdata$$["Multiplier"]];

SHP ⊨ {1/4 (-1 + η) η (-1 + ξ) ξ, 
   1/4 (-1 + η) η ξ (1 + ξ), 
   1/4 η (1 + η) ξ (1 + ξ), 
   1/4 η (1 + η) (-1 + ξ) ξ, -(1/
     2) (-1 + η) η (-1 + ξ^2), -(1/
     2) (-1 + η^2) ξ (1 + ξ), -(1/
     2) η (1 + η) (-1 + ξ^2), -(1/
     2) (-1 + η^2) (-1 + ξ) ξ, (-1 + η^2) (-1 + \
ξ^2)};
SMSFreeze[X, SHP.XI];

rξ ⊨ SMSD[X, ξ]; rη ⊨ 
 SMSD[X, η]; rζ ⊨ rξ\[Cross]rη;
e1r ⊨ 
 rξ/SMSSqrt[rξ.rξ]; e2r ⊨ 
 rη/SMSSqrt[rη.rη]; e3r ⊨ 
 rζ/SMSSqrt[rζ.rζ];
detJ ⊨ SMSSqrt[rζ.rζ];

u ⊨ SHP.uI;
x ⊨ X + u;
aξ ⊨ SMSD[x, ξ]; aη ⊨ 
 SMSD[x, η]; aζ ⊨ aξ\[Cross]aη;
e1 ⊨ 
 aξ/SMSSqrt[aξ.aξ]; e2 ⊨ 
 aη/SMSSqrt[aη.aη]; e3 ⊨ 
 aζ/SMSSqrt[aζ.aζ];
Qr ⊨ 
  e1\[TensorProduct]e1r + e2\[TensorProduct]e2r + 
   e3\[TensorProduct]e3r;
force ⊨ (Qr.\[DoubleStruckT]);
Π ⊨ λload force.u;

δΠ ⊨ 
  SMSD[Π, DOFVector, "Constant" -> {force}];
 ΔδΠ ⊨ 
  SMSD[δΠ, DOFVector];
SMSExport[detJ ω δΠ , p$$, "AddIn" -> True];
SMSExport[detJ ω ΔδΠ, s$$, 
  "AddIn" -> True];

SMSEndDo[];

SMSWrite[];

And here the modified problem concerning the boundary mesh and the rigid solid on the end. I also switched to an element more suitable for finite rotations.

 << AceFEM`

{Lx, Ly, Lz} = {5, 5, 100};
{Nx, Ny, Nz} = {4, 4, 20};
points = {{0, 0, 0}, {Lx, 0, 0}, {Lx, Ly, 0}, {0, Ly, 0}, {0, 0, 
    Lz}, {Lx, 0, Lz}, {Lx, Ly, Lz}, {0, Ly, Lz}};

SMTInputData[];
SMTAddDomain[{"A", 
   "OL:SED3H2DFHYH2NeoHooke", {"E *" -> 35000, "\[Nu] *" -> 0.}}];
SMTAddDomain[{"Rigid", 
   "OL:SED3H2DFHYH2NeoHooke", {"E *" -> 10^10, "\[Nu] *" -> 0}}];
SMTAddDomain[{"B+", "S2SurfaceLoad", {"tz *" -> 10^3}}];
SMTAddDomain[{"B-", "S2SurfaceLoad", {"tz *" -> -10^3}}];
SMTAddMesh[Hexahedron[points], "A", "H2", {Nx, Ny, Nz}];
SMTAddMesh[
  Hexahedron[
   Join[# + {0, 0, Lz} & /@ points[[;; 4]], # + {0, 0, 10^0} & /@ 
     points[[5 ;;]]]], "A", "H2", {Nx, Ny, Nz}];
SMTAddMesh[
  Polygon[{{0, 0, Lz}, {Lx/2, 0, Lz}, {Lx/2, Ly, Lz}, {0, Ly, Lz}}], 
  "B+", "S2", {Nx, Ny}/2];
SMTAddMesh[
  Polygon[{{Lx/2, 0, Lz}, {Lx, 0, Lz}, {Lx, Ly, Lz}, {Lx/2, Ly, Lz}}],
   "B-", "S2", {Nx, Ny}/2];
SMTAddEssentialBoundary[Polygon[points[[;; 4]], "D"], 1 -> 0, 2 -> 0, 
 3 -> 0]; SMTAnalysis[];
Monitor[
 Do[
  SMTNextStep[.005, .005];
  While[SMTConvergence[], SMTNewtonIteration[]];
  post = SMTShowMesh["DeformedMesh" -> True, ViewPoint -> Front];
  , {step, 2000}]
 , post]
post

Answer 4

So after a wide range of trial and error and non-convergence issues, I developed the following element that suitably satisfies my boundary conditions.

<< AceGen`

SMSInitialize["BoundaryConstraint" , "Environment" -> "AceFEM", 
  "Mode" -> "Optimal"];
nNodes = 9; nLNodes = 8;
SMSTemplate["SMSTopology" -> "S2", "SMSSymmetricTangent" -> True, 
  "SMSNoNodes" -> 
   nNodes + 
    nLNodes(*9 nodes for displacements plus 8 nodes for the \
lagrangian multiplier (without node number 1)*), 
  "SMSDefaultIntegrationCode" -> 0, 
  "SMSDOFGlobal" -> Join[Table[3, nNodes], Table[1, nLNodes]], 
  "SMSAdditionalNodes" -> "{#2,#3,#4,#5,#6,#7,#8,#9}&", 
  "SMSNodeID" -> Join[Table["D", nNodes], Table["LR", nLNodes]]];

SMSStandardModule["Tangent and residual"];

Xi \[DoubleRightTee] Array[SMSReal[nd$$[#, "X", 1]] &, nNodes];
Yi \[DoubleRightTee] Array[SMSReal[nd$$[#, "X", 2]] &, nNodes];
Zi \[DoubleRightTee] Array[SMSReal[nd$$[#, "X", 3]] &, nNodes];

ui \[DoubleRightTee] 
  SMSReal[Array[nd$$[#1, "at", 1] &, nNodes], -0.1, 0.1];
vi \[DoubleRightTee] 
  SMSReal[Array[nd$$[#1, "at", 2] &, nNodes], -0.1, 0.1];
wi \[DoubleRightTee] 
  SMSReal[Array[nd$$[#1, "at", 3] &, nNodes], -0.1, 0.1];

\[Lambda]i \[DoubleRightTee] 
  SMSReal[Array[nd$$[nNodes + #1, "at", 1] &, nLNodes]];

\[DoubleStruckA]t \[DoubleRightTee] 
  Join[Flatten[Transpose[{ui, vi, wi}]], 
   Flatten[Transpose[{\[Lambda]i}]]];

xi \[DoubleRightTee] Xi + ui; 
yi \[DoubleRightTee] Yi + vi;
zi \[DoubleRightTee] Zi + wi;

multiplier \[DoubleRightTee] SMSReal[rdata$$["Multiplier"]];
\[Theta]t \[DoubleRightTee] multiplier*(\[Pi]/2) // N;
\[Rho] = 0;

gi \[DoubleRightTee] 
  MapThread[(xi[[1]] - #2) Cos[\[Theta]t] - 
     Sin[\[Theta]t] (zi[[1]] - #1) &, {zi[[2 ;; 9]], xi[[2 ;; 9]]}];

\[CapitalPi] \[DoubleRightTee] \[Lambda]i.gi + 1/2 \[Rho] gi.gi;

SMSDo[i, 1, Length[\[DoubleStruckA]t]];
    Rg \[DoubleRightTee] SMSD[\[CapitalPi], \[DoubleStruckA]t, i];
    SMSExport[Rg, p$$[i], "AddIn" -> True];
	SMSDo[j, 1, Length[\[DoubleStruckA]t]];
		Kg \[DoubleRightTee] SMSD[Rg, \[DoubleStruckA]t, j];
		SMSExport[Kg, s$$[i, j], "AddIn" -> True];
    SMSEndDo[];
SMSEndDo[];

SMSWrite[];
SMTMakeDll[]; 

Note that since the tangent trig function is not defined in 90 degree, I have changed it to sin and cos so that the above problem will be avoided. The simulation is as following:

<< AceFEM`
SMTInputData[];
\[Phi]M = 2;
L = 80; y = 4; z = 4;
SMTAddDomain[{"A", 
   "OL:SED3H2DFHYH2NeoHooke" , {"E *" -> 25000, "\[Nu] *" -> 0.35}}];
SMTAddDomain[{"B", "BoundaryConstraint", {}}];
{Nx, Ny, Nz} = {10, 10, 80};
SMTAddMesh[
  Raster3D[{{{{0, -y, -z}, {0, y, -z}}, {{0, -y, z}, {0, y, 
       z}}}, {{{L, -y, -z}, {L, y, -z}}, {{L, -y, z}, {L, y, z}}}}], 
  "A", "H1", {Nx, Ny, Nz}];
SMTAddMesh[
  Raster[{{{L, -y, -z}, {L, y, -z}}, {{L, -y, z}, {L, y, z}}}], "B", 
  "S2", {Nx, Ny}];
SMTAddEssentialBoundary[
  Polygon[{{0, -y, -z}, {0, y, -z}, {0, y, z}, {0, -y, z}}], 1 -> 0, 
  2 -> 0, 3 -> 0];
(*SMTAddEssentialBoundary[Point[{L,0,-z},"LR"],1\[Rule]0];*)
SMTAnalysis["Output" -> "test.out"];

SMTNextStep["\[Lambda]" -> \[Phi]M/50.];
M\[Phi] = {{0, 0}};

While[
 While[step = 
   SMTConvergence[10^-8, 
    16, {"Adaptive BC", 12, 0.000000001, \[Phi]M/10, \[Phi]M}],
  SMTNewtonIteration[];];
 If[step[[4]] == "MinBound" , SMTStatusReport["Analyze"]; 
  SMTStepBack[];];
 If[Not[step[[1]]],
  M0 = SMTResidual["X" == 0 &];
  M1 = Max[Length /@ M0];
  M2 = Select[M0, Length[#] == M1 &];
  Mi = M2[[;; , 3]];
  Zi = SMTNodeData[
      Intersection[SMTFindNodes["X" == 0 &], SMTFindNodes["D"]], 
      "X"][[;; , 3]] // Abs;
  M = Mi.Zi;
  \[Phi] = SMTRData["Multiplier"];
  AppendTo[M\[Phi], {\[Phi], M}];
  SMTShowMesh["DeformedMesh" -> True, "Show" -> "Window" , 
   "Field" -> "Mises stress", ViewVertical -> {1, 0, 0}];
  ];
 step[[3]]
 , If[step[[1]], SMTStepBack[];];
 SMTNextStep["\[CapitalDelta]\[Lambda]" -> step[[2]]];
 ]

And as always, your comments are welcomed.