Finite element engineering applications based on FEMPlayer and Mathematica Player Pro

Joze Korelc and Teja Melink

University of Ljubljana
Jamova 2, SI-1000, Ljubljana, Slovenia, EU

Matija Gams
Slovenian National Building and Civil Engineering Institute
Dimičeva 12, SI-1000 Ljubljana, Slovenia, EU

Contents

Customized finite element based engineering applications?

AceFEM ⟺ FEMPlayer  and why we need yet another version?

Building simple FEMPlayer application

Application: linear buckling analysis

Application: Verification of the seismic resistance of the masonry building structure

Customized finite element based engineering applications

AceFEM ⟺ FEMPlayer

AceFEM ⟺ FEMPlayer

Typical FEMPlayer application bundle consists of:
1	FEMPlayer installation (32 and 64 version available for all platforms Windows/Mac/Linux)
2	dll files of all elements (32 and 64 version available for all platforms Windows/Mac/Linux)
3	main application notebook (incorporates all application specific definitions)

Model Application: Pnevmatic Tyre

Simulation require:
=>	large strain hyperelastic Neo-Hooke type solid finite element capable of describing rubber material
=>	gas pressure finite element with the surface traction perpendicular to the deformed surface of the domain
=>	procedure to find initial undeformed shape of the tyre

The first task is to find out the shape of the empty tyre.  We only know that the shape of the weightless tyre is circular. The shape of the empty tyre can be obtained by starting with the weightless tyre and then increasing the volume force to its final value. The gas pressure is then increased in order to obtain the shape of the tyre for various gas pressure levels.

Solid finite element - hyperelastic rubber material

Generate two-dimensional, four node finite element for the analysis of the steady state problems in the mechanics of solids. The element has the following characteristics:

=>	quadrilateral topology,
=>	4 node element,
=>	isoparametric mapping from the reference to the actual frame,
=>	global unknowns are displacements of the nodes,
=>	the element should allow arbitrary large displacements and rotations,
=>	the problem is defined by the hyperelastic Neo-Hooke type strain energy potential,
	Π=(∫_(Subscript[Ω, 0])^( ))((λ/2)(((det F-1))^2)+μ(((Tr[C]-3)/2)-Log[det F])-u .Q )d(Subscript[Ω, 0]),
	where C=(F^T)F is right Cauchy-Green tensor, F=I+▽u  is deformation gradient, 
	u is displacements field, Q is prescribed body force, t is prescribed surface traction, 
	(Subscript[Ω, 0]) is the initial domain of the problem and λ, μ are the first and the second Lame's material 		constants.

<<AceGen`;
SMSInitialize["Rubber","Environment"->"AceFEM"];
SMSTemplate["SMSTopology"->"Q1",
 "SMSGroupDataNames"-> {"E - Elastic modulus", "ν - Poisson ratio",
      "t - Thickness of the element","Qx -Volume load X","Qy -Volume load Y"}
,"SMSDefaultData"->{21000,0.3,1,0,0}
]
SMSStandardModule["Tangent and residual"];
    Xi⊨Array[SMSReal[nd$$[#,"X",1],{{-1,1,1,-1}[[#]]+SMSRandom[]}]&,4];
    Yi⊨Array[SMSReal[nd$$[#,"X",2],{{-1,-1,1,1}[[#]]+SMSRandom[]}]&,4];
    ui⊨Array[SMSReal[nd$$[#,"at",1]]&,4];
    vi⊨Array[SMSReal[nd$$[#,"at",2]]&,4];
    {Em, ν, thick,Qx,Qy}⊨Array[SMSReal[es$$["Data",#]] &, 5];
    SMSDo[IpIndex,1,SMSInteger[es$$["id","NoIntPoints"]]];
    {ξ,η,ζ,wGauss}⊢Array[SMSReal[es$$["IntPoints",#1,IpIndex]]&,4];
	Ni⊨1/4 {(1-ξ)(1-η),(1+ξ)(1-η),(1+ξ)(1+η),(1-ξ)(1+η)};
    X⊨SMSFreeze[Ni.Xi];
    Y⊨SMSFreeze[Ni.Yi];
    Jm⊨SMSD[{X,Y},{ξ,η}];Jd⊨Det[Jm];
    SMSDefineDerivative[{ξ,η},{X,Y},SMSSimplify[SMSInverse[Jm]]];
    u⊨Ni.ui;v⊨Ni.vi;
    𝔻⊨SMSD[{u,v},{X,Y}];
    𝔽⊨IdentityMatrix[2]+𝔻;
    ℂ⊨Transpose[𝔽].𝔽;
  𝕁⊨Det[𝔽];
  {λ,μ}⊨SMSHookeToLame[Em,ν];
  mult⊨SMSReal[rdata$$["Multiplier"]];
    Π⊨λ/2(𝕁-1)^2+μ((Tr[ℂ]-2)/2-Log[𝕁])-mult {Qx,Qy}.{u,v}; 
    𝕒={ui,vi}//Transpose//Flatten;
    SMSDo[i,1,𝕒//Length];
    Ψi⊨thick Jd wGauss  SMSD[Π,𝕒,i];
    SMSExport[SMSResidualSign Ψi,p$$[i],"AddIn"->True];
    SMSDo[j,1,𝕒//Length];
    Kij⊨SMSD[Ψi,𝕒,j];
    SMSExport[Kij,s$$[i,j],"AddIn"->True]; 
    SMSEndDo[];
    SMSEndDo[];
    SMSEndDo[];
SMSWrite[];

SMTMakeDll[]

Gas pressure element

=>	surface is composed of linear segments,
=>	traction is perpendicular to the deformed surface of the domain,
=>	global unknowns are displacements of the nodes,
=>	contribution of the surface pressure to the virtual work of the problem is defined by
	Ψ=(∫_Γ^( ))tδudΓ
	where δu is variation of the displacement field, t is prescribed surface traction and Γ is the  		deformed boundary of the problem.
=>	Due to the surface pressure acting perpendicular to the deformed mesh the surface traction 		element contributes also to the global 		tangent matrix and makes the global tangent unsymmetrical.

SMSInitialize["GasPressure","Environment"->"AceFEM"];
SMSTemplate["SMSTopology"->"L1"
,"SMSSymmetricTangent"->False
,"SMSDefaultIntegrationCode"->0
,"SMSGroupDataNames"-> {"p -pressure","t -thickness"}
,"SMSDefaultData"->{0,1}];
SMSStandardModule["Tangent and residual"];
Xi⊨SMSReal[{nd$$[1,"X",1],nd$$[1,"X",2]}];
Xj⊨SMSReal[{nd$$[2,"X",1],nd$$[2,"X",2]}];
ui⊨SMSReal[{nd$$[1,"at",1],nd$$[1,"at",2]}];
uj⊨SMSReal[{nd$$[2,"at",1],nd$$[2,"at",2]}];

{p,thick}⊨Array[SMSReal[es$$["Data",#]] &, 2];
xi⊨Xi+ui;
xj⊨Xj+uj;
L⊨SMSSqrt[(xj-xi).(xj-xi)];
{cosϕ,sinϕ}⊨(xj-xi)/L;
ξ⊨SMSFictive[];
{u,v}=(1/2){ (1-ξ) ,(1+ξ) }.{ui,uj};
p={{cosϕ,-sinϕ},{sinϕ,cosϕ}}.{0,p};
Ψi⊨thick Integrate[-(L/2) p.SMSD[{u,v},Join[ui,uj]],{ξ,-1,1}];
Kij⊨SMSD[Ψi,Join[ui,uj]];
SMSExport[Ψi,p$$,"AddIn"->True];
SMSExport[Kij,s$$,"AddIn"->True];
SMSWrite[];

SMTMakeDll[]

Application definitions

FEPInitialize

GUI related definition:
=>	FEPGUIOpen
=>	FEPGUIInitialize
=>	FEPGUClose
=>	FEPMessage  (languages)
=>	FEPAbortApplication

Simulation related definition:
=>	FEPShowMesh
=>	FEPInputData
=>	FEPRunSimulation

FEPInitialize[]:=(
FEPLanguage=1;
GUISize=200;
FEPRunning=False;
FEPThickness=0.1;
FEPPressure=30;
)

FEPGUIOpen[]:=(
Column[{
Row[{Dynamic[Dynamic[FEPMessage[11]]],InputField[Dynamic[FEPThickness],Number,ImageSize->50]}]
,Row[{Dynamic[Dynamic[FEPMessage[12]]],InputField[Dynamic[FEPPressure],Number,ImageSize->50]}]
,Button[Dynamic[FEPMessage[9]],SMTInputData[];FEPShowMesh[],Method->"Queued"]
,Button[Dynamic[FEPMessage[10]],FEPRunSimulation[],Method->"Queued"]
,Button[Dynamic[FEPMessage[5]],FEPAbortApplication[FEPMessage[8],False],Method->"Queued"]
,Row[{Dynamic[FEPMessage[3]]," ",PopupMenu[Dynamic[FEPLanguage],{1->"English",2->"Slovenščina"}]}]
,Dynamic[Refresh[
Framed[Graphics[{If[FEPRunning,Red,Green],Rectangle[{0,0},{5,1},RoundingRadius->0.5]},ImageSize->{Automatic,20}],Background->Blue]
,TrackedSymbols:>{FEPRunning},UpdateInterval->0.5]
]
},Center]
)

FEPShowMesh[]:=(
FEPRunning=True;
FEPInputData[];
SMTAnalysis[];
MeshWindowContents=Column[{
Grid[{{Dynamic[FEPMessage[6]],SMTIData["NoElements"]},
{Dynamic[FEPMessage[7]],SMTIData["NoEquations"]}}]
,SMTShowMesh["BoundaryConditions"->True,ImageSize->{{ MeshSize},{ MeshSize}}
]}];
FEPRunning=False;
)

FEPInputData[]:=(
SMTInputData[];
SMTAddDomain[{"tyre","Rubber",{"E *"->1000.,"ν *"->.3}},{"gas","GasPressure",{}}];
ne=40;
SMTAddEssentialBoundary[ "Y"<-1.02&,1->0,2->0];
SMTMesh["tyre","Q1",{ne,4},Table[{(1+FEPThickness ){Cos[ϕ],Sin[ϕ]},{ Cos[ϕ],Sin[ϕ]}}, {ϕ,0,2. π,2. π/20}]//Transpose];SMTMesh["gas","L1",{ne},Table[{Cos[ϕ],Sin[ϕ]}, {ϕ,0,2. π,2. π/20}]//Reverse];
)

FEPRunSimulation[]:=Module[{},
FEPRunning=True;
FEPInputData[];
SMTAnalysis[];
undef=SMTShowMesh["BoundaryConditions"->True,"Marks"->False,"Show"->False,"Mesh"->False,"Elements"->RGBColor[1,0,0]];SMTNextStep[1,1];

(*find the shape of the empty tyre*)
Table[
SMTNextStep[0,0];
SMTDomainData["tyre","Data","Qy *"->-force];
While[SMTConvergence[10^-9,10],SMTNewtonIteration[]];
MeshWindowContents=SMTShowMesh["DeformedMesh"->True,"Show"->False,PlotRange->{{-1.5,1.5},{-1.4,1.2}}];
,{force,1,5,1}];

(*increase pressure*)
graph={};
Map[(SMTNextStep[0,0];
SMTDomainData["gas","Data","p *"->#];
While[SMTConvergence[10^-9,10],SMTNewtonIteration[]];
AppendTo[graph,MeshWindowContents=SMTShowMesh["DeformedMesh"->True,"Show"->False,"Domains"->"tyre"]];
)&,label=Join[Range[0,2,.2],Range[3,FEPPressure,FEPPressure/5.]]];
PostWindowContents=GraphicsGrid[Partition[MapThread[Show[undef,#,
DisplayFunction->Identity,PlotLabel->#2,
PlotRange->{{-1.5,1.5},{-1.4,2}}]&,{graph,label}],4]
,Spacings->0,ImageSize->450];
FEPRunning=False;
]

FEPAbortApplication[m_,force_]:=(
If[DialogInput[
Column[{m,
"",Row[{DefaultButton[DialogReturn[True]]
,"   ",CancelButton[DialogReturn[False]]}]}
,Center] || force,
WindowFloating->True,WindowSize->{200,All}]
,FEPGUClose[];
Quit[];
DialogReturn[]
]
)

FEPGUIInitialize[]:=(
MeshSize=400.;
MeshWindowContents=Graphics[Table[{Hue[t/20],Circle[{Cos[2Pi t/20],Sin[2Pi t/20]}]},{t,20}],ImageSize->{MeshSize,MeshSize}];
MeshWindow=CreateDocument[
Dynamic[MeshWindowContents]
,WindowMargins->{{GUISize+20,Automatic},{Automatic,0}}
,WindowElements->{},ShowCellBracket->False
,WindowSize->All,WindowFrameElements->{"CloseBox"}
,Saveable -> False,WindowToolbars->{}
,Background->RGBColor[1,1,1],WindowTitle->"Tyre mesh"
,WindowFrame->"Palette",ShowStringCharacters->False
,Selectable->True,Editable->False
];
PostSize=600;
PostWindowContents=Graphics[Table[{Hue[t/20],Circle[{Cos[2Pi t/20],Sin[2Pi t/20]}]},{t,20}],ImageSize->{PostSize,PostSize}];
PostWindow=CreateDocument[
Dynamic[PostWindowContents]
,WindowMargins->{{GUISize+20+MeshSize+20,Automatic},{Automatic,0}}
,WindowElements->{},ShowCellBracket->False
,WindowSize->All,WindowFrameElements->{"CloseBox"}
,Saveable -> False,WindowToolbars->{}
,Background->RGBColor[1,1,1],WindowTitle->"Tyre results"
,WindowFrame->"Palette",ShowStringCharacters->False
,Selectable->True,Editable->False
];
)

FEPGUClose[]:=(
NotebookClose[PostWindow];
NotebookClose[MeshWindow];
NotebookClose[FEPMainWindow];
)

FEPMessage[i_]:={
(*1*){"Cannot open finite element analysis module.","Napaka pri odpiranju modula za končne elemente.",}
(*2*),{"Analysis","Izračun"}
(*3*),{"Language","Jezik"}
(*4*),{"Application has been aboted","Napaka v aplikaciji"}
(*5*),{"Close application","Zapri program"}
(*6*),{"Number of elements: ","Število elementov"}
(*7*),{"Number of equations: ","Število anačb"}
(*8*),{"Close application?","Zapri program?"}
(*9*),{"Show FE mesh","Prika(Overscript[z, ⋁])i mre(Overscript[z, ⋁])o elementov"}
(*10*),{"Run simulation","Izračun"}
(*11*),{"Rubber thickness","Debelina gume"}
(*12*),{"Gus pressure","Zračni pritisk"}
}[[i,FEPLanguage]];

Create main FEMPlayer Notebook

include application definitions into notebook

<<AceFEM`;
<<AceGen`;
NotebookFind[EvaluationNotebook[],"Definitions",All,CellTags,AutoScroll->False];alldefinitions=NotebookRead[EvaluationNotebook[]];
inputdefinitions=Cases[alldefinitions,Cell[i_,"Input",___]:>i,Infinity];
Map[ToExpression,inputdefinitions];

create executable FEMPLayer notebook "PneumaticTyre.nb"

With[{inputdefinitions=inputdefinitions},
NotebookSave[

CreatePalette[
DynamicModule[{}

,Framed[Column[{"Loading Pneumatic Tyre Application","Please wait ..."},
Center,BaseStyle->{Large,Bold}],Background->Green,FrameStyle->Directive[Red,Thickness[100]],
RoundingRadius->20,FrameMargins->50]

,Initialization:>(
SetDirectory[NotebookDirectory[]];
Check[
Get["AceFEM`Player`",""];
,Print["Cannot open finite element analysis module."];Abort[];
];

Map[ToExpression,inputdefinitions];

FEPInitialize[];
FEPMainWindow=CreateDialog[
FEPGUIOpen[]
,Deployed->True,WindowTitle->"Tyre analysis",Background->LightGreen
,WindowSize->All,WindowMargins->{{0,Automatic},{Automatic,0}}
];
FEPGUIInitialize[];
NotebookClose[Notebooks["PneumaticTyre.nb"][[1]]];
)

,SynchronousInitialization->False
]
,WindowTitle->"Tyre application loader",WindowSize->All,WindowMargins->Automatic
]
,ToFileName[NotebookDirectory[],"PneumaticTyre.nb"]
]
]

PneumaticTyre bundle FEMPlayer application for Windows

PneumaticTyre application bundle consists of:
1	FEMPlayer installation (32 and 64 version for all platforms Windows/Mac/Linux)
2	dll files of solid ("Rubber.W32.dll","Rubber.W64.dll") and gas pressure ("GasPressure.W32.dll","GasPressure.W64.dll") elements 
3	main application notebook "PneumaticTyre.nb" (incorporates all application specific definitions)

Bundle size: 3.4 MB.

Application: Buckling analysis of stell members

NotebookOpen["D:\\Delo\\Shell\\Buckling\\Buckling.nb"];

Application: Masonry building structure - earthquake failure

Verification of the seismic resistance - EUROCODE

H - horizotal force
G - vertical force (weight)
d - displacement of the center of mass

Seismic resistance is governed by the factor:    ((Subscript[H, max])/(Subscript[G, total]))≳0.45

EUROCODE Method

NotebookOpen["D:\\Delo\\Gams\\Sremb\\SrembApplication.nb"];

Two step solution procedure

Conclusions

=>	The FEMPlayer package is Mathematica package for running interactive AceFEM based finite element 
	applications using Mathematica Player Pro
 
 =>	dedicated to solve only specific engineering task
 
 =>	FEMPlayer package and Mathematica Player Pro combined provide the 
 	power of Mathematica and AceFEM for the fraction of the cost
 
 =>	development of high-performance and low-cost interactive finite element based engineering applications

=>	Other fields of application:
	FEMPlayer + exact snsitivity analysis + Mathematica optimization algorithms => Large-scale engineering optimizations