# Double shell wall contact with AceGen

I have a shell with double walls and some space between (Picture 1). The 2 walls are connected with three edges. If I apply force to outer wall and if the force is big enough, outer shell wall will drill through the inner shell wall (Picture 2). I would like to avoid that using contact elements. How would I do that?

Here is the code to generate problem (my attempt to add contact elements is in the comments):

<< NDSolveFEM

mesh1 = ToBoundaryMesh["Coordinates" -> {
{-1, -1, 0}, {-0.5, -1, 0}, {0, -1, 0}, {0.5, -1, 0}, {1, -1, 0},
{-1, 0, 0}, {0, 0, 0}, {1, 0, 0},
{-1, 1, 0}, {-0.5, 1, 0}, {0, 1, 0}, {0.5, 1, 0}, {1, 1, 0}},
"BoundaryElements" -> {QuadElement[{{1, 3, 11, 9, 2, 7, 10, 6},
{3, 5, 13, 11, 4, 8, 12, 7}}]}
];
mesh2 = ToBoundaryMesh["Coordinates" -> {
{-1, -1, 0}, {-0.5, -1, 0}, {0, -1, 0}, {0.5, -1, 0}, {1, -1, 0},
{-1, 0, 0}, {0, 0, 0.141}, {1, 0, 0},
{-1, 1, 0}, {-0.5, 1, 0.1}, {0, 1, 0.141}, {0.5, 1, 0.1}, {1, 1, 0}},
"BoundaryElements" -> {QuadElement[{{1, 3, 11, 9, 2, 7, 10, 6},
{3, 5, 13, 11, 4, 8, 12, 7}}]}
];

<< AceFEM

SMTInputData[];
{"OuterShell", "ML:SEMSS2SDFHYS2SP6DFiSVenant",
{"E *" -> 10^6, "t *" -> 10^-3, "\[Nu] *" -> 0.3}},
{"InnerShell", "ML:SEMSS2SDFHYS2SP6DFiSVenant",
{"E *" -> 10^7, "t *" -> 10^-3, "\[Nu] *" -> 0.3}},
{"Contact", "ExamplesCTD3V1S1DN2", {"\[Rho] *" -> 1}}
];

SMTAddMesh[mesh1, "InnerShell"(*, "BodyID" -> "B1", "BoundaryDomainID" -> "Contact"*)];
SMTAddMesh[mesh2, "OuterShell"(*, "BodyID" -> "B2", "BoundaryDomainID" -> "Contact"*)];

SMTAddEssentialBoundary[Point[{-1, -1, 0}, "D"], 1 -> 0, 2 -> 0, 3 -> 0];
SMTAddEssentialBoundary[Point[{1, -1, 0}, "D"], 1 -> 0, 2 -> 0, 3 -> 0];
SMTAddEssentialBoundary[Point[{1, 1, 0}, "D"], 1 -> 0, 2 -> 0, 3 -> 0];
SMTAddEssentialBoundary[Point[{-1, 1, 0}, "D"], 1 -> 0, 2 -> 0, 3 -> 0];

SMTAddNaturalBoundary[{"Y" == 1 && "X" == 0 & && "D" && {"DomainID", "OuterShell"}, 3 -> -8}];

SMTAnalysis[];

nstep = 10;
\[CapitalDelta]\[Lambda] = 1/nstep;
tolNR = 10^-8;
maxNR = 15;

Do[
SMTNextStep["\[CapitalDelta]\[Lambda]" -> \[CapitalDelta]\[Lambda]];
While[
SMTConvergence[tolNR, maxNR],
SMTNewtonIteration[]
],
{i, 1, nstep}
]

• Which version of AceFEM are you using? This is the beginning of the output showing after running SMTAnalysis[]: Direct definition of input data (SMTNodes=..., SMTElements=...) is not supported anymore. That is old, undocumented input data format from the time before the first release of AceFEM in year 2006. Please use only documented commands (SMTAddNode, SMTAddMesh, SMTMesh, SMTAddElement, SMTAddEssentialBoundary SMTAddNaturalBoundary, SMTAddInitialBoundary, SMTAddSensitivity).  Aug 13 '19 at 11:31
• I am using version 6.923 for MacOS. I load element code using AceShare -> Main Library. Everything else is default I believe. I don't get any errors. Aug 13 '19 at 12:28

The contact search in AceFEM only covers 3D bodies, not surfaces. In your case, you probably don't need full contact formulation. Instead, you can insert "gap" elements between the surfaces.

To @JozeK's suggestion I managed to write gap element similar to the one in this question: Prescribe the rotation as the boundary condition and the plane section constraint in AceFEM

I post the code for element here if anyone else will find it helpful.

<<AceGen

SMSInitialize["GapElement", "Environment" -> "AceFEM"];

SMSTemplate[
"SMSTopology" -> "C1" ,
"SMSNoNodes" -> 3,
"SMSDOFGlobal" -> {3, 3, 1},
"SMSDefaultIntegrationCode" -> 0,
"SMSNodeID" -> {"D", "D", "Lagrange -LP -L"},
"SMSDomainDataNames" -> {"k -spring constant"},
"SMSDefaultData" -> {10^6}
];

SMSStandardModule["Tangent and residual"];

X ⊢ SMSIO["All coordinates"][[{1, 2}]];
u ⊢ SMSIO["All DOFs"][[{1, 2}]];
λ ⊢ SMSIO["All DOFs"][[3, 1]];
\[DoubleStruckP] = SMSIO["Nodal DOFs"] // Flatten;
{k} ⊢ SMSIO["All domain data"];

x ⊨ X + u;
\[DoubleStruckN] ⊨ X[[2]] - X[[1]];
\[DoubleStruckT] ⊨ x[[2]] - x[[1]];
gap ⊨ \[DoubleStruckT].\[DoubleStruckN] * \[DoubleStruckT].\[DoubleStruckT];

constraint = λ + k gap;
Πconstraint ⊨ SMSIf[
constraint < -10^-8,
(λ + k/2 gap) gap,
-(1/(2 k)) λ^2
];

\[DoubleStruckCapitalR] ⊨ SMSD[Πconstraint, \[DoubleStruckP]];
\[DoubleStruckCapitalK] = SMSD[\[DoubleStruckCapitalR], \[DoubleStruckP]];

SMSWrite[];
SMTMakeDll[];


In the example in the question I implemented this element in the code this way

<< NDSolveFEM

crdsTop = {{-1, -1, 0}, {-0.5, -1, 0}, {0, -1, 0}, {0.5, -1, 0},
{1, -1, 0}, {-1, 0, 0}, {0, 0, 0}, {1, 0, 0}, {-1, 1, 0},
{-0.5, 1, 0}, {0, 1, 0}, {0.5, 1, 0}, {1, 1, 0}};
crdsBottom = {{-1, -1, 0}, {-0.5, -1, 0}, {0, -1, 0}, {0.5, -1, 0},
{1, -1, 0}, {-1, 0, 0}, {0, 0, 0.141}, {1, 0, 0}, {-1, 1, 0},
{-0.5, 1, 0.1}, {0, 1, 0.141}, {0.5, 1, 0.1}, {1, 1, 0}};
mesh1 = ToBoundaryMesh[
"Coordinates" -> crdsTop,
"BoundaryElements" -> {QuadElement[{{1, 3, 11, 9, 2, 7, 10, 6}, {3, 5, 13, 11, 4, 8, 12, 7}}]}
];
mesh2 = ToBoundaryMesh[
"Coordinates" -> crdsBottom,
"BoundaryElements" -> {QuadElement[{{1, 3, 11, 9, 2, 7, 10, 6},{3, 5, 13, 11, 4, 8, 12, 7}}]}
];
crdsGap = Flatten[DeleteCases[
MapThread[List, {crdsTop, crdsBottom}], {t_, b_} /; t == b
], 1];
eleGap = Partition[Range[Length@crdsGap], 2];

<< AceFEM

SMTInputData[];
{"OuterShell", "ML:SEMSS2SDFHYS2SP6DFiSVenant",{"E *" -> 10^6, "t *" -> 10^-3, "ν *" -> 0.3}},
{"InnerShell", "ML:SEMSS2SDFHYS2SP6DFiSVenant",{"E *" -> 10^7, "t *" -> 10^-3, "ν *" -> 0.3}},
{"Gap", "GapElement", {"k *" -> 10^6}}
];

SMTAddEssentialBoundary[Point[{-1, -1, 0}, "D"], 1 -> 0, 2 -> 0, 3 -> 0];
SMTAddEssentialBoundary[Point[{1, -1, 0}, "D"], 1 -> 0, 2 -> 0, 3 -> 0];
SMTAddEssentialBoundary[Point[{1, 1, 0}, "D"], 1 -> 0, 2 -> 0, 3 -> 0];
SMTAddEssentialBoundary[Point[{-1, 1, 0}, "D"], 1 -> 0, 2 -> 0, 3 -> 0];
SMTAddNaturalBoundary[{"Y" == 1 && "X" == 0 & && "D" && {"DomainID", "OuterShell"}, 3 -> -9}];

SMTAnalysis[];

nstep = 10;
Δλ = 1/nstep;
tolNR = 10^-6;
maxNR = 15;

Do[
SMTNextStep["Δλ" -> Δλ];
While[
SMTConvergence[tolNR, maxNR],
SMTNewtonIteration[]
],
{i, 1, nstep}
]
`