diff --git a/Elasticity/wrench2D.geo b/Elasticity/wrench2D.geo new file mode 100644 index 0000000000000000000000000000000000000000..782b7012ea859cd1cc7a8d35ff253340384d735e --- /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 0000000000000000000000000000000000000000..2291330453e0bf7a1ac06f2f2e789d64d7f94c83 --- /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 0000000000000000000000000000000000000000..578f64f094a83241d0b0e761e073ff1a49ef11ea --- /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); + + +