tutorial 3: linear elasticity

Elasticity/wrench2D.geo

+
+Include "wrench2D_common.pro";
+
+Solver.AutoMesh = 2;
+
+lc = Refine; 
+lc2 = lc; // nut region
+lc3 = lc/5; // edge grip region
+
+Ri = 0.494*in;
+Re = 0.8*in;
+L = LLength;
+LL = 1*in;
+E = 0.625*in;
+F = Width;
+D = 0.45*in; // Nut
+
+theta = Inclination;
+
+Point(1) = { 0, 0, 0, lc2};
+Point(2) = { -Ri, 0, 0, lc2};
+th1 = Asin[E/2./Ri];
+Point(3) = { -Ri + Ri*Cos[th1], Ri*Sin[th1], 0, lc2};
+th2 = Asin[E/2./Re];
+Point(4) = { -Re*Cos[th2]+D, Re*Sin[th2], 0, lc3};
+Point(5) = { -Re*Cos[th2], Re*Sin[th2], 0, lc2};
+th3 = th2 - theta;
+Point(6) = { Re*Cos[th3], Re*Sin[th3], 0, lc};
+Point(7) = {  (L-LL)*Cos[theta]+F/2.*Sin[theta],
+	     -(L-LL)*Sin[theta]+F/2.*Cos[theta], 0, lc};
+Point(8) = { L*Cos[theta]+F/2.*Sin[theta], 
+	     -L*Sin[theta]+F/2.*Cos[theta], 0, lc};
+Point(9) = { L*Cos[theta]-F/2.*Sin[theta], 
+	     -L*Sin[theta]-F/2.*Cos[theta], 0, lc};
+th4 = -th2 - theta;
+Point(10) = { Re*Cos[th4], Re*Sin[th4], 0, lc};
+Point(11) = { -Re*Cos[th2], -Re*Sin[th2], 0, lc2};
+Point(12) = { -Re*Cos[th2]+D, -Re*Sin[th2], 0, lc3};
+Point(13) = { -Ri + Ri*Cos[th1], -Ri*Sin[th1], 0, lc2};
+
+contour[] = {};
+contour[] += newl; Circle(newl) = {1,2,3};
+contour[] += newl; Line(newl) = {3,4};
+contour[] += newl; cl1=newl; Line(newl) = {4,5}; 
+contour[] += newl; Circle(newl) = {5,1,6};
+contour[] += newl; Line(newl) = {6,7};
+contour[] += newl; cl2=newl; Line(newl) = {7,8}; 
+contour[] += newl; Line(newl) = {8,9};
+contour[] += newl; Line(newl) = {9,10};
+contour[] += newl; Circle(newl) = {10,1,11};
+contour[] += newl; cl3=newl; Line(newl) = {11,12}; 
+contour[] += newl; Line(newl) = {12,13};
+contour[] += newl; Circle(newl) = {13,2,1};
+ll=newll; Line Loop(ll) = { contour[] };
+
+wrench=news; Plane Surface(wrench) = {-ll}; 
+
+/* Using quadrangular elements instead of triangular elements is pretty simple.
+   You just have to invoke the Recombine command at this place.
+   GetDP deals will all the rest by itself. 
+ */
+If(Recomb==1)
+  Recombine Surface{wrench};
+EndIf
+
+Physical Surface(1)= wrench;
+Physical Line(2)= {cl1, cl3};
+Physical Line(3)= {cl2};
+
+
+
+

Elasticity/wrench2D.pro

+/* -------------------------------------------------------------------
+   Tutorial 3 : linear elastic model of a wrench
+
+   Features:
+   - "Grad u" GetDP specific formulation for linear elasticity
+   - first and second order elements 
+   - triangular and quadrangular elements
+
+   To compute the solution interactively from the Gmsh GUI:
+       File > Open > wrench.pro
+       Run (button at the bottom of the left panel)
+   ------------------------------------------------------------------- */
+
+/* Linear elasticity with GetDP:
+
+   GetDP has a peculiar way to deal with linear elasticity.
+   Instead of vector field "u = Vector[ ux, uy, uz ]",
+   the displacement field is regarded as two (2D case) or 3 (3D case) scalar fields.
+   Unlike conventional formulations, GetDP's formulation is written 
+   in terms of the gradient "Grad u" of the displacement field,
+   which is a non-symmetric tensor, and the needed symmetrization 
+   (to define the strain tensor and relate it to the stress tensor) 
+   is done through the constitutive relationship (Hooke law). 
+   The reason for this unusual formulation is to be able to use also for elastic problems
+   the powerful geometrical and homological kernel of GetDP,
+   which relies on the operators Grad, Curl and Div.
+
+   The "Grad u" formulation entails a small increase of assembly work
+   but makes in counterpart lots of geometrical features
+   implemented in GetDP (change of coordinates, function spaces, etc...)
+   applicable to elastic problems out-of-the-box, 
+   since the scalar fields { ux, uy, uz } have exactly the same geometrical properties
+   as, e.g. a scalar electrip potential or a temperature field. 
+*/
+
+
+Include "wrench2D_common.pro";
+
+Young = 1e9 * 
+  DefineNumber[ 200, Name "Material/Young modulus [GPa]"];
+Poisson = 
+  DefineNumber[ 0.3, Name "Material/Poisson coefficient []"];
+AppliedForce = 
+  DefineNumber[ 100, Name "Material/Applied force [N]"];
+
+// Approximation of the maximum deflection by an analytical model: 
+// Deflection = PL^3/(3EI) with I = Width^3*Thickness/12
+Deflection = 
+  DefineNumber[4*AppliedForce*((LLength-0.018)/Width)^3/(Young*Thickness)*1e3, 
+	       Name "Solution/Deflection (analytical) [mm]", ReadOnly 1];
+
+Group {
+  // Physical regions: Give explicit labels to the regions defined in the .geo file 
+  Wrench = Region[ 1 ];
+  Grip = Region[ 2 ];
+  Force = Region[ 3 ];
+
+  // Abstract regions:
+  Vol_Elast_Mec = Region[ { Wrench } ];
+  Vol_Force_Mec = Region[ { Wrench } ];
+  Sur_Clamp_Mec = Region[ { Grip } ];
+  Sur_Force_Mec = Region[ { Force } ];
+/* 
+  Signification of the abstract regions:
+  Vol_Elast_Mec       Elastic domain
+  Vol_Force_Mec       Region with imposed volumic force
+  Sur_Force_Mec       Surface with imposed surface traction
+  Sur_Clamp_Mec       Surface with imposed zero displacements (all components)
+*/
+}
+
+Function {
+  /* Material coefficients. 
+     No need to define them regionwise here ( E[{Wrench}] = ... ; )
+     as there is only one region in this model. */
+  E[] = Young;
+  nu[] = Poisson;
+  /* Components of the volumic force applied to the region "Vol_Force_Mec"
+     Gravity could be defined here ( force_y[] = 7000*9.81; ) ; */
+  force_x[] = 0;
+  force_y[] = 0;
+  /* Components of the surface traction force applied to the region "Sur_Force_Mec" */
+  pressure_x[] = 0;
+  pressure_y[] = AppliedForce/(SurfaceArea[]*Thickness);
+}
+
+
+/* Hooke law
+
+   The material law 
+
+   sigma_ij = C_ijkl epsilon_ij 
+
+   is represented in 2D by 4 2x2 tensors C_ij[], i,j=1,2 
+   depending on the Lamé coefficients of the isotropic linear material,
+
+   lambda = E[]*nu[]/(1.+nu[])/(1.-2.*nu[]);
+   mu = E[]/2./(1.+nu[]);
+
+   as follows
+
+   EPC:  a[] = E/(1-nu^2)        b[] = mu     c[] = E nu/(1-nu^2)
+   EPD:  a[] = lambda + 2 mu     b[] = mu     c[] = lambda
+    3D:  a[] = lambda + 2 mu     b[] = mu     c[] = lambda
+
+   respectively for the 2D plane strain (EPD), 2D plane stress (EPS) and 3D cases.
+*/
+
+Function {
+  Flag_EPC = 1; 
+  If(Flag_EPC) // Plane stress
+    a[] = E[]/(1.-nu[]^2);
+    c[] = E[]*nu[]/(1.-nu[]^2);
+  Else // Plane strain or 3D
+    a[] = E[]*(1.-nu[])/(1.+nu[])/(1.-2.*nu[]);
+    c[] = E[]*nu[]/(1.+nu[])/(1.-2.*nu[]);
+  EndIf
+  b[] = E[]/2./(1.+nu[]);
+
+  C_xx[] = Tensor[ a[],0  ,0  ,    0  ,b[],0  ,    0  ,0  ,b[] ];
+  C_xy[] = Tensor[ 0  ,c[],0  ,    b[],0  ,0  ,    0  ,0  ,0   ];
+
+  C_yx[] = Tensor[ 0  ,b[],0  ,    c[],0  ,0  ,    0  ,0  ,0   ];
+  C_yy[] = Tensor[ b[],0  ,0  ,    0  ,a[],0  ,    0  ,0  ,b[] ];
+}
+
+/* Clamping boundary condition */
+Constraint {
+  { Name Displacement_x;
+    Case {
+      { Region Sur_Clamp_Mec ; Type Assign ; Value 0; }
+    }
+  }
+  { Name Displacement_y;
+    Case {
+       { Region Sur_Clamp_Mec ; Type Assign ; Value 0; }
+    }
+  }
+}
+
+/* As explained above, the displacement field is discretized 
+   as two scalar fields "ux" and "uy", which are the spatial components
+   of the vector field "u" in a fixed Cartesian coordinate system. 
+
+   Boundary conditions like
+
+   ux = ... ;
+   uy = ... ;
+
+   translate naturally into Dirichlet constraints 
+   on the scalar "ux" and "uy" FunctionSpaces. 
+   Conditions like, however,
+
+   u . n = ux Cos [th] + uy Sin [th] = ... ;
+
+   are less naturally accounted for within the "Grad u" formulation.
+   Fortunately, they are rather uncommon. 
+
+   Finite element shape (triangles or quadrangles) makes no difference
+   in the definition of the FunctionSpaces. The appropriate shape functions to be used
+   are determined by GetDP at a much lower level on basis of the information
+   contained in the *.msh file. 
+
+   Second order elements, on the other hand, are implemented in the hierarchical fashion
+   by adding to the first order node-based shape functions 
+   a set of second order edge-based functions
+   to complete a basis for 2d order polynomials on the reference element. 
+ */
+
+
+// Domain of definition of the "ux" and "uy" FunctionSpaces
+Group {
+  Dom_H_u_Mec = Region[ { Vol_Elast_Mec, 
+			  Sur_Force_Mec, 
+			  Sur_Clamp_Mec} ];
+}
+
+Flag_Degree = 
+  DefineNumber[ 0, Name "Geometry/Use degree 2 (hierarch.)", Choices{0,1}, Visible 1];
+FE_Degree = ( Flag_Degree == 0 ) ? 1 : 2; // Convert flag value into polynomial degree
+
+FunctionSpace {
+  { Name H_ux_Mec ; Type Form0 ;
+    BasisFunction {
+      { Name sxn ; NameOfCoef uxn ; Function BF_Node ;
+        Support Dom_H_u_Mec ; Entity NodesOf[ All ] ; }
+     If ( FE_Degree == 2 )
+        { Name sxn2 ; NameOfCoef uxn2 ; Function BF_Node_2E ;
+          Support Dom_H_u_Mec; Entity EdgesOf[ All ] ; }
+     EndIf
+    }
+    Constraint {
+      { NameOfCoef uxn ;
+        EntityType NodesOf ; NameOfConstraint Displacement_x ; }
+      If ( FE_Degree == 2 )
+         { NameOfCoef uxn2 ;
+	   EntityType EdgesOf ; NameOfConstraint Displacement_x ; }
+     EndIf
+    }
+  }
+  { Name H_uy_Mec ; Type Form0 ;
+    BasisFunction {
+      { Name syn ; NameOfCoef uyn ; Function BF_Node ;
+        Support Dom_H_u_Mec ; Entity NodesOf[ All ] ; }
+     If ( FE_Degree == 2 )
+        { Name syn2 ; NameOfCoef uyn2 ; Function BF_Node_2E ;
+          Support Dom_H_u_Mec; Entity EdgesOf[ All ] ; }
+     EndIf
+    }
+    Constraint {
+      { NameOfCoef uyn ;
+        EntityType NodesOf ; NameOfConstraint Displacement_y ; }
+      If ( FE_Degree == 2 )
+      { NameOfCoef uyn2 ;
+        EntityType EdgesOf ; NameOfConstraint Displacement_y ; }
+      EndIf
+    }
+  }
+}
+
+
+Jacobian {
+  { Name Vol;
+    Case {
+      { Region All; Jacobian Vol; }
+    }
+  }
+  { Name Sur;
+    Case {
+      { Region All; Jacobian Sur; }
+    }
+  }
+}
+
+/* Adapt the number of Gauss points to the polynomial degree of the finite elements
+is as simple as this: */
+If (FE_Degree == 1)
+Integration {
+  { Name Gauss_v;
+    Case {
+      {Type Gauss;
+        Case {
+	  { GeoElement Line       ; NumberOfPoints  3; }
+          { GeoElement Triangle   ; NumberOfPoints  3; }
+          { GeoElement Quadrangle ; NumberOfPoints  4; }
+        }
+      }
+    }
+  }
+}
+ElseIf (FE_Degree == 2)
+Integration {
+  { Name Gauss_v;
+    Case {
+      {Type Gauss;
+        Case {
+	  { GeoElement Line       ; NumberOfPoints  5; }
+          { GeoElement Triangle   ; NumberOfPoints  7; }
+          { GeoElement Quadrangle ; NumberOfPoints  7; }
+        }
+      }
+    }
+  }
+}
+EndIf
+
+Formulation {
+  { Name Elast_u ; Type FemEquation ;
+    Quantity {
+      { Name ux  ; Type Local ; NameOfSpace H_ux_Mec ; }
+      { Name uy  ; Type Local ; NameOfSpace H_uy_Mec ; }
+    }
+    Equation {
+      Galerkin { [ -C_xx[] * Dof{d ux}, {d ux} ] ;
+        In Vol_Elast_Mec ; Jacobian Vol ; Integration Gauss_v ; }
+      Galerkin { [ -C_xy[] * Dof{d uy}, {d ux} ] ;
+        In Vol_Elast_Mec ; Jacobian Vol ; Integration Gauss_v ; }
+      Galerkin { [ -C_yx[] * Dof{d ux}, {d uy} ] ;
+        In Vol_Elast_Mec ; Jacobian Vol ; Integration Gauss_v ; }
+      Galerkin { [ -C_yy[] * Dof{d uy}, {d uy} ] ;
+        In Vol_Elast_Mec ; Jacobian Vol ; Integration Gauss_v ; }
+
+      Galerkin { [ force_x[] , {ux} ];
+        In Vol_Force_Mec ; Jacobian Vol ; Integration Gauss_v ; }
+      Galerkin { [ force_y[] , {uy} ];
+        In Vol_Force_Mec ; Jacobian Vol ; Integration Gauss_v ; }
+      
+      Galerkin { [ pressure_x[] , {ux} ];
+        In Sur_Force_Mec ; Jacobian Sur ; Integration Gauss_v ; }
+      Galerkin { [ pressure_y[] , {uy} ];
+        In Sur_Force_Mec ; Jacobian Sur ; Integration Gauss_v ; }
+    }
+  }
+}
+
+Resolution {
+  { Name Elast_u ;
+    System {
+      { Name Sys_Mec ; NameOfFormulation Elast_u; }
+    }
+    Operation {
+      InitSolution [Sys_Mec];
+      Generate[Sys_Mec];
+      Solve[Sys_Mec];
+      SaveSolution[Sys_Mec] ;
+    }
+  }
+}
+
+PostProcessing {
+  { Name Elast_u ; NameOfFormulation Elast_u ;
+    PostQuantity {
+      { Name u ; Value { Term { [ Vector[ {ux}, {uy}, 0 ]]; 
+	      In Vol_Elast_Mec ; Jacobian Vol ; } } }
+      { Name uy ; Value { Term { [ 1e3*{uy} ]; 
+	      In Vol_Elast_Mec ; Jacobian Vol ; } } }
+      { Name sig_xx ; Value { Term { 
+	    [ CompX[  C_xx[]*{d ux} + C_xy[]*{d uy} ] ]; 
+	    In Vol_Elast_Mec ; Jacobian Vol ; } } }
+      { Name sig_xy ; Value { Term { 
+	    [ CompY[  C_xx[]*{d ux} + C_xy[]*{d uy} ] ]; 
+	    In Vol_Elast_Mec ; Jacobian Vol ; } } }
+      { Name sig_yy ; Value { Term { 
+	    [  CompY [ C_yx[]*{d ux} + C_yy[]*{d uy} ] ]; 
+	    In Vol_Elast_Mec ; Jacobian Vol ; } } }
+    }
+  }
+}
+
+
+
+PostOperation {
+  { Name pos; NameOfPostProcessing Elast_u;
+    Operation {
+
+      If(FE_Degree == 1) 
+	Print[ sig_xx, OnElementsOf Wrench, File "sigxx.pos" ];
+        Print[ u, OnElementsOf Wrench, File "u.pos" ];
+      Else
+	Print[ sig_xx, OnElementsOf Wrench, File "sigxx2.pos" ];
+        Print[ u, OnElementsOf Wrench, File "u2.pos" ];
+      EndIf
+      Echo[ StrCat["l=PostProcessing.NbViews-1; ", 
+      		   "View[l].VectorType = 5; ",
+      		   "View[l].ExternalView = l; ",
+      		   "View[l].DisplacementFactor = 100; ",
+      		   "View[l-1].IntervalsType = 3; "
+		   ],
+             File "tmp.geo", LastTimeStepOnly] ;
+      //Print[ sig_yy, OnElementsOf Wrench, File "sigyy.pos" ];
+      //Print[ sig_xy, OnElementsOf Wrench, File "sigxy.pos" ];
+      Print[ uy, OnPoint{probe_x, probe_y, 0}, 
+	     File > "deflection.pos", Format TimeTable,
+	     SendToServer "Solution/Deflection (computed) [mm]", Color "AliceBlue" ];
+    }
+  }
+}
+
+
+// Tell Gmsh which GetDP commands to execute when running the model.
+DefineConstant[
+  R_ = {"Elast_u", Name "GetDP/1ResolutionChoices", Visible 0},
+  P_ = {"pos", Name "GetDP/2PostOperationChoices", Visible 0},
+  C_ = {"-solve -pos -v2", Name "GetDP/9ComputeCommand", Visible 0}
+];
+

Elasticity/wrench2D_common.pro

+// Some useful conversion coefficients
+mm = 1.e-3;     // millimeters to meters
+cm = 1.e-2;     // centimeters to meters
+in = 0.0254;    // inches to meters
+deg = Pi/180.;  // degrees to radians
+
+Refine = 
+  mm*DefineNumber[ 2, Name "Geometry/2Main characteristic length"];
+Recomb = 
+  DefineNumber[ 0, Name "Geometry/1Recombine", Choices{0,1}];
+Thickness = 
+  mm*DefineNumber[ 10, Name "Geometry/4Thickness (mm)"];
+Width = 
+  mm*DefineNumber[ 0.625*in/mm, Name "Geometry/5Arm Width (mm)"];
+LLength = 
+  cm*DefineNumber[ 6.0*in/cm, Name "Geometry/6Arm Length (cm)"];
+
+// Definition of the coordinates of the arm end
+// at which the maximal deflection will be evaluated.
+Inclination = 14*deg; // Arm angle with x-axis
+eps = 1e-6;
+probe_x = (LLength-eps)*Cos[Inclination];
+probe_y =-(LLength-eps)*Sin[Inclination];
+
+// a useful Print command to debug a model or communicate with the user
+Printf("Maximal deflection calculated at point (%f,%f)", probe_x, probe_y);
+       
+
+