From e4ce8affd03b2b00160a569dffec3ce3f4e9b1cc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Henrotte?= <francois.henrotte@uclouvain.be>
Date: Tue, 3 Oct 2017 09:17:50 +0200
Subject: [PATCH] tutorial 3: linear elasticity
---
Elasticity/wrench2D.geo | 72 +++++++
Elasticity/wrench2D.pro | 366 +++++++++++++++++++++++++++++++++
Elasticity/wrench2D_common.pro | 29 +++
3 files changed, 467 insertions(+)
create mode 100644 Elasticity/wrench2D.geo
create mode 100644 Elasticity/wrench2D.pro
create mode 100644 Elasticity/wrench2D_common.pro
diff --git a/Elasticity/wrench2D.geo b/Elasticity/wrench2D.geo
new file mode 100644
index 0000000..782b701
--- /dev/null
+++ b/Elasticity/wrench2D.geo
@@ -0,0 +1,72 @@
+
+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};
+
+
+
+
diff --git a/Elasticity/wrench2D.pro b/Elasticity/wrench2D.pro
new file mode 100644
index 0000000..2291330
--- /dev/null
+++ b/Elasticity/wrench2D.pro
@@ -0,0 +1,366 @@
+/* -------------------------------------------------------------------
+ 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}
+];
+
diff --git a/Elasticity/wrench2D_common.pro b/Elasticity/wrench2D_common.pro
new file mode 100644
index 0000000..578f64f
--- /dev/null
+++ b/Elasticity/wrench2D_common.pro
@@ -0,0 +1,29 @@
+// 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);
+
+
+
--
GitLab