2
$\begingroup$

I'm trying to prestress a sample by a compressive strain of 10 percent and then stretch it by 30 percent. The way I intend to do this is apply the prestress in one step and then have the code run and strain the sample by 30 percent as usual. What I'm having issues with is that in every iteration my code both applies the prestress and the 30 percent strain, whereas I only want the prestress applied once in the begining of the code and then let the code apply the 30 percent stress using the increments outlined in the code below:

CompressionStrainX = -0.3;
CompressionStrainY = 0;
Prestress = -0.1;
AnalysisSwitch = 1; (* UNIFORM COMPRESSION Case *)

SolverSwitch = 1; (* 1 -in core; 2 -out-of-core *)
SolverMatrixType = -2;
For[AnalysisIndex = 1, AnalysisIndex <=  Length[AnalysisBatch], 
 AnalysisIndex++,

 Analysis name;
 AnalysisName = AnalysisBatch[[AnalysisIndex, 1]];
 E0StratumCorneum = AnalysisBatch[[AnalysisIndex, 2]];
 AnalysisSwitch = AnalysisBatch[[AnalysisIndex, 3]];

 Start analysis;
 SMTInputData[];

 SMTAddDomain[{
   {"Stratum-Corneum", EFormulationStratumCorneum,
    {"\[Nu]\[DoubleStruckG] *" -> \[Nu]\[DoubleStruckG], 
     "cSC*" -> cSC, "bX *" -> bX, "bY *" -> bY, "bZ *" -> bZ, 
     "\[Rho] *" -> \[Rho] , "\[Beta] *" -> \[Beta] , 
     "\[Delta] *" -> \[Delta] , "Eg" -> Eg}},
   {"Living-Epidermis", EFormulationLivingEpidermis,
    {"\[Nu]\[DoubleStruckG] *" -> \[Nu]\[DoubleStruckG], "c*" -> c, 
     "k1 *" -> k1, "k2 *" -> k2, "q *" -> q,
     "bX *" -> bX, "bY *" -> bY, "bZ *" -> bZ, "\[Rho] *" -> \[Rho] , 
     "\[Beta] *" -> \[Beta] , "\[Delta] *" -> \[Delta] , 
     "Eg" -> Eg}},
   {"Papillary-Dermis", EFormulationPapillaryDermis,
    {"\[Nu]\[DoubleStruckG] *" -> \[Nu]\[DoubleStruckG], "c*" -> c, 
     "k1 *" -> k1, "k2 *" -> k2, "q *" -> q,
     "bX *" -> bX, "bY *" -> bY, "bZ *" -> bZ, "\[Rho] *" -> \[Rho] , 
     "\[Beta] *" -> \[Beta] , "\[Delta] *" -> \[Delta] , 
     "Eg" -> Eg}},
   {"Reticular-Dermis", EFormulationReticularDermis,
    {"\[Nu]\[DoubleStruckG] *" -> \[Nu]\[DoubleStruckG], "c*" -> c, 
     "k1 *" -> k1, "k2 *" -> k2, "q *" -> q,
     "bX *" -> bX, "bY *" -> bY, "bZ *" -> bZ, "\[Rho] *" -> \[Rho] , 
     "\[Beta] *" -> \[Beta] , "\[Delta] *" -> \[Delta] , "Eg" -> Eg}}
   }];
 (*{"Living-Epidermis",EFormulationLivingEpidermis,
 {"\[Nu]\[DoubleStruckG] *"\[Rule] \
\[Nu]\[DoubleStruckG],"c*"\[Rule]c,"k1 *"\[Rule] k1,"k2 *"\[Rule] \
k2,"q *"\[Rule] q,
 "bX *"\[Rule] bX,"bY *"\[Rule] bY,"bZ *"\[Rule] bZ,"\[Rho] *"\[Rule]\
\[Rho] ,"\[Beta] *"\[Rule]\[Beta] ,"\[Delta] *"\[Rule]\[Delta] ,"Eg"\
\[Rule]Eg}},
 {"Papillary-Dermis",EFormulationPapillaryDermis,
 {"\[Nu]\[DoubleStruckG] *"\[Rule] \
\[Nu]\[DoubleStruckG],"c*"\[Rule]c,"k1 *"\[Rule] k1,"k2 *"\[Rule] \
k2,"q *"\[Rule] q,
 "bX *"\[Rule] bX,"bY *"\[Rule] bY,"bZ *"\[Rule] bZ,"\[Rho] *"\[Rule]\
\[Rho] ,"\[Beta] *"\[Rule]\[Beta] ,"\[Delta] *"\[Rule]\[Delta] ,"Eg"\
\[Rule]Eg}},
 {"Reticular-Dermis",EFormulationReticularDermis,
 {"\[Nu]\[DoubleStruckG] *"\[Rule] \
\[Nu]\[DoubleStruckG],"c*"\[Rule]c,"k1 *"\[Rule] k1,"k2 *"\[Rule] \
k2,"q*"\[Rule] q,
 "bX *"\[Rule] bX,"bY *"\[Rule] bY,"bZ *"\[Rule] bZ,"\[Rho] *"\[Rule]\
\[Rho] ,"\[Beta] *"\[Rule]\[Beta] ,"\[Delta] *"\[Rule]\[Delta] ,"Eg"\
\[Rule]Eg}}*)

 SMTMesh["Stratum-Corneum", "H1", {Xnd, Ynd, ZndStratumCorneum}, 
  StratumCorneumMastermesh, "BodyID" -> "a", 
  "InterpolationOrder" -> 3];
 SMTMesh["Living-Epidermis", "H1", {Xnd, Ynd, ZndLivingEpidermis}, 
  LivingEpidermisMasterMesh, "BodyID" -> "a", 
  "InterpolationOrder" -> 3];
 SMTMesh["Papillary-Dermis", 
  "H1-rf", {Xnd/3, Ynd/3, ZndPapillaryDermis}, 
  PapillaryDermisMasterMesh, "BodyID" -> "a", 
  "InterpolationOrder" -> 3];
 SMTMesh["Reticular-Dermis", 
  "H1-rf", {Xnd/9, Ynd/9, ZndReticularDermis + 1}, 
  ReticularDermisMasterMesh, "BodyID" -> "a", 
  "InterpolationOrder" -> 3];


 Extract boundaries of the mesh;
 {MeshXmin, MeshXmax} = {Min[SMTNodes[[All, 2]]], 
   Max[SMTNodes[[All, 2]]]};
 {MeshYmin, MeshYmax} = {Min[SMTNodes[[All, 3]]], 
   Max[SMTNodes[[All, 3]]]};
 {MeshZmin, MeshZmax} = {Min[SMTNodes[[All, 4]]], 
   Max[SMTNodes[[All, 4]]]};

 Switch[AnalysisSwitch, 1,
  UNIAXIAL COMPRESSION X;
   Essential boundary conditions;
   SMTAddEssentialBoundary["X" == MeshXmin &, 1 -> 0];
  SMTAddEssentialBoundary["Y" == MeshYmin &, 2 -> 0];
  (*SMTAddEssentialBoundary["Y"\[Equal]MeshYmax&,2\[Rule] 0];*)

  SMTAddEssentialBoundary["Z" == MeshZmin &, 3 -> 0];

  Displacement condition X;
  Xdisplacement = -Prestress*(MeshXmax - MeshXmin);
  SMTAddEssentialBoundary["X" == MeshXmax &, 1 -> Xdisplacement];
  ,
  2,(*This is the step that I only want applied once in one session*)
  EQUIBIAXIAL COMPRESSION;
   Essential boundary conditions;
   SMTAddEssentialBoundary["X" == MeshXmin &, 1 -> 0];
  SMTAddEssentialBoundary["Y" == MeshYmin &, 2 -> 0];
  SMTAddEssentialBoundary["Z" == MeshZmin &, 3 -> 0];

  Displacement condition X;
  Xdisplacement = -CompressionStrainX*(MeshXmax - MeshXmin);
  SMTAddEssentialBoundary["X" == MeshXmax &, 1 -> Xdisplacement];

  Displacement condition Y;
  Ydisplacement = -CompressionStrainY*(MeshYmax - MeshYmin);
  SMTAddEssentialBoundary["Y" == MeshYmax &, 2 -> Ydisplacement];
  ,
  3,
  SHEAR ALONG Y;
   Essential boundary conditions;
   SMTAddEssentialBoundary["X" == MeshXmin &, 1 -> 0, 2 -> 0];
  SMTAddEssentialBoundary["Z" == MeshZmin &, 3 -> 0];

  Displacement condition Y;
  ShearDisplacement = ShearStrainY*(MeshYmax - MeshYmin);
  SMTAddEssentialBoundary["X" == MeshXmax &, 1 -> 0, 
   2 -> ShearDisplacement];

  ,
  4,
  SHEAR ALONG X;
   Essential boundary conditions;
   SMTAddEssentialBoundary["Y" == MeshYmin &, 1 -> 0, 2 -> 0];
  SMTAddEssentialBoundary["Z" == MeshZmin &, 3 -> 0];

  Displacement condition Y;
  ShearDisplacement = ShearStrainX*(MeshXmax - MeshXmin);
  SMTAddEssentialBoundary["Y" == MeshYmax &, 1 -> ShearDisplacement, 
   2 -> 0];
  ,
  5,
  SHEAR AND COMPRESSION;
   Essential boundary conditions;
   SMTAddEssentialBoundary["X" == MeshXmin &, 1 -> 0, 2 -> 0];
  SMTAddEssentialBoundary["Z" == MeshZmin &, 3 -> 0];

  Displacement condition Y;
  ShearDisplacement = -ShearStrainY*(MeshYmax - MeshYmin);
  SMTAddEssentialBoundary["X" == MeshXmax &, 2 -> ShearDisplacement];

  Displacement condition X;
  Xdisplacement = -CompressionStrainX*(MeshXmax - MeshXmin);
  SMTAddEssentialBoundary["X" == MeshXmax &, 1 -> Xdisplacement];

  ];

 Switch[SolverSwitch,
  In core PARDISO solver;
  1,
  SMTAnalysis["Output" -> AnalysisName <> ".out", 
   "OptimizeDll" -> True];
  SMTSetSolver[5, SolverMatrixType, {{60, 0}}];
  ,
  2,
  Out of core PARDISO solver;
  SMTAnalysis["Output" -> AnalysisName <> ".out", 
   "OptimizeDll" -> True];
  SMTSetSolver[5, SolverMatrixType, {{60, 2}}];
  ];
 EXTRACTION OF NODES AND ASSOCIATED DATA;
 Identify nodes of each skin layer;
 StratumCorneumNodes = 
  SMTFindNodes[{"DomainID", SMTFindDomains[{"Stratum-Corneum"}]}];
 LivingEpidermisNodes = 
  SMTFindNodes[{"DomainID", SMTFindDomains[{"Living-Epidermis"}]}];
 PapillaryDermisNodes = 
  SMTFindNodes[{"DomainID", SMTFindDomains[{"Papillary-Dermis"}]}];
 ReticularDermisNodes = 
  SMTFindNodes[{"DomainID", SMTFindDomains[{"Reticular-Dermis"}]}];

 Select nodes belonging to SC top lines;
 NodesSelectionYmax = 
  Table[4*i + (Xnd + 1)*4*(Xnd), {i, 1, Xnd + 1}];
 (*NodesSelectionYmid=Table[4*i+(Xnd+1)*4*(Xnd/2),{i,1,Xnd+1}]*) (* \
EVEN number of divisions *)

 NodesSelectionYmid = 
  Table[4*i + (Xnd + 1)*4*(Xnd + 1)/2, {i, 1, Xnd + 1}];
 NodesSelectionYmin = Table[4*i, {i, 1, Xnd + 1}];

 Initialise converged step counter;
 ic = 0;

 Output result file at zero step ;
 SMTDump[AnalysisName <> "_" <> ToString[ic]];


 Leading parameter;
 V0 = Velocity0; (* velocity *)

 Analysis characteristics;
 \[Gamma]0 = 0.01;
 \[CapitalDelta]\[Gamma]Min = 10^-6 V0;
 \[CapitalDelta]\[Gamma]Max = 0.1;
 \[Gamma]Max = 1;

 tolNR = 10^-8;
 tolNR0 = 10^-4; (* relaxed tolerance to go through critical points *)

 maxNR = 25; targetNR = 12;

 \[CapitalDelta]\[Lambda]Min = 10^-6 V0;
 \[CapitalDelta]\[Lambda]LineSearch = 
  100*\[CapitalDelta]\[Lambda]Min;
 stressgraph = {{0, 0}};
 SMTNextStep["t[\[Gamma]]" -> (#/V0 &), 
  "\[Lambda][\[Gamma]]" -> (# &), 
  "\[CapitalDelta]\[Gamma]" -> \[Gamma]0];
 plot\[Sigma]e = {{0, 0}};
 While[
  While[
   step = SMTConvergence[
     Conditional tolerance;
     Conditional tolerance;
     If[SMTRData["Multiplier"] > 0.9975*\[Gamma]Max, tolNR, tolNR0],
     maxNR, {"Adaptive \[Gamma]", 
      targetNR, \[CapitalDelta]\[Gamma]Min, \
\[CapitalDelta]\[Gamma]Max, \[Gamma]Max}]
   , 
   (*If[SMTIData["Iteration"]>4 || SMTRData[
   "MultiplierIncrement"]>\[CapitalDelta]\[Lambda]LineSearch,SMTRData[
   "LineSearchStepLength",1],SMTRData["LineSearchStepLength",0.1]];*)

      SMTNewtonIteration[];
   ];
  If[step[[4]] === "MinBound", 
   SMTStatusReport[
     "\[CapitalDelta]\[Gamma] < \[CapitalDelta]\[Gamma]Min"];];
  step[[3]] 
  , If[step[[1]]
   , SMTStepBack[]
   ];

  Display converged configuration;

  SMTShowMesh["DeformedMesh" -> True, Axes -> True, 
   AxesLabel -> {"X", "Y", "Z"}, "Field" -> "w", "Show" -> "Window"];
  (*AppendTo[plot\[Sigma]e,{SMTPostData["\[Lambda]11",{0,0,MeshZmax}],
  SMTPostData["\[Sigma]11(MPa)",Point[{0,0,MeshZmax}]]}];*)

  Save current converged increment configuration;
  ic = ic + 1;
  SMTDump[AnalysisName <> "_" <> ToString[ic]];

  SMTNextStep["t[\[Gamma]]" -> (#/V0 &), 
   "\[Lambda][\[Gamma]]" -> (# &), 
   "\[CapitalDelta]\[Gamma]" -> step[[2]]]
  ];


 SMTNextStep[100, 0];
 SMTNewtonIteration[];

 SMTDump[AnalysisName <> "_" <> ToString[ic] <> 
   ".After.Dissipation"];

 SMTSimulationReport[];

 ]

This is just the loop that I have at the moment, the actual code is much longer than this. I need step two within AnalysisSwitch to only run once and I'm not really sure how to get this to work. I was thinking of something along the lines of SMSIf[switch == 2 && (AnalysisSwitch = 1), SMSReturn[]]; but I have tried that and its not really working. Any ideas would be much appreciated. I apologize for not posting the whole code but it's quite long.

$\endgroup$
  • $\begingroup$ I have edited the title of the question to be more descriptive. I hope I got it right, since I am not sure if I understood the question? Next time, please, try to make the contents of the question shorter and more clear. Code example that you provided doesn't work (again) because some symbol definitions and AceGen libraries are missing. $\endgroup$ – Pinti Jan 19 at 21:03
  • $\begingroup$ Regarding the last paragraph, the code with SMSIf doesn't work because that function belongs to AceGen package and should only be used for finite element library generation. $\endgroup$ – Pinti Jan 19 at 21:08
6
+50
$\begingroup$

I am not sure if I have understood the question but let me give a minimal example of how to change the type of boundary conditions during AceFEM analysis. We will squeeze a box in one direction and then stretch it in the other direction.

First we define functions which will make the whole process more readable and organized.

Get["AceFEM`"]

setupPhase[noElements_Integer] := Module[
  {},
  SMTInputData[];
  (* Element libraries are downloaded from AceShare server. *)
  SMTAddDomain["Domain1", "OL:SED3H1DFHYH1NeoHooke", {}];
  SMTAddMesh[
   Raster3D@{
     {{{0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {1, 1, 0}}},
     {{{0, 0, 1}, {1, 0, 1}}, {{0, 1, 1}, {1, 1, 1}}}
     },
   "Domain1", "H1", noElements*{1, 1, 1}
   ];
  SMTAddEssentialBoundary[{{"X" == 0 &, 1 -> 0}, {"Y" == 0 &, 2 -> 0}, {"Z" == 0 &, 3 -> 0}}];
  SMTAnalysis[]
];

(* Handy function for visualization. *)
showMesh[imageSize_] := SMTShowMesh[
  "DeformedMesh" -> True,
  "BoundaryConditions" -> True,
  "Field" -> "u",
  "Contour" -> {False, 0, 0.5, 4},
  Axes -> True,
  AxesStyle -> 12,
  ImageSize -> imageSize,
  PlotRange -> {{-0.1, 2}, {-0.1, 1.1}, {-0.1, 1.1}}
];

(* Increase multiplier from 0 to 1 in specified number of steps.*)
analysisLoop[noSteps_Integer] := Module[
  {Δλ, results},
  Δλ = 1./noSteps;
  results = {};
  Do[
   SMTNextStep["Δλ" -> Δλ];
   While[SMTConvergence[10^-8, 15], SMTNewtonIteration[]];
   (* Collect any results here. *)
   AppendTo[results, showMesh[Medium]],
   {i, 1, noSteps}
   ];
  (* Results are returned from the function. *)
  results
];

Then we call the functions for setup phase and analysis loop. Between two loops we change the boundary conditions and reset multiplier (lambda) to start again from zero. Keep in mind that this is only one of many ways to achieve the same result.

setupPhase[4]

(* Set intensity of essential BC (displacement). *)
SMTAddEssentialBoundary[{"Y" == 1 &, 2 -> -0.25}];
(* Run analysis and collect results. *)
results = analysisLoop[5];
(* Reset intensity of the same essential BC (displacement is zero). *)
SMTAddEssentialBoundary[{"Y" == 1 &, 2 -> 0}, "Set" -> True];

(* Reset BC multiplier value. *)
SMTRData["Multiplier", 0.];
(* Set intensity of essential BC in other direction. *)
SMTAddEssentialBoundary[{"X" == 1 &, 1 -> 0.5}];
results = Join[results, analysisLoop[10]];

ListAnimate[results]

animation

|improve this answer|||||
$\endgroup$
  • $\begingroup$ thank you for your comprehensive answer. I have tried it and it did not work because there were a lot of details and steps in which node data needs to be gathered for successive steps. I was wondering if there was any way I could simply format my code to skip Switch[AnalysisSwitch ==1] at some section. What I mean is that just have one loop instead of two in the example you have presented. $\endgroup$ – IfIcantdoithomieitcantbedone Jan 21 at 12:28
  • $\begingroup$ @IfIcantdoithomieitcantbedone Sorry, but I don't understand your question. I suggest that you try to experiment with BC manipulation on a simple case like mine above and then apply that knowledge to more complex cases. I think that it is possible to have a single analysis loop and manipulate all BC values from within, but I don't see a good reason for that without knowing the full background of your example. But that is your task... $\endgroup$ – Pinti Jan 21 at 20:08

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.