0
$\begingroup$

hope everybody is well and safe. So, I'm creating an AceGen element for the inflation of an anisotropic hyperelastic tube (2-fibre families). I have created the element in a strip formulation and it was successful. With the tube, the solution is a little more complex as I have to keep the fibre directions tangent to the surface of the tube in any position. I know what I need to do, but I'm not sure exactly how to do it.

Here is a snippet of my code that ends with the definition of the strain energy density:

<< AceGen`; << AceFEM`;
GLmphEnvironment = "AceFEM";
GLmphElementName = "Strip_Test_GOH";
GLmphSymmetricTangent = True;
Default material properties;
E\[DoubleStruckG] = 50000;
\[Nu]\[DoubleStruckG] = 0.45;
\[Rho]0 = 7500;
c1 = 7.64;
\[Kappa]g = 100000;
k1 = 996.6; 
k2 = 524.6;
\[Kappa] = 0.226;
q = 0.837758;
Initialization Procedure;
SMSInitialize[GLmphElementName, "Environment" -> "AceFEM", 
  "Mode" -> "Optimal"];
SetOptions[SMSPrint, "Output" -> "File", 
  "Condition" -> "DebugElement"];
SMSTemplate["SMSTopology" -> "H2", "SMSSymmetricTangent" -> True, 
  "SMSDomainDataNames" -> {
    "E -elastic modulus",
    "\[Nu] -poisson ratio",
    "\[Rho] -density", 
    "c1 -first MR parameter",
    "\[Kappa]g -bulk modulus",
    "k1 -first GOH parameter",
     "k2 -secondv GOH parameter",
     "\[Kappa] -dispersion parameter",
    "q -mean orientation angle of collagen fibres"} , 
  "SMSDefaultData" -> {E\[DoubleStruckG], \[Nu]\[DoubleStruckG], \
\[Rho]0, c1, \[Kappa]g, k1, k2, \[Kappa], q }
  ];
Constitutive equations;
ElementDefinitions[task_] := (
  \[CapitalXi] = {\[Xi], \[Eta], \[Zeta]} \[RightTee] 
    Table[SMSReal[es$$["IntPoints", i, Ig]], {i, 3}];
  wgp \[RightTee] SMSReal[es$$["IntPoints", 4, Ig]];
  \[DoubleStruckCapitalX]IO \[RightTee] 
   Table[SMSReal[nd$$[i, "X", j]], {i, SMSNoNodes}, {j, 
 SMSNoDimensions}];
  n = 2; zero = Table[{s, 0}, {s, -1, 1, 2/n}];
  \[DoubleStruckCapitalN]\[Xi] \[DoubleRightTee] 
   Simplify[
Table[InterpolatingPolynomial[
  ReplacePart[zero, {i, 2} -> 1], \[Xi]], {i, n + 1}]];
  \[DoubleStruckCapitalN]\[Eta] \[DoubleRightTee] 
   Simplify[
Table[InterpolatingPolynomial[
  ReplacePart[zero, {i, 2} -> 1], \[Eta]], {i, n + 1}]];
  \[DoubleStruckCapitalN]\[Zeta] \[DoubleRightTee] 
   Simplify[
Table[InterpolatingPolynomial[
  ReplacePart[zero, {i, 2} -> 1], \[Zeta]], {i, n + 1}]];
  \[DoubleStruckCapitalN]2 \[DoubleRightTee] 
   Table[\[DoubleStruckCapitalN]\[Xi][[
  i]] \[DoubleStruckCapitalN]\[Eta][[
  k]] \[DoubleStruckCapitalN]\[Zeta][[l]], {i, 1, n + 1}, {k, 1, 
 n + 1}, {l, 1, n + 1}];
  \[DoubleStruckCapitalN]h \[DoubleRightTee] 
   Extract[\[DoubleStruckCapitalN]2, {{1, 1, 1}, {3, 1, 1}, {1 + 2, 
  1 + 2, 1}, {1, 1 + 2, 1}, {1, 1, 1 + 2}, {1 + 2, 1, 
  1 + 2}, {1 + 2, 1 + 2, 1 + 2}, {1, 1 + 2, 1 + 2}, {1 + 1, 1, 
  1}, {1 + 2, 1 + 1, 1}, {1 + 1, 1 + 2, 1}, {1, 1 + 1, 1}, {1, 1, 
  1 + 1}, {1 + 2, 1, 1 + 1}, {1 + 2, 1 + 2, 1 + 1}, {1, 1 + 2, 
  1 + 1}, {1 + 1, 1, 1 + 2}, {1 + 2, 1 + 1, 1 + 2}, {1 + 1, 1 + 2,
   1 + 2}, {1, 1 + 1, 1 + 2}, {1 + 1, 1 + 1, 1}, {1 + 1, 1, 
  1 + 1}, {1 + 2, 1 + 1, 1 + 1}, {1 + 1, 1 + 2, 1 + 1}, {1, 1 + 1,
   1 + 1}, {1 + 1, 1 + 1, 1 + 2}, {1 + 1, 1 + 1, 1 + 1}}];
  Position;
  \[DoubleStruckCapitalX] \[RightTee] 
   SMSFreeze[\[DoubleStruckCapitalN]h.\[DoubleStruckCapitalX]IO];
  Jm \[DoubleRightTee] SMSD[\[DoubleStruckCapitalX], \[CapitalXi]];
  Jmi \[DoubleRightTee] SMSInverse[Jm];
  Jd \[DoubleRightTee] SMSDet[Jm];
  Displacement at current increment;
  \[DoubleStruckU]IO \[RightTee] 
   SMSReal[Table[
 nd$$[i, "at", j], {i, SMSNoNodes}, {j, SMSDOFGlobal[[i]]}]];
  Flattened vectors of degrees of freedom;
  (**\[DoubleStruckP]e\[DoubleRightTee]Flatten[\[DoubleStruckU]IO];**)
\

  \[DoubleStruckP]e \[DoubleRightTee] Flatten[\[DoubleStruckU]IO];
  \[DoubleStruckU] \[DoubleRightTee] \[DoubleStruckCapitalN]h.\
\[DoubleStruckU]IO;(*Displacement field *)
  Current position;
  (**\[DoubleStruckX]\[DoubleRightTee]\[DoubleStruckCapitalX]+\
\[DoubleStruckU];**)
  Displacement gradient;
  \[DoubleStruckCapitalD] \[DoubleRightTee] 
   SMSD[\[DoubleStruckU], \[DoubleStruckCapitalX], 
    "Dependency" -> {\[CapitalXi], \[DoubleStruckCapitalX], Jmi}];
  \[DoubleStruckCapitalF] \[DoubleRightTee] 
   IdentityMatrix[3] + \[DoubleStruckCapitalD];
  fGauss \[DoubleRightTee] Jd;
  Kinematics;
  SMSFreeze[\[DoubleStruckCapitalF]f, \[DoubleStruckCapitalF], 
   "Ignore" -> NumberQ];
  \[DoubleStruckCapitalJ] \[DoubleRightTee] 
   SMSDet[\[DoubleStruckCapitalF]f];
  \[DoubleStruckCapitalC] \[DoubleRightTee] 
   Transpose[\[DoubleStruckCapitalF]f].\[DoubleStruckCapitalF]f;
  SMSFreeze[\[DoubleStruckCapitalC]iso, 
   SMSPower[\[DoubleStruckCapitalJ] , -2/3] \[DoubleStruckCapitalC], 
   "Symmetric" -> True, "Ignore" -> NumberQ];
  I1iso \[DoubleRightTee] Tr[\[DoubleStruckCapitalC]iso];
  (*I2iso\[DoubleRightTee](1/2) (I1iso^2-
  Tr[\[DoubleStruckCapitalC]iso.\[DoubleStruckCapitalC]iso]);*)
  q = 0.837758;
  a4 \[DoubleRightTee] {Cos[q], Sin[q], 0}; (*These are the fibre directions*)
  a6 \[DoubleRightTee] {Cos[q], -Sin[q], 0};
  I4iso \[DoubleRightTee] 
   SMSPower[\[DoubleStruckCapitalJ], (-2/
       3)] (a4.(\[DoubleStruckCapitalC].a4));
  I6iso \[DoubleRightTee] 
   SMSPower[\[DoubleStruckCapitalJ], (-2/
       3)] (a6.(\[DoubleStruckCapitalC].a6));
  (*positive powers greater than 1*)
  Assignment of material properties;
  {E\[DoubleStruckG], \[Nu]\[DoubleStruckG], \[Rho]0, c1,  \[Kappa]g, 
    k1, k2, \[Kappa], q} \[RightTee] 
   SMSReal[Table[es$$["Data", i], {i, Length[SMSDomainDataNames]}]];
  Strain energy density;
  Volumetric Component;
  Wvol \[DoubleRightTee] (\[Kappa]g/2) (\[DoubleStruckCapitalJ] - 1)^2 
  (*U=(1/4) K\[DoubleStruckG] ((JFe-1)^2+Log[JFe]^2);*);
  Decoupled Mooney - Rivlin;
  Wiso \[DoubleRightTee] c1 (I1iso - 3)(*+c2 (I2iso-3);*);
  GOH Model;
  Wfibre4 \[DoubleRightTee] (k1/(2*k2))*((Exp[
        k2 (\[Kappa] I1iso + ((1 - 3 \[Kappa]) I4iso) - 1)^2]) - 1);
  Wfibre6 \[DoubleRightTee] (k1/(2*k2))*((Exp[
        k2 (\[Kappa] I1iso + ((1 - 3 \[Kappa]) I6iso) - 1)^2]) - 1);
 

With the fibre directions being a4 and a6. I've had a look at Automation of the Finite Element Method by Korelc and Wriggers, and I know I should do something that kind of looks like this, which is essentially using convective coordinates:

    {ui, vi, wi} \[DoubleRightTee] 
  Transpose[
   Table[{\[DoubleStruckU]eIO[[i, 1]], \[DoubleStruckU]eIO[[i, 
      2]], \[DoubleStruckU]eIO[[i, 3]]}, {i, SMSNoNodes}]];
u \[DoubleRightTee] SMSFreeze[\[DoubleStruckCapitalN]h.ui];
v \[DoubleRightTee] SMSFreeze[\[DoubleStruckCapitalN]h.vi];
w \[DoubleRightTee] SMSFreeze[\[DoubleStruckCapitalN]h.wi];
X \[RightTee] SMSFreeze[\[DoubleStruckCapitalN]h.Xi];
Y \[RightTee] SMSFreeze[\[DoubleStruckCapitalN]h.Yi];
Z \[RightTee] SMSFreeze[\[DoubleStruckCapitalN]h.Zi];

x \[DoubleRightTee] {u + X, v + Y, w + Z};

r\[Xi]x \[DoubleRightTee] SMSD[x, \[Xi]];
r\[Eta]x \[DoubleRightTee] SMSD[x, \[Eta]];


r\[Xi]nx \[DoubleRightTee] r\[Xi]x/SMSSqrt[r\[Xi]x];
r\[Eta]nx \[DoubleRightTee] r\[Eta]x/SMSSqrt[r\[Eta]x];

I have two main questions here:

  1. Is the basis of my code even correct? I have a feeling that I should be using the SMSCovariantBase to transform the x,y,z reference vectors to cylindrical coordinates. The cylindrical mesh is only outlined in the AceFem file so I'm wondering if I should also use cylindrical coordinates for the acegen element.
  2. If I could keep the fibre directions locally tangent to element surface, what is the mathematical relation between the tangent vector and the vector of fibre directions that I have outlined? If this is possible, I must be missing something so simple, like a multiplication operation or something.

Thanks for your help and time.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.