3
$\begingroup$

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?

Picture 1: Undeformed state Picture 2: Deformed state

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

<< NDSolve`FEM`

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[];
SMTAddDomain[
 {"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}
]
$\endgroup$
  • $\begingroup$ 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). $\endgroup$ – Karol Frydrych Aug 13 at 11:31
  • 2
    $\begingroup$ 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. $\endgroup$ – Jan12 Aug 13 at 12:28
3
$\begingroup$

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.

$\endgroup$
4
$\begingroup$

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,
    "SMSAdditionalNodes" -> Hold[{Null} &],
    "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]];
SMSIO[\[DoubleStruckCapitalR], "Add to", "Residual"];
\[DoubleStruckCapitalK] = SMSD[\[DoubleStruckCapitalR], \[DoubleStruckP]];
SMSIO[\[DoubleStruckCapitalK], "Add to", "Tangent"];

SMSWrite[]; 
SMTMakeDll[];

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

<< NDSolve`FEM`

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[];
SMTAddDomain[
  {"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}}
];

SMTAddMesh[mesh1, "InnerShell"];
SMTAddMesh[mesh2, "OuterShell"];
SMTAddMesh[crdsGap, {"Gap" -> eleGap}];
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}
]

enter image description here

$\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.