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

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

## 4 Answers

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[], "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[] - xI[]; 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. • +1, Nice example of constrain element. – Pinti May 22 '18 at 8:12 • Thanks for your answer. Actually, I think the method that I used before does not give the correct boundary condition. The method that you have proposed here is one of the suitable ways of prescribing this boundary condition. I have thought about another way of doing it and I will post it as soon as having it ready. By the way, I think your method works. – KratosMath May 22 '18 at 20:33 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, #[] (Cos[ϕ] - 1) - #[] Sin[ϕ], #[[ 3]] (Cos[ϕ] - 1) + #[] 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[] == "MinBound" , SMTStatusReport["Analyze"]; SMTStepBack[];]; If[Not[step[]], 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[] , If[step[], SMTStepBack[];]; SMTNextStep["Δλ" -> step[]]; ] 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. • First of all, congratulations as this works out fine for you. It is worth mentioning though, that your "ExtremeSurface" may contract freely in one direction, but the other in plane direction is still rigid. However, this may cause effects which may (or may not) be neglect able, depending on your goal. To be honest, for me, as I was looking for a validation of a beam formulation with rigid cross sections as assumption, it was too. The boundary condition you prescribed at the fixed end, also restricts the in plane contraction of the front section here, just to mention it. – Sascha Maassen May 21 '18 at 11:41 • Thanks for your comment. By and large, I should examine both answers and figure out which one does the job better. – KratosMath May 21 '18 at 12:03 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 • Thanks for your answer. I will post an answer to discuss some points and issues. – KratosMath May 20 '18 at 10:06 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[] - #2) Cos[\[Theta]t] - Sin[\[Theta]t] (zi[] - #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[] == "MinBound" , SMTStatusReport["Analyze"];
SMTStepBack[];];
If[Not[step[]],
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[]
, If[step[], SMTStepBack[];];
SMTNextStep["\[CapitalDelta]\[Lambda]" -> step[]];
] And as always, your comments are welcomed.