Mixed 3D Solid Finite element

Generate the three-dimensional, eight node finite element for the analysis of hyperelastic solid problems. The element has the following characteristics:
=>	hexahedral topology,
=>	8 nodes,
=>	isoparametric mapping from the reference to the actual frame,

global unknowns are displacements of the nodes,
	u=(Subscript[u, i])(Subscript[N, i]), v=(Subscript[v, i])(Subscript[N, i]), w=(Subscript[w, i])(Subscript[N, i])
	enhanced strain formulation to improve shear and volumetric locking response,

	Δu=(({{(Subscript[u, ,X]), (Subscript[u, ,Y]), (Subscript[u, ,Z])}, {(Subscript[v, ,X]), (Subscript[v, ,Y]), (Subscript[v, ,Z])}, {(Subscript[w, ,X]), (Subscript[w, ,Y]), (Subscript[w, ,Z])}}))
	
	D=Δu+((Det[(Subscript[J, 0])])/(Det[J]))(({{ξ (Subscript[α, 1]), η (Subscript[α, 2]), ζ (Subscript[α, 3])}, {ξ (Subscript[α, 4]), η (Subscript[α, 5]), ζ (Subscript[α, 6])}, {ξ (Subscript[α, 7]), η (Subscript[α, 8]), ζ (Subscript[α, 9])}})).((Subscript[J, 0])^(-1))
	
	where α={(Subscript[α, 1]),(Subscript[α, 2]),...,(Subscript[α, 9])} are internal degrees of freedom eliminated at the element level.
	the classical hyperelastic Neo-Hooke's potential energy,

	Π=(∫_(Subscript[Ω, 0])^( ))((λ/2)(((det F-1))^2)+μ(((Tr[C]-2)/2)-Log[det F])-u .Q )d(Subscript[Ω, 0]),
	
	where C=(F^T)F is right Cauchy-Green tensor, F=I+D  is enhanced deformation gradient.

AceFEM

<<"AceGen`";
SMSInitialize["ExamplesHypersolid3D","Environment"->"AceFEM"];
SMSTemplate["SMSTopology"->"H1","SMSNoDOFCondense"->9]

SMSStandardModule["Tangent and residual"];
SMSDo[IpIndex,1,SMSInteger[es$$["id","NoIntPoints"]]];

{Xi,Yi,Zi}⊨Array[SMSReal[nd$$[#2,"X",#1]]&,{3,8}];
  {ut,vt,wt}⊨Array[SMSReal[nd$$[#2,"at",#1]]&,{3,8}];
  SMSGroupDataNames={"Elastic modulus","Poisson ratio"};
  SMSDefaultData={21000,0.3};
  {Em,ν}⊨SMSReal[Array[es$$["Data",#1]&,2]];
  {ξ,η,ζ,wGauss}⊢Array[SMSReal[es$$["IntPoints",#1,IpIndex]]&,4];
  {ξi,ηi,ζi}={{-1,1,1,-1,-1,1,1,-1},{-1,-1,1,1,-1,-1,1,1},{-1,-1,-1,-1,1,1,1,1}};
  Ni⊨MapThread[ 1/8(1+ξ #1)(1+η #2)(1+ζ #3)&,{ξi,ηi,ζi}];
  {X,Y,Z}⊨SMSFreeze[{Ni.Xi,Ni.Yi,Ni.Zi}];
  Jm⊨SMSD[{X,Y,Z},{ξ,η,ζ}];Jd⊨Det[Jm];
  {u,v,w}⊨{Ni.ut,Ni.vt,Ni.wt};
  Δu⊨SMSD[{u,v,w},{X,Y,Z},"Implicit"->{{{ξ,η,ζ},{X,Y,Z},SMSInverse[Jm]}}];
  Jm0⊨SMSReplaceAll[Jm,{ξ->0,η->0,ζ->0}];
  α⊨SMSReal[Array[ed$$["ht",#]&,9]];
  H0={{ξ α[[1]],η α[[2]],ζ α[[3]]},{ξ α[[4]],η α[[5]],ζ α[[6]]},{ξ α[[7]],η α[[8]],ζ α[[9]]}};
  H=((Det[Jm0] )/Jd)H0.SMSInverse[Jm0];
  𝔻⊨Δu+H;
  𝔽⊨SMSFreeze[IdentityMatrix[3]+𝔻];
  ℂ⊨Transpose[𝔽].𝔽;
  𝕁⊨Det[𝔽];
  {λ,μ}⊨SMSHookeToLame[Em,ν];
  Π⊨(1/2) λ (((𝕁-1))^2)+μ ((1/2) (Tr[ℂ]-2)-Log[𝕁]);
  𝕒=Flatten[{Transpose[{ut,vt,wt}],α}];

SMSDo[i,1,Length[𝕒]];
    Ψi⊨Jd wGauss  SMSD[Π,𝕒,i];
    SMSExport[SMSResidualSign Ψi,p$$[i],"AddIn"->True];
    SMSDo[j,i,Length[𝕒]];
      Kij⊨SMSD[Ψi,𝕒,j];
      SMSExport[Kij,s$$[i,j],"AddIn"->True];
    SMSEndDo[]; 
  SMSEndDo[]; 
SMSEndDo[];

SMSStandardModule["Postprocessing"];
SMSDo[IpIndex,1,SMSInteger[es$$["id","NoIntPoints"]]];
  SMSEvaluateCellsWithTag["ExamplesHypersolid3D","Session"];
  SMSNPostNames={"DeformedMeshX","DeformedMeshY","DeformedMeshZ","u","v","w"};
  SMSExport[{ut,vt,wt,ut,vt,wt}//Transpose,npost$$];
  SMSGPostNames={"Sxx","Sxy","Sxz","Syx","Syy","Syz","Szx","Szy","Szz",
"Exx","Exy","Exz","Eyx","Eyy","Eyz","Ezx","Ezy","Ezz"};
  Cauchy=(1/𝕁)SMSD[Π,𝔽] .Transpose[𝔽];
  GreenLagrange=1/2(ℂ-IdentityMatrix[3]);
  SMSExport[{Cauchy,GreenLagrange}//Flatten,gpost$$[IpIndex,#]&];
SMSEndDo[];

SMSWrite[];

FEAP

Regenerate the three-dimensional, eight node finite element from chapter Mixed 3D Solid FE for AceFEM for FEAP environment.

<<"AceGen`";
SMSInitialize["ExamplesHypersolidFEAP","Environment"->"FEAP"];
SMSTemplate["SMSTopology"->"H1","SMSNoDOFCondense"->9]

SMSStandardModule["Tangent and residual"];
SMSDo[IpIndex,1,SMSInteger[es$$["id","NoIntPoints"]]];

{Xi,Yi,Zi}⊨Array[SMSReal[nd$$[#2,"X",#1]]&,{3,8}];
  {ut,vt,wt}⊨Array[SMSReal[nd$$[#2,"at",#1]]&,{3,8}];
  SMSGroupDataNames={"Elastic modulus","Poisson ratio"};
  {Em,ν}⊨SMSReal[Array[es$$["Data",#1]&,2]];
  {ξ,η,ζ,wGauss}⊢Array[SMSReal[es$$["IntPoints",#1,IpIndex]]&,4];
  {ξi,ηi,ζi}={{-1,1,1,-1,-1,1,1,-1},{-1,-1,1,1,-1,-1,1,1},{-1,-1,-1,-1,1,1,1,1}};
  Ni⊨MapThread[ 1/8(1+ξ #1)(1+η #2)(1+ζ #3)&,{ξi,ηi,ζi}];
  {X,Y,Z}⊨SMSFreeze[{Ni.Xi,Ni.Yi,Ni.Zi}];
  Jm⊨SMSD[{X,Y,Z},{ξ,η,ζ}];Jd⊨Det[Jm];
  {u,v,w}⊨{Ni.ut,Ni.vt,Ni.wt};
  Δu⊨SMSD[{u,v,w},{X,Y,Z},"Implicit"->{{{ξ,η,ζ},{X,Y,Z},SMSInverse[Jm]}}];
  Jm0⊨SMSReplaceAll[Jm,{ξ->0,η->0,ζ->0}];
  α⊨SMSReal[Array[ed$$["ht",#]&,9]];
  H0={{ξ α[[1]],η α[[2]],ζ α[[3]]},{ξ α[[4]],η α[[5]],ζ α[[6]]},{ξ α[[7]],η α[[8]],ζ α[[9]]}};
  H=((Det[Jm0] )/Jd)H0.SMSInverse[Jm0];
  𝔻⊨Δu+H;
  𝔽⊨IdentityMatrix[3]+𝔻;
  ℂ⊨Transpose[𝔽].𝔽;𝕁⊨Det[𝔽];
  {λ,μ}⊨SMSHookeToLame[Em,ν];
  Π⊨(1/2) λ (((𝕁-1))^2)+μ ((1/2) (Tr[ℂ]-2)-Log[𝕁]);
  𝕒=Flatten[{Transpose[{ut,vt,wt}],α}];

SMSDo[i,1,SMSNoAllDOF];
    Ψi⊨Jd wGauss  SMSD[Π,𝕒,i];
    SMSExport[SMSResidualSign Ψi,p$$[i],"AddIn"->True];
    SMSDo[j,i,SMSNoAllDOF];
      Kij⊨SMSD[Ψi,𝕒,j];
      SMSExport[Kij,s$$[i,j],"AddIn"->True];
    SMSEndDo[]; 
  SMSEndDo[]; 
SMSEndDo[];

SMSWrite[];

Here is the FEAP input data file for the test example from the chapter  Mixed 3D Solid FE for AceFEM . You need to install FEAP environment in order to run the example.

feap
0,0,0,3,3,8

block
cart,6,15,6,1,1,1,10
1,10.,0.,0.
2,10.,2.,0.
3,0.,2.,0.
4,0.,0.,0.
5,10.,0.,2.
6,10.,2.,2.
7,0.,2.,3.
8,0.,0.,3.

ebou
1,0,1,1,1
1,10.,,,1

edisp,add
1,10.,,,-1.

mate,1
user,10
1000,0.3

end

macr
tol,,1e-9
prop,,1
dt,,1
loop,,5
time
loop,,10
tang,,1
next
disp,,340
next
end

stop

Here is the generated element compiled and linked into the FEAP's Visual Studio project. See Install.txt for details. The SMSFEAPRun function then starts FEAP with a beam.inp file as a standard FEAP input file and a beam.out file as output file.

SMSFEAPMake["ExamplesHypersolidFEAP"]

SMSFEAPRun["feap.inp"]

ReadList["feap.out","Record"][[-4]]

ELFEN

Regenerate the three-dimensional, eight node finite element from chapter Mixed 3D Solid FE for AceFEM for ELFEN environment.

The AceGen input presented in previous example can be used again with the "Environment"->"ELFEN" option to produce Elfen's source code file. However, due to the non-standard approach to the implementation of the Newton-Raphson loop in ELFEN  result would not be the most efficient. More efficient implementation is obtained if the evaluation of the tangent matrix and residual vector are separated. The procedure is controlled by the values of  environment constants "SkipTangent", "SkipResidual" and "SubIterationMode". 
When the tangent matrix is required the variables are set to
idata$$["SkipTangent"]=0,
idata$$["SkipResidual"]=1,
idata$$["SubIterationMode"]=1
and when the residual is required the variables are set to
idata$$["SkipTangent"]=1,
idata$$["SkipResidual"]=0,
idata$$["SubIterationMode"]=0.
Additionally, the non-standard evaluation of the Newton-Raphson loop makes implementation of the mixed FE models difficult. Thus only displacement element is generated.
The generated code is then incorporated into ELFEN.

<<"AceGen`";
SMSInitialize["ExamplesHypersolidELFEN","Environment"->"ELFEN"];
SMSTemplate["SMSTopology"->"H1"]

SMSStandardModule["Tangent and residual"];
SMSDo[IpIndex,1,SMSInteger[es$$["id","NoIntPoints"]]];
  {Xi,Yi,Zi}⊨Array[SMSReal[nd$$[#2,"X",#1]]&,{3,8}];
  {ut,vt,wt}⊨Array[SMSReal[nd$$[#2,"at",#1]]&,{3,8}];
  SMSGroupDataNames={"Elastic modulus","Poisson ratio"};
  {Em,ν}⊨SMSReal[Array[es$$["Data",#1]&,2]];
  {ξ,η,ζ,wGauss}⊢Array[SMSReal[es$$["IntPoints",#1,IpIndex]]&,4];
  {ξi,ηi,ζi}={{-1,1,1,-1,-1,1,1,-1},{-1,-1,1,1,-1,-1,1,1},{-1,-1,-1,-1,1,1,1,1}};
  Ni⊨MapThread[ 1/8(1+ξ #1)(1+η #2)(1+ζ #3)&,{ξi,ηi,ζi}];
  {X,Y,Z}⊨SMSFreeze[{Ni.Xi,Ni.Yi,Ni.Zi}];
  Jm⊨SMSD[{X,Y,Z},{ξ,η,ζ}];Jd⊨Det[Jm];
  {u,v,w}⊨{Ni.ut,Ni.vt,Ni.wt};
  Δu⊨SMSD[{u,v,w},{X,Y,Z},"Implicit"->{{{ξ,η,ζ},{X,Y,Z},SMSInverse[Jm]}}];
  Jm0⊨SMSReplaceAll[Jm,{ξ->0,η->0,ζ->0}];
  𝔽⊨IdentityMatrix[3]+Δu;
  ℂ⊨Transpose[𝔽].𝔽;𝕁⊨Det[𝔽];
  {λ,μ}⊨SMSHookeToLame[Em,ν];
  Π⊨(1/2) λ (((𝕁-1))^2)+μ ((1/2) (Tr[ℂ]-2)-Log[𝕁]);
  𝕒=Flatten[Transpose[{ut,vt,wt}]];
  SMSIf[SMSInteger[idata$$["SkipTangent"]]==1];
    SMSDo[i,1,SMSNoAllDOF];
      Ψi⊨Jd wGauss  SMSD[Π,𝕒,i];
      SMSExport[SMSResidualSign Ψi,p$$[i],"AddIn"->True];
    SMSEndDo[]; 
  SMSElse[];
    SMSDo[i,1,SMSNoAllDOF];
      Ψi⊨Jd wGauss  SMSD[Π,𝕒,i];
      SMSDo[j,i,SMSNoAllDOF];
        Kij⊨SMSD[Ψi,𝕒,j];
        SMSExport[Kij,s$$[i,j],"AddIn"->True];
      SMSEndDo[]; 
    SMSEndDo[]; 
  SMSEndIf[];
SMSEndDo[];

SMSWrite[];

Here is the generated element compiled and linked into the ELFEN's Visual Studio project. See Install.txt for details. The SMSELFENRun function then starts ELFEN with a ELFENExample.dat file as a input file and a tmp.res file as output file. The ELFEN input data file for the one element test example is available in a $BaseDirectory/Applications/AceGen/Include/ELFEN/ directory.

SMSELFENMake["ExamplesHypersolidELFEN"]

SMSELFENRun["ELFEN.dat"]

ABAQUS

Regenerate the three-dimensional, eight node finite element from chapter Mixed 3D Solid FE for AceFEM for ABAQUS environment.

The AceGen input presented in previous example can be used again with the "Environment"->"ABAQUS" option to produce ABAQUS's source code file. The current ABAQUS interface does not support internal degrees of freedom. Consequently, the mixed deformation modes are skipped. The generated code is then incorporated into ABAQUS.

<<"AceGen`";
SMSInitialize["ExamplesHypersolidABAQUS","Environment"->"ABAQUS"];
SMSTemplate["SMSTopology"->"H1"]

SMSStandardModule["Tangent and residual"];
SMSDo[IpIndex,1,SMSInteger[es$$["id","NoIntPoints"]]];
  {Xi,Yi,Zi}⊨Array[SMSReal[nd$$[#2,"X",#1]]&,{3,8}];
  {ut,vt,wt}⊨Array[SMSReal[nd$$[#2,"at",#1]]&,{3,8}];
  SMSGroupDataNames={"Elastic modulus","Poisson ratio"};
  {Em,ν}⊨SMSReal[Array[es$$["Data",#1]&,2]];
  {ξ,η,ζ,wGauss}⊢Array[SMSReal[es$$["IntPoints",#1,IpIndex]]&,4];
  {ξi,ηi,ζi}={{-1,1,1,-1,-1,1,1,-1},{-1,-1,1,1,-1,-1,1,1},{-1,-1,-1,-1,1,1,1,1}};
  Ni⊨MapThread[ 1/8(1+ξ #1)(1+η #2)(1+ζ #3)&,{ξi,ηi,ζi}];
  {X,Y,Z}⊨SMSFreeze[{Ni.Xi,Ni.Yi,Ni.Zi}];
  Jm⊨SMSD[{X,Y,Z},{ξ,η,ζ}];Jd⊨Det[Jm];
  {u,v,w}⊨{Ni.ut,Ni.vt,Ni.wt};
  Δu⊨SMSD[{u,v,w},{X,Y,Z},"Implicit"->{{{ξ,η,ζ},{X,Y,Z},SMSInverse[Jm]}}];
  Jm0⊨SMSReplaceAll[Jm,{ξ->0,η->0,ζ->0}];
  𝔽⊨IdentityMatrix[3]+Δu;
  ℂ⊨Transpose[𝔽].𝔽;𝕁⊨Det[𝔽];
  {λ,μ}⊨SMSHookeToLame[Em,ν];
  Π⊨(1/2) λ (((𝕁-1))^2)+μ ((1/2) (Tr[ℂ]-2)-Log[𝕁]);
  𝕒=Flatten[Transpose[{ut,vt,wt}]];
  SMSDo[i,1,SMSNoAllDOF];
    Ψi⊨Jd wGauss  SMSD[Π,𝕒,i];
    SMSDo[j,i,SMSNoAllDOF];
      Kij⊨SMSD[Ψi,𝕒,j];
      SMSExport[Kij,s$$[i,j],"AddIn"->True];
    SMSEndDo[]; 
  SMSEndDo[]; 
SMSEndDo[];

SMSWrite[];

Here is the generated element compiled and linked into the ABAQUS's. The SMSABAQUSRun function then starts ABAQUS with a ABAQUSExample.inp file as a input file and a tmp.res file as output file. The ABAQUS input data file for the one element test example is available in a $BaseDirectory/Applications/AceGen/Include/ABAQUS/ directory.

SMSABAQUSRun["ABAQUSExample","UserElement"->"ExamplesHypersolidABAQUS"]

AceFEM Test example

You need to install AceFEM package in order to run the example.

<<AceFEM`;
SMTInputData[];
SMTAddDomain["A","ExamplesHypersolid3D",{"E*"->1000.,"Poi*"->.3}];
SMTAddEssentialBoundary[{ "X"==0&,0,0,0},{ "X"==10 && "Z"==0 &,,,-1}];
SMTMesh["A","H1",{20,8,8},
{
{{{0,0,0},{5,0,-1},{10,0,0}},{{0,2,0},{5,2,-1},{10,2,0}}},{{{0,0,3},{5,0,1},{10,0,4}},{{0,2,3},{5,2,1},{10,2,4}}}
}
];
SMTAnalysis[];

graph={};
SMTNextStep[1,1];
While[
While[step=SMTConvergence[10^-8,10,{"Adaptive",8,.001,1,10}], SMTNewtonIteration[];];
SMTStatusReport[];
If[Not[step[[1]]],
SMTShowMesh["DeformedMesh"->True,"Show"->"Window"|{"Animation","Bending"},"Field"->"Sxx",
"User"->SMTShowMesh["Elements"->False,"RawOutput"->True][[5]]];
AppendTo[graph,{SMTPostData["w",{10,1,2}],SMTPostData["Sxx",{0,1,2}]}]
];
If[step[[4]]==="MinBound",Print["Error: Δλ < Δλmin"]];
step[[3]] 
,If[step[[1]],SMTStepBack[];];
SMTNextStep[1,step[[2]]]
];

SMTPostData[{"u","v","w","Sxx","Sxz"},{10,1,1}]

ListLinePlot[graph,Filling->Axis,AxesLabel->{"w","(Subscript[σ, xx])"}]

Show[SMTShowMesh["Elements"->False],
SMTShowMesh["DeformedMesh"->True,"Mesh"->False,"Field"->Map[Norm,SMTNodeData["at"]]]]

SMTMakeAnimation["Bending"]