DefineConstant[
visu = {1, Choices{0, 1}, AutoCheck 0,
Name "Input/Solver/Visu", Label "Real-time visualization"}
];
Group {
Left = Region[{1}];
Top = Region[{2}];
Right = Region[{3}];
Domain_Mecha = Region[{100}]; // Mechanical domain
}
Flag_ExternalForce = 0;
Function {
DefineFunction[ FT_F, PK2, C_Tot_xx, C_Tot_xy, C_Tot_xz, C_Tot_yx, C_Tot_yy, C_Tot_yz, C_Tot_zx, C_Tot_zy, C_Tot_zz, E_g, nu_g];
// 1. Parameters of the problem
//=============================
E = 771.0e6;
nu = 0.29;
//Lame parameters
//===============
lambda = E/(2.0 * (1.0 + nu)); // First Lame parameter
mu = (nu * E)/((1.0 + nu) * (1.0 - 2 * nu)); // Second Lame parameter
//Other material properties
//=========================
rho_0 = 2700; // for Alu
cp_0 = 897; // for Alu
// Parameters used during resolution
//==================================
t_0 = 0.0;
t_end = 4.0e-0;
period = (t_end - t_0);
Freq = 1.0/period;
tol_abs[] = 1e-8;
tol_rel[] = 1e-8;
iter_max = 50;
dt_max[] = $var_Dt_Max;
dtime[] = $var_DTime;
v_TS[] = $var_TS;
num_steps = 32;
theta_value = 1.0;
// Definition of parameters of the Dirichlet BCs, the dynamic problem and the external force
//==========================================================================================
force_amplitude = 1.8e8;
force_ext[] = Vector[0., force_amplitude, 0.0];
}
//Include "smp_stent_KMat_Elasticity.pro";
Function {
// 3. Kinetic information : second Piola Kirchhoff [PK2] stresses
//===============================================================
// Green Lagrange strain measure GL = E = 0.5 * (C - I) and the 2nd Piola-Kirchoff stress
//=======================================================================================
FT_F[] = Transpose[$1] * $1; // C = (F^T * F)
GL[] = 0.5 * (FT_F[$1] - TensorDiag[1.0, 1.0, 1.0]);
GL_Trace[] = (CompXX[GL[$1]] + CompYY[GL[#1]] + CompZZ[GL[#1]]) ;
PK2[] = (lambda * GL_Trace[$1] * TensorDiag[1.0, 1.0, 1.0] + 2 * mu * GL[$1]);
// 4. The nonlinear right hand side
//=================================
//RHS[] = PK2[$1] * Transpose[$1]; // The first Piola Kirchhoff : P = S * F^T
RHS[] = $1 * PK2[$1]; // The first Piola Kirchhoff : P = F S
// x
P~{1}[] = Vector[CompXX[RHS[$1#1]], CompXY[RHS[#1]], CompXZ[RHS[#1]] ];
// y
P~{2}[] = Vector[CompYX[RHS[$1#1]], CompYY[RHS[#1]], CompYZ[RHS[#1]] ];
// z
P~{3}[] = Vector[CompZX[RHS[$1#1]], CompZY[RHS[#1]], CompZZ[RHS[#1]] ];
// 5. The stiffness matrix : C_Tot_ij = C_Mat_ij + C_Geo_ij
//==========================================================
// 5.1. Geometric stiffness "C_Geo_ij" (see formula (4.69) of "Wriggers, P.,
// 2008. Nonlinear finite element methods. Springer Science & Business
// Media.")
//==========================================================
// xx
C_Geo~{1}~{1}[] = PK2[$1];
// yx
C_Geo~{2}~{1}[] = Tensor[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
// zx
C_Geo~{3}~{1}[] = Tensor[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
// xy
C_Geo~{1}~{2}[] = Tensor[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
// yy
C_Geo~{2}~{2}[] = PK2[$1];
// zy
C_Geo~{3}~{2}[] = Tensor[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
// xz
C_Geo~{1}~{3}[] = Tensor[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
// yz
C_Geo~{2}~{3}[] = Tensor[ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
// zz
C_Geo~{3}~{3}[] = PK2[$1];
// The material stiffness
//========================
// xx
C_Mat_mu~{1}~{1}[] = mu * Tensor[ 2 * CompXX[$1#1] * CompXX[#1] + CompXY[#1] * CompXY[#1] + CompXZ[#1] * CompXZ[#1],
CompXX[#1] * CompXY[#1],
CompXX[#1] * CompXZ[#1],
CompXY[#1] * CompXX[#1],
CompXX[#1] * CompXX[#1] + 2 * CompXY[#1] * CompXY[#1] + CompXZ[#1] * CompXZ[#1],
CompXY[#1] * CompXZ[#1],
CompXZ[#1] * CompXX[#1],
CompXZ[#1] * CompXY[#1],
CompXX[#1] * CompXX[#1] + CompXY[#1] * CompXY[#1] + 2 * CompXZ[#1] * CompXZ[#1] ];
// yx
C_Mat_mu~{2}~{1}[] = mu * Tensor[ 2 * CompYX[$1#1] * CompXX[#1] + CompYY[#1] * CompXY[#1] + CompYZ[#1] * CompXZ[#1],
CompYX[#1] * CompXY[#1],
CompYX[#1] * CompXZ[#1],
CompYY[#1] * CompXX[#1],
CompYX[#1] * CompXX[#1] + 2 * CompYY[#1] * CompXY[#1] + CompYZ[#1] * CompXZ[#1],
CompYY[#1] * CompXZ[#1],
CompYZ[#1] * CompXX[#1],
CompYZ[#1] * CompXY[#1],
CompYX[#1] * CompXX[#1] + CompYY[#1] * CompXY[#1] + 2 * CompYZ[#1] * CompXZ[#1] ];
// zx
C_Mat_mu~{3}~{1}[] = mu * Tensor[ 2 * CompZX[$1#1] * CompXX[#1] + CompZY[#1] * CompXY[#1] + CompZZ[#1] * CompXZ[#1],
CompZX[#1] * CompXY[#1],
CompZX[#1] * CompXZ[#1],
CompZY[#1] * CompXX[#1],
CompZX[#1] * CompXX[#1] + 2 * CompZY[#1] * CompXY[#1] + CompZZ[#1] * CompXZ[#1],
CompZY[#1] * CompXZ[#1],
CompZZ[#1] * CompXX[#1],
CompZZ[#1] * CompXY[#1],
CompZX[#1] * CompXX[#1] + CompZY[#1] * CompXY[#1] + 2 * CompZZ[#1] * CompXZ[#1] ];
// xy
C_Mat_mu~{1}~{2}[] = mu * Tensor[ 2 * CompXX[$1#1] * CompYX[#1] + CompXY[#1] * CompYY[#1] + CompXZ[#1] * CompYZ[#1],
CompXX[#1] * CompYY[#1],
CompXX[#1] * CompYZ[#1],
CompXY[#1] * CompYX[#1],
CompXX[#1] * CompYX[#1] + 2 * CompXY[#1] * CompYY[#1] + CompXZ[#1] * CompYZ[#1],
CompXY[#1] * CompYZ[#1],
CompXZ[#1] * CompYX[#1],
CompXZ[#1] * CompYY[#1],
CompXX[#1] * CompYX[#1] + CompXY[#1] * CompYY[#1] + 2 * CompXZ[#1] * CompYZ[#1] ];
// yy
C_Mat_mu~{2}~{2}[] = mu * Tensor[ 2 * CompYX[$1#1] * CompYX[#1] + CompYY[#1] * CompYY[#1] + CompYZ[#1] * CompYZ[#1],
CompYX[#1] * CompYY[#1],
CompYX[#1] * CompYZ[#1],
CompYY[#1] * CompYX[#1],
CompYX[#1] * CompYX[#1] + 2 * CompYY[#1] * CompYY[#1] + CompYZ[#1] * CompYZ[#1],
CompYY[#1] * CompYZ[#1],
CompYZ[#1] * CompYX[#1],
CompYZ[#1] * CompYY[#1],
CompYX[#1] * CompYX[#1] + CompYY[#1] * CompYY[#1] + 2 * CompYZ[#1] * CompYZ[#1] ];
// zy
C_Mat_mu~{3}~{2}[] = mu * Tensor[ 2 * CompZX[$1#1] * CompYX[#1] + CompZY[#1] * CompYY[#1] + CompZZ[#1] * CompYZ[#1],
CompZX[#1] * CompYY[#1],
CompZX[#1] * CompYZ[#1],
CompZY[#1] * CompYX[#1],
CompZX[#1] * CompYX[#1] + 2 * CompZY[#1] * CompYY[#1] + CompZZ[#1] * CompYZ[#1],
CompZY[#1] * CompYZ[#1],
CompZZ[#1] * CompYX[#1],
CompZZ[#1] * CompYY[#1],
CompZX[#1] * CompYX[#1] + CompZY[#1] * CompYY[#1] + 2 * CompZZ[#1] * CompYZ[#1] ];
// xz
C_Mat_mu~{1}~{3}[] = mu * Tensor[ 2 * CompXX[$1#1] * CompZX[#1] + CompXY[#1] * CompZY[#1] + CompXZ[#1] * CompZZ[#1],
CompXX[#1] * CompZY[#1],
CompXX[#1] * CompZZ[#1],
CompXY[#1] * CompZX[#1],
CompXX[#1] * CompZX[#1] + 2 * CompXY[#1] * CompZY[#1] + CompXZ[#1] * CompZZ[#1],
CompXY[#1] * CompZZ[#1],
CompXZ[#1] * CompZX[#1],
CompXZ[#1] * CompZY[#1],
CompXX[#1] * CompZX[#1] + CompXY[#1] * CompZY[#1] + 2 * CompXZ[#1] * CompZZ[#1] ];
// yz
C_Mat_mu~{2}~{3}[] = mu * Tensor[ 2 * CompYX[$1#1] * CompZX[#1] + CompYY[#1] * CompZY[#1] + CompYZ[#1] * CompZZ[#1],
CompYX[#1] * CompZY[#1],
CompYX[#1] * CompZZ[#1],
CompYY[#1] * CompZX[#1],
CompYX[#1] * CompZX[#1] + 2 * CompYY[#1] * CompZY[#1] + CompYZ[#1] * CompZZ[#1],
CompYY[#1] * CompZZ[#1],
CompYZ[#1] * CompZX[#1],
CompYZ[#1] * CompZY[#1],
CompYX[#1] * CompZX[#1] + CompYY[#1] * CompZY[#1] + 2 * CompYZ[#1] * CompZZ[#1] ];
// zz
C_Mat_mu~{3}~{3}[] = mu * Tensor[ 2 * CompZX[$1#1] * CompZX[#1] + CompZY[#1] * CompZY[#1] + CompZZ[#1] * CompZZ[#1],
CompZX[#1] * CompZY[#1],
CompZX[#1] * CompZZ[#1],
CompZY[#1] * CompZX[#1],
CompZX[#1] * CompZX[#1] + 2 * CompZY[#1] * CompZY[#1] + CompZZ[#1] * CompZZ[#1],
CompZY[#1] * CompZZ[#1],
CompZZ[#1] * CompZX[#1],
CompZZ[#1] * CompZY[#1],
CompZX[#1] * CompZX[#1] + CompZY[#1] * CompZY[#1] + 2 * CompZZ[#1] * CompZZ[#1] ];
// 5.2.2. Material stiffness - the part resulting from "lambda" : C_Mat_lambda_ij
//===============================================================================
// xx
C_Mat_lambda~{1}~{1}[] = lambda * Tensor[ CompXX[$1#1] * CompXX[#1], CompXX[#1] * CompXY[#1], CompXX[#1] * CompXZ[#1],
CompXY[#1] * CompXX[#1], CompXY[#1] * CompXY[#1], CompXY[#1] * CompXZ[#1],
CompXZ[#1] * CompXX[#1], CompXZ[#1] * CompXY[#1], CompXZ[#1] * CompXZ[#1]];
// yx
C_Mat_lambda~{2}~{1}[] = lambda * Tensor[ CompYX[$1#1] * CompXX[#1], CompYX[#1] * CompXY[#1], CompYX[#1] * CompXZ[#1],
CompYY[#1] * CompXX[#1], CompYY[#1] * CompXY[#1], CompYY[#1] * CompXZ[#1],
CompYZ[#1] * CompXX[#1], CompYZ[#1] * CompXY[#1], CompYZ[#1] * CompXZ[#1]];
// zx
C_Mat_lambda~{3}~{1}[] = lambda * Tensor[ CompZX[$1#1] * CompXX[#1], CompZX[#1] * CompXY[#1], CompZX[#1] * CompXZ[#1],
CompZY[#1] * CompXX[#1], CompZY[#1] * CompXY[#1], CompZY[#1] * CompXZ[#1],
CompZZ[#1] * CompXX[#1], CompZZ[#1] * CompXY[#1], CompZZ[#1] * CompXZ[#1]];
// xy
C_Mat_lambda~{1}~{2}[] = lambda * Tensor[ CompXX[$1#1] * CompYX[#1], CompXX[#1] * CompYY[#1], CompXX[#1] * CompYZ[#1],
CompXY[#1] * CompYX[#1], CompXY[#1] * CompYY[#1], CompXY[#1] * CompYZ[#1],
CompXZ[#1] * CompYX[#1], CompXZ[#1] * CompYY[#1], CompXZ[#1] * CompYZ[#1]];
// yy
C_Mat_lambda~{2}~{2}[] = lambda * Tensor[ CompYX[$1#1] * CompYX[#1], CompYX[#1] * CompYY[#1], CompYX[#1] * CompYZ[#1],
CompYY[#1] * CompYX[#1], CompYY[#1] * CompYY[#1], CompYY[#1] * CompYZ[#1],
CompYZ[#1] * CompYX[#1], CompYZ[#1] * CompYY[#1], CompYZ[#1] * CompYZ[#1]];
// zy
C_Mat_lambda~{3}~{2}[] = lambda * Tensor[ CompZX[$1#1] * CompYX[#1], CompZX[#1] * CompYY[#1], CompZX[#1] * CompYZ[#1],
CompZY[#1] * CompYX[#1], CompZY[#1] * CompYY[#1], CompZY[#1] * CompYZ[#1],
CompZZ[#1] * CompYX[#1], CompZZ[#1] * CompYY[#1], CompZZ[#1] * CompYZ[#1]];
// xz
C_Mat_lambda~{1}~{3}[] = lambda * Tensor[ CompXX[$1#1] * CompZX[#1], CompXX[#1] * CompZY[#1], CompXX[#1] * CompZZ[#1],
CompXY[#1] * CompZX[#1], CompXY[#1] * CompZY[#1], CompXY[#1] * CompZZ[#1],
CompXZ[#1] * CompZX[#1], CompXZ[#1] * CompZY[#1], CompXZ[#1] * CompZZ[#1]];
// yz
C_Mat_lambda~{2}~{3}[] = lambda * Tensor[ CompYX[$1#1] * CompZX[#1], CompYX[#1] * CompZY[#1], CompYX[#1] * CompZZ[#1],
CompYY[#1] * CompZX[#1], CompYY[#1] * CompZY[#1], CompYY[#1] * CompZZ[#1],
CompYZ[#1] * CompZX[#1], CompYZ[#1] * CompZY[#1], CompYZ[#1] * CompZZ[#1]];
// zz
C_Mat_lambda~{3}~{3}[] = lambda * Tensor[ CompZX[$1#1] * CompZX[#1], CompZX[#1] * CompZY[#1], CompZX[#1] * CompZZ[#1],
CompZY[#1] * CompZX[#1], CompZY[#1] * CompZY[#1], CompZY[#1] * CompZZ[#1],
CompZZ[#1] * CompZX[#1], CompZZ[#1] * CompZY[#1], CompZZ[#1] * CompZZ[#1]];
// 5.2. Total material stiffness : C_Mat_ij = C_Mat_lambda_ij + C_Mat_mu_ij
//==========================================================================
For i In {1:3}
For j In {1:3}
C_Mat~{i}~{j}[] = C_Mat_lambda~{i}~{j}[$1#1] + Transpose[C_Mat_mu~{i}~{j}[#1]] ;
EndFor
EndFor
}
Function {
// 5.3 Total stiffness : C_Mat_ij + C_Geo_ij
//===========================================
For i In {1:3}
For j In {1:3}
C_Tot~{i}~{j}[] = C_Mat~{i}~{j}[$1] + C_Geo~{i}~{j}[$1];
EndFor
EndFor
}
Constraint {
{ Name Displacement_1 ; // x
Case {
{ Region Left; Value 0.0 ; }
}
}
{ Name Displacement_2 ; // y
Case {
{ Region Left; Value 0.0 ; }
}
}
{ Name Displacement_3 ; // z
Case {
{ Region Left; Value 0.0 ; }
}
}
}
Include "Jacobian_Lib.pro"
Include "Integration_Lib.pro"
Include "Total_Lagrangian.pro"
Resolution {
{ Name Mec_SMP ;
System {
{ Name Sys_Mec ; NameOfFormulation Total_Lagrangian ; }
}
Operation {
SetGlobalSolverOptions["-petsc_prealloc 400"];
CreateDir["res/"];
Evaluate[ $syscount = 0 ];
InitSolution[Sys_Mec];
Evaluate[ $var_DTime = period/(num_steps)];
Evaluate[ $var_Dt_Max = period/(num_steps)];
TimeLoopTheta[t_0, 1.0 * t_end, dtime[], theta_value]{ //TimeLoopNewmark[t_0, t_end, dtime[], beta, gamma] {
Generate[Sys_Mec]; Solve[Sys_Mec]; Evaluate[ $syscount = $syscount + 1 ];
Generate[Sys_Mec]; GetResidual[Sys_Mec, $res0]; Evaluate[ $res = $res0, $iter = 0 ];
Print[{$iter, $res, $res / $res0}, Format "Residual %03g: abs %14.12e rel %14.12e"];
While[$res > tol_abs[] && $res / $res0 > tol_rel[] && $res / $res0 <= 1 && $iter < iter_max]{
Solve[Sys_Mec]; Evaluate[ $syscount = $syscount + 1 ];
Generate[Sys_Mec]; GetResidual[Sys_Mec, $res]; Evaluate[ $iter = $iter + 1 ];
Print[{$iter, $res, $res / $res0}, Format "Residual %03g: abs %14.12e rel %14.12e"];
}
Test[ ($iter < iter_max) && ($res / $res0 <= 1) ]{
SaveSolution[Sys_Mec];
Test[ GetNumberRunTime[visu]{"Input/Solver/Visu"} ]{
PostOperation[Mec];
}
Test[ ($iter < iter_max / 10) && ($DTime < dt_max[]) ]{
Evaluate[ $dt_new = Min[$DTime * 2.0, dt_max[]] ];
Print[{$dt_new}, Format "*** Fast convergence: increasing time step to %g"];
SetDTime[$dt_new];
}
}
{
Evaluate[ $dt_new = $DTime/2.0 ];
Print[{$iter, $dt_new}, Format "*** Non convergence (iter %g): recomputing with reduced step %g"];
SetTime[$Time - $DTime];
SetTimeStep[$TimeStep - 1];
RemoveLastSolution[Sys_Mec];
SetDTime[$dt_new];
Evaluate[ $var_DTime = $dt_new ];
}
}
Print[{$syscount}, Format "Total number of linear systems solved: %g"];
}
}
}
PostOperation {
{ Name Mec ; NameOfPostProcessing Total_Lagrangian;
Operation {
Print[ u, OnElementsOf Domain_Mecha, File "res/u.pos"] ;
Print[ PK2_xx, OnElementsOf Domain_Mecha, File "res/PK2_xx.pos"] ;
Print[ PK2_xy, OnElementsOf Domain_Mecha, File "res/PK2_xy.pos"] ;
Print[ PK2_xz, OnElementsOf Domain_Mecha, File "res/PK2_xz.pos"] ;
Print[ PK2_yx, OnElementsOf Domain_Mecha, File "res/PK2_yx.pos"] ;
Print[ PK2_yy, OnElementsOf Domain_Mecha, File "res/PK2_yy.pos"] ;
Print[ PK2_yz, OnElementsOf Domain_Mecha, File "res/PK2_yz.pos"] ;
Print[ PK2_zx, OnElementsOf Domain_Mecha, File "res/PK2_zx.pos"] ;
Print[ PK2_zy, OnElementsOf Domain_Mecha, File "res/PK2_zy.pos"] ;
Print[ PK2_zz, OnElementsOf Domain_Mecha, File "res/PK2_zz.pos"] ;
Print[ F_tot_xx, OnElementsOf Domain_Mecha, File "res/F_tot_xx.pos"] ;
Print[ F_tot_xy, OnElementsOf Domain_Mecha, File "res/F_tot_xy.pos"] ;
Print[ F_tot_xz, OnElementsOf Domain_Mecha, File "res/F_tot_xz.pos"] ;
Print[ F_tot_yx, OnElementsOf Domain_Mecha, File "res/F_tot_yx.pos"] ;
Print[ F_tot_yy, OnElementsOf Domain_Mecha, File "res/F_tot_yy.pos"] ;
Print[ F_tot_yz, OnElementsOf Domain_Mecha, File "res/F_tot_yz.pos"] ;
Print[ F_tot_zx, OnElementsOf Domain_Mecha, File "res/F_tot_zx.pos"] ;
Print[ F_tot_zy, OnElementsOf Domain_Mecha, File "res/F_tot_zy.pos"] ;
Print[ F_tot_zz, OnElementsOf Domain_Mecha, File "res/F_tot_zz.pos"] ;
}
}
}
DefineConstant[
R_ = {"Mec_SMP", Name "GetDP/1ResolutionChoices", Visible 0},
C_ = {"-solve -bin -v 3 -v2 -pos", Name "GetDP/9ComputeCommand", Visible 0},
P_ = { "", Name "GetDP/2PostOperationChoices", Visible 0}
];