Finite Element Generation

Quadrilateral Finite element

(*Element is a part of AceShare online Library*)

Line Load element

<<AceGen`AndAceFEM`;driver=2;
SMSInitialize["SurfaceLoadGConst","Language"->{"Mathematica","C"}[[driver]],"Mode"->"Optimal"];
SMSTemplate[{"MDriver","CDriver"}[[driver]],
"SMSTopology"->"L1",
"SMSDefaultIntegrationCode"->0,
"SMSNoNodes"->3,
"SMSDOFGlobal"->{2,2,1},
"SMSSegments"->{{1,2}},
"SMSAdditionalNodes"->Hold[{{0,0}}&],
"SMSNodeID"->{"D","D","load -F -AL -S"},
"SMSSymmetricTangent"->False
];
Map[(
SMSStandardModule[#];
dof⊨SMSReal[Table[nd$$[i,"at",j],{i,SMSNoNodes},{j, SMSDOFGlobal[[i]]}],-0.1,0.1];
X⊨SMSReal[Table[nd$$[i,"X",j],{i,SMSNoNodes}, {j,SMSNoDimensions}],-100,100];
𝕒t⊨Flatten[dof];
ri⊨X[[1]];
rj⊨X[[2]];
ui⊨dof[[1]];
uj⊨dof[[2]];
λ⊨dof[[3,1]];
SMSGroupDataNames = {"X pressure","Y pressure","thickness"};
{pX,pY,thick}⊨Array[SMSReal[es$$["Data",#]] &, 3];
L⊨SMSSqrt[(rj-ri).(rj-ri)];
ξ⊨SMSFictive[];
{u,v}=(1/2){ (1-ξ) ,(1+ξ) }.{ui,uj};
p=λ  {pX,pY};
Ψi⊨Join[thick Integrate[-(L/2) p.SMSD[{u,v},Join[ui,uj]],{ξ,-1,1}],{0}];
Switch[#
,"Tangent and residual",
SMSExport[Ψi,p$$,"AddIn"->True];
K=SMSD[Ψi,𝕒t];
SMSExport[K,s$$,"AddIn"->True];
,"Sensitivity pseudo-load",
SensIndex⊨SMSInteger[idata$$["SensIndex"]];
SensType⊨SMSInteger[es$$["SensType",SensIndex]];
SensTypeIndex⊨SMSInteger[es$$["SensTypeIndex",SensIndex]];
ϕsen⊨SMSFictive[];
sXYZ⊨SMSReal[Table[nd$$[i,"sX",SensTypeIndex,j],{i,SMSNoNodes}, {j,SMSNoDimensions}]];
MapThread[SMSDefineDerivative[#1,ϕsen,#2]&,{X//Flatten,sXYZ//Flatten}];
Si⊨SMSD[Ψi,ϕsen];
SMSExport[SMSResidualSign Si,p$$,"AddIn"->True];
];
)&,{"Tangent and residual","Sensitivity pseudo-load"}];
SMSStandardModule["Dependent sensitivity"];
SMSWrite[];

Element for including the additional equation

<<AceGen`AndAceFEM`;driver=2;
SMSInitialize["PrescribedDispl","Language"->{"Mathematica","C"}[[driver]],"Mode"->"Optimal"];
SMSTemplate[{"MDriver","CDriver"}[[driver]],
"SMSTopology"->"LX",
"SMSDefaultIntegrationCode"->0,
"SMSNoNodes"->2,
"SMSDOFGlobal"->{2,1},
"SMSSegments"->{},
"SMSAdditionalNodes"->Hold[{{0,0}}&],
"SMSNodeID"->{"D","load -F -AL -S"},
"SMSSymmetricTangent"->False
];
Map[(
SMSStandardModule[#];
dof⊨SMSReal[Table[nd$$[i,"at",j],{i,SMSNoNodes},{j, SMSDOFGlobal[[i]]}],-0.1,0.1];
X⊨SMSReal[Table[nd$$[i,"X",j],{i,SMSNoNodes}, {j,SMSNoDimensions}],-100,100];
SMSGroupDataNames = {"Prescribed displ"};
{up}⊨Array[SMSReal[es$$["Data",#]] &, 1];
𝕒t⊨Flatten[dof];
λ⊨dof[[2,1]];
λu⊨SMSReal[rdata$$["Multiplier"]];
Ψi⊨{0,0,dof[[1,2]]-(λu up)};
Switch[#
,"Tangent and residual",
SMSExport[Ψi,p$$,"AddIn"->True];
K=SMSD[Ψi,𝕒t];
SMSExport[K,s$$,"AddIn"->True];
,"Sensitivity pseudo-load",
SensIndex⊨SMSInteger[idata$$["SensIndex"]];
SensType⊨SMSInteger[es$$["SensType",SensIndex]];
SensTypeIndex⊨SMSInteger[es$$["SensTypeIndex",SensIndex]];
ϕsen⊨SMSFictive[];
sXYZ⊨SMSReal[Table[nd$$[i,"sX",SensTypeIndex,j],{i,SMSNoNodes}, {j,SMSNoDimensions}]];
MapThread[SMSDefineDerivative[#1,ϕsen,#2]&,{X//Flatten,sXYZ//Flatten}];
Si⊨SMSD[Ψi,ϕsen];
SMSExport[SMSResidualSign Si,p$$,"AddIn"->True];
];
)&,{"Tangent and residual","Sensitivity pseudo-load"}];
SMSStandardModule["Dependent sensitivity"];
SMSWrite[];

Data

<<AceFEM`;
L=200; (*Length*)
H=2; (* Width *)
T=1; (* Thickness *)
fy=23.5;(* yield stress *)
Em=21000;(* elastic modulus *)
Nλ=20;(* Number of considered shapes *)
e0=L/250; (* Amplitude of the equivalent geometrical imperfection *)

BoundaryConditions={{ "X"==0&,0,0}};
AlternativeBoundaryConditions={{ "X"==0&&"Y"==H/2&,0,0},{ "X"==L&&"Y"==H/2&,,0}};

nx=200;(* number of finite elements in x-direction*)
ny=2; (* number of finite elements in y-direction*)

domains={{"Domain1","OL:SEPSQ1DIJCQ1TSNHookeAMisses",{Em,0.3,fy,T}},{"Load","SurfaceLoadGConst",{-0.001,-0.001/1000,1.}},{"PrescDispl","PrescribedDispl",{-2}}};

Generation of the shape base

Buckling shapes

(* -------------------------------------------- Finite element mesh  -------------------------------------------- *)
SMTInputData["CDriver"];
SMTAddDomain[{domains[[1]]}];
SMTAddEssentialBoundary[BoundaryConditions];
SMTAddNaturalBoundary[{{ "X"==L&&"Y"==H/2&,-1,0}}]; (* Axial load for buckling *)
SMTMesh["Domain1","Q1",{nx,ny},{{{0,0},{L,0}},{{0,H},{L,H}}}];
SMTAnalysis["Output"->"uklonobl.out","Solver"->{5,{11,0,0,9,1},{}},"OptimizeDll"-> True];

(* -------------------------------------------- Buckling Shape evaluation  -------------------------------------------- *)
K1=SMTData["TangentMatrix"];
SMTNextStep[1,1];While[SMTConvergence[10^-8,10],SMTNewtonIteration[];];
K2=SMTData["TangentMatrix"];
{BucklLoad,BucklShape}=Eigensystem[{K2,K1},-Nλ,Method -> {"Arnoldi","MaxIterations"->10000}];
BucklLoad=Map[(1/(1-#))&,BucklLoad];
ΓA=Map[(normed=#1/Max[Abs[#1]];Transpose[Array[Map[If[#<1,0,normed[[#]]]&,Take[SMTNodeData[#,"DOF"]+1,SMTSpatialDimension]]&,SMTIData["NoNodes"]]])&,BucklShape//Reverse];
Map[Print[Grid[{{#." buckling shape",SMTShowMesh["DeformedMesh"->ΓA[[#]],"ImageSize"->{{400},{400}},"BoundaryConditions"->True,ViewPoint->{0.487, 1.861, 1.579},"Scale"->10]}}]]&,Range[Nλ]];
Γ=Map[(#//Transpose)&,ΓA];Γ>>Gama.dat; (*Only ΓA(buckling shapes) are taken for the shape base*)

Determination of the most unfavorable imperfection

Definition of the Optimization Function

λm:=λm0+(∑_(k=1)^(Nλ))(δλδα[[k]]*ToExpression[αtxt][[k]]);
LPcoefficientF:=Coefficient[λm,ToExpression[αtxt]];

Definition of Constraints

LPC1AllDOF:=Flatten[Map[(#//Transpose)&,(Map[(Γ[[#,points1]])&,Range[Length[Γ]]])//Transpose],1];
LPC1:=Table[LPC1AllDOF[[i+1]],{i,1,Length[LPC1AllDOF],2}] ;(*center line, y-direction*)
LPC2AllDOF:=Flatten[Map[(#//Transpose)&,(Map[(Γ[[#,points2]])&,Range[Length[Γ]]])//Transpose],1];
LPC2:=Table[LPC2AllDOF[[i]],{i,1,Length[LPC2AllDOF],2}];(*upper boundary line, x-direction*)
LPC3AllDOF:=Flatten[Map[(#//Transpose)&,(Map[(Γ[[#,points3]])&,Range[Length[Γ]]])//Transpose],1];
LPC3:=Table[LPC3AllDOF[[i]],{i,1,Length[LPC3AllDOF],2}];(*lower boundary line, x-direction*)

Iterative optimization process

Γ=<<Gama.dat;
Clear["α*"];
αtxt=MapThread[(ToString[#])&,{Array[ToExpression["α"<>ToString[#]]&,Length[Γ]]}];
αsymb=ToExpression[αtxt];

e01=e0;e02=e0/20;e03=e0/20;
Tol=(10^(-6)); (*Tolerance for Norm[Δα]*)
Xtotal=-Γ[[1]]*e0; (*aproximation for the first iteration*)

XtotalHistory={Xtotal};
αHistory={};
Sumα=Table[0,{Length[αtxt]}];
Δα=Table[1,{Length[Γ]}];

While[Norm[Δα]>Tol,

(*----------Direct and Sensitivity analysis--------------*)
SMTInputData["CDriver"];
SMTAddDomain[domains];
SMTAddEssentialBoundary[BoundaryConditions];
SMTMesh["Domain1","Q1",{nx,ny},{{{0,0},{L,0}},{{0,H},{L,H}}}];SMTMesh["Load","L1",{ny},{{L,0},{L,H}}];
SMTAddElement["PrescDispl",{{L,0}}];SMTAddSensitivity[Table[{"α"<>ToString[i],0,Map[domains[[#1,1]]->{2,i}&,Range[Length[domains]]]}//Flatten,{i,Length[Γ]}]];
SMTAnalysis["Output"->"Analysis.out","Solver"->{5,{11,0,0,9,1},{}},"OptimizeDll"-> True];
δXδϕ=Map[Flatten,Transpose[Map[Partition[Flatten[Append[#,Table[0,{1},{2}]]],2]&,Γ]]];
points1=SMTFindNodes[{"Find","Y"==H/2&}];points2=SMTFindNodes[{"Find","Y"==H&}];points3=SMTFindNodes[{"Find","Y"==0&}];

Print[Style["Starting direct and sensitivity analysis...   ",Blue //Darker,18,Bold]];
wnode1=SMTFindNodes[("X"==L&&"Y"==0)&][[1]];
λnode=SMTFindNodes["ID"=="load"&][[1]];
Print["Initial imperfection (scale f. = 10) :",SMTShowMesh["DeformedMesh"->(Append[Xtotal,{0,0}]//Transpose),"Scale"->10,"BoundaryConditions"->True]];
SMTNodeData["X",SMTNodeData["X"]+Append[Xtotal,{0,0}]];SMTUpdatePostData[];SMTNodeData["Data",N[δXδϕ]];
graph={{0,0}};SensFromSteps={Table[0,{Length[Γ]}]};

SMTNextStep[1,0.01];
While[
While[step=SMTConvergence[10^-8,10,{"Adaptive",8,0.00000000001,1,1.}], SMTNewtonIteration[];"Alternate"->Automatic ];
If[step[[4]]==="MinBound",SMTStatusReport[" Error: Δλ < Δλmin"];SMTStepBack[];Break[];];
If[Not[step[[1]]],SMTSensitivity[];AppendTo[SensFromSteps,SMTNodeData[λnode,"st"]];AppendTo[graph,{-SMTNodeData[wnode1,"at"][[2]],SMTNodeData[λnode,"at"][[1]]}];
SMTShowMesh["DeformedMesh"->True,"Scale"->10,"Show"->"Window"];
If[graph[[-1,2]]<graph[[-2,2]],SMTStepBack[];graph=Delete[graph,-1];SensFromSteps=Delete[SensFromSteps,-1];Print[Style["limit load reached !",15,Bold,RGBColor[1,0,0]]];Break[]];
If[SMTIData["DiagonalSign"]>0,SMTStepBack[];Print[Style["Bifurcation point!",15,Bold,RGBColor[1,0,0]]];Break[]];
If[SMTPost["Plastic deformation"];SMTRData["PostMaxValue"]>(50 fy/Em),SMTStepBack[];
Print[Style["Maximum plastic deformation reached",15,Bold,RGBColor[1,0,0]]];Break[];];
];step[[3]] ,If[step[[1]],SMTStepBack[]];SMTNextStep[1,step[[2]]];
];
Print[ListLinePlot[graph,PlotRange->All,PlotStyle->{{Hue[0.99],Thickness[0.002]},{Hue[0.6],Thickness[0.001]}},AxesLabel->{"displacement [cm]","λ"},PlotMarkers->{Automatic,8}]];
δλδα=SensFromSteps[[-1]];
Print["δλδα= ",δλδα];
(*---------------------------------------Optimization-------------------------------------------------*)
Print[Style["Starting Optimization...   ",Red //Darker//Darker,18,Bold]];
λm0=graph[[-1,2]];
sol=LinearProgramming[
LPcoefficientF,
Flatten[{LPC1,LPC1,LPC2,LPC2,LPC3,LPC3},1],
Flatten[{
Map[{(e01-#[[2]]),-1}&,(Xtotal[[points1]])],Map[{(-e01-#[[2]]),1}&,(Xtotal[[points1]])],
Map[{(e02-#[[1]]),-1}&,(Xtotal[[points2]])],Map[{(-e02-#[[1]]),1}&,(Xtotal[[points2]])],
Map[{(e03-#[[1]]),-1}&,(Xtotal[[points3]])],Map[{(-e03-#[[1]]),1}&,(Xtotal[[points3]])]
},1],
Table[{-∞,∞},{Length[Γ]}],Method->{"InteriorPoint","Automatic"}[[2]]];
AppendTo[αHistory,Δα=αsymb/.(MapThread[Rule[#1,#2]&,{ToExpression[αtxt],sol}])];Sumα=Sumα+Δα;ΔXtotal=Δα.Γ;AppendTo[XtotalHistory,Xtotal=Xtotal+ΔXtotal];

Print[Style["Maximum component of Max[Abs[Xtotal]] = ",10,Blue//Darker],(Xtotal//Abs//Max)];
Print[Style["Norm[Δα] = ",10,Blue//Darker],Norm[Δα]];
];

If[Norm[Δα]<Tol,Print[Style["Convergence reached",Green//Darker//Darker,20,Bold]];]Print["The most unfavorable imperfection (scale f. = 10):",SMTShowMesh["DeformedMesh"->(Append[(XtotalHistory[[-1]]),{0,0}]//Transpose),"Scale"->10,"BoundaryConditions"->True,Ticks->{Automatic,Automatic},ViewPoint->{0.129, 1.875, 1.617},"Elements"->False,Axes->True,AxesStyle->{Directive[Black,20],Directive[Black,20]},"ImageSize"->{{800},{400}}]];

RESULTS

Most unfavorable imperfection (X) (scale factor = 10)

For[i=1,i≤Length[XtotalHistory],i++,
Print[i,". iteration",SMTShowMesh["DeformedMesh"->(Append[(XtotalHistory[[i]]),{0,0}]//Transpose),"Scale"->10,"BoundaryConditions"->True]];]

Evaluated most unfavorable Imperfection increament (Δα.Γ)  (scale factor = 20)

For[i=1,i≤Length[αHistory],i++,
Print[i,".iteration ΔX",SMTShowMesh["DeformedMesh"->(Append[((αHistory)[[i]]).Γ,{0,0}]//Transpose),"Scale"->30,"BoundaryConditions"->True]];]

Contribution of base shapes to the evaluated most unfavorable imperfection shape

Needs["BarCharts`"];
BarChart[Sumα,BarStyle->Darker[LightBlue,0.3],AxesStyle->{Directive[Red//Darker,12],Directive[Blue//Darker,12]},Ticks->Automatic,AxesLabel->{"Base shape","αi"}]

Convergence plot

ListLinePlot[{Table[i,{i,Length[αHistory]}],Map[Norm[#]&,αHistory]}//Transpose]