Commit 511d7f1b by Jeremy BLEYER

 .. Numerical tours of continuum mechanics using FEniCS documentation master file, created by sphinx-quickstart on Wed Jun 8 21:25:10 2016. You can adapt this file completely to your liking, but it should at least contain the root toctree directive. Linear problems in solid mechanics ================================== Contents: .. toctree:: :maxdepth: 1 demo/elasticity/2D_elasticity.py.rst demo/modal_analysis_dynamics/cantilever_modal.py.rst
 ... @@ -78,10 +78,9 @@ Injecting into :eq:constitutive_3D, we have for the 2D plane stress relation: ... @@ -78,10 +78,9 @@ Injecting into :eq:constitutive_3D, we have for the 2D plane stress relation: \boldsymbol{\sigma} = \lambda^* \text{tr}(\boldsymbol{\varepsilon})\mathbf{1} + 2\mu\boldsymbol{\varepsilon} \boldsymbol{\sigma} = \lambda^* \text{tr}(\boldsymbol{\varepsilon})\mathbf{1} + 2\mu\boldsymbol{\varepsilon} where :math:\boldsymbol{\sigma}, \boldsymbol{\varepsilon}, \mathbf{1} are 2D tensors and with where :math:\boldsymbol{\sigma}, \boldsymbol{\varepsilon}, \mathbf{1} are 2D tensors and with :math:\lambda^* = \dfrac{\lambda^2}{\lambda+\mu}. Hence, the 2D constitutive relation :math:\lambda^* = \dfrac{2\lambda\mu}{\lambda+2\mu}. Hence, the 2D constitutive relation is identical to the plane strain case by changing only the value of the Lamé coefficient :math:\lambda is identical to the plane strain case by changing only the value of the Lamé coefficient :math:\lambda. (equivalently, this corresponds to using a pseudo-Poisson coefficient :math:\nu^*=\dfrac{\nu}{1-\nu} We can then have:: instead of :math:\nu when defining :math:\lambda in :eq:Lame_coeff). We can then have:: E = Constant(1e5) E = Constant(1e5) nu = Constant(0.3) nu = Constant(0.3) ... @@ -90,7 +89,7 @@ instead of :math:\nu when defining :math:\lambda in :eq:Lame_coeff). We ca ... @@ -90,7 +89,7 @@ instead of :math:\nu when defining :math:\lambda in :eq:Lame_coeff). We ca mu = E/2/(1+nu) mu = E/2/(1+nu) lmbda = E*nu/(1+nu)/(1-2*nu) lmbda = E*nu/(1+nu)/(1-2*nu) if model == "plane_stress": if model == "plane_stress": lmbda = lmbda**2/(lmbda+mu) lmbda = 2*mu*lmbda/(lmbda+2*mu) def sigma(v): def sigma(v): return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v) return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v) ... @@ -160,5 +159,5 @@ Euler-Bernoulli beam theory which is here :math:w_{beam} = \dfrac{qL^4}{8EI}:: ... @@ -160,5 +159,5 @@ Euler-Bernoulli beam theory which is here :math:w_{beam} = \dfrac{qL^4}{8EI}:: print("Maximal deflection:", -u(L,H/2.)[1]) print("Maximal deflection:", -u(L,H/2.)[1]) print("Beam theory deflection:", float(3*rho_g*L**4/2/E/H**3)) print("Beam theory deflection:", float(3*rho_g*L**4/2/E/H**3)) One finds :math:w_{FE} = 5.8172\text{e-3} against :math:w_{beam} = 5.8594\text{e-3} One finds :math:w_{FE} = 5.8638\text{e-3} against :math:w_{beam} = 5.8594\text{e-3} that is a 0.72% difference. that is a 0.07% difference. \ No newline at end of file \ No newline at end of file
 ... @@ -109,7 +109,7 @@ Hermitian matrices. By default, SLEPc computes the largest eigenvalues, here ... @@ -109,7 +109,7 @@ Hermitian matrices. By default, SLEPc computes the largest eigenvalues, here we instead look for the smallest eigenvalues (they should all be real). To we instead look for the smallest eigenvalues (they should all be real). To improve convergence of the eigensolver for finding the smallest eigenvalues improve convergence of the eigensolver for finding the smallest eigenvalues (by default it computes the largest ones), a spectral transform is performed (by default it computes the largest ones), a spectral transform is performed using the keyword shift-invert i.e. the original problem is transformed into using the keyword shift-invert i.e. the original problem is transformed into an equivalent problem with eigenvalues given by :math:\dfrac{1}{\lambda - \sigma} an equivalent problem with eigenvalues given by :math:\dfrac{1}{\lambda - \sigma} instead of :math:\lambda where :math:\sigma is the value of the spectral shift. instead of :math:\lambda where :math:\sigma is the value of the spectral shift. It is therefore much easier to compute eigenvalues close to :math:\sigma i.e. It is therefore much easier to compute eigenvalues close to :math:\sigma i.e. ... @@ -118,7 +118,7 @@ transformed back by SLEPc to their original value :math:\lambda. ... @@ -118,7 +118,7 @@ transformed back by SLEPc to their original value :math:\lambda. We now ask SLEPc to extract the first 6 eigenvalues by calling its solve function We now ask SLEPc to extract the first 6 eigenvalues by calling its solve function and extract the corresponding eigenpair (first two arguments of get_eigenpair and extract the corresponding eigenpair (first two arguments of get_eigenpair correspond to the real and complex part of the eigenvalue, the last two to the correspond to the real and complex part of the eigenvalue, the last two to the real and complex part of the eigenvector):: real and complex part of the eigenvector):: ... @@ -168,8 +168,8 @@ where :math:S=BH is the beam section, :math:I the bending inertia and ... @@ -168,8 +168,8 @@ where :math:S=BH is the beam section, :math:I the bending inertia and the solution of which can be well approximated by :math:(2n+1)\pi/2 for :math:n\geq 3. the solution of which can be well approximated by :math:(2n+1)\pi/2 for :math:n\geq 3. Since the beam possesses two bending axis, each solution to the previous equation is Since the beam possesses two bending axis, each solution to the previous equation is associated to two frequencies with bending along the weak axis (:math:I=I_{\text{weak}} = HB^3/12) associated with two frequencies, one with bending along the weak axis (:math:I=I_{\text{weak}} = HB^3/12) and along the strong axis (:math:I=I_{\text{strong}} = BH^3/12). Since :math:I_{\text{strong}} = 4I_{\text{weak}} and the other along the strong axis (:math:I=I_{\text{strong}} = BH^3/12). Since :math:I_{\text{strong}} = 4I_{\text{weak}} for the considered numerical values, the strong axis bending frequency will be twice that corresponsing for the considered numerical values, the strong axis bending frequency will be twice that corresponsing to bending along the weak axis. The solution :math:\alpha_n are computed using the to bending along the weak axis. The solution :math:\alpha_n are computed using the scipy.optimize.root function with initial guess given by :math:(2n+1)\pi/2. scipy.optimize.root function with initial guess given by :math:(2n+1)\pi/2. ... ...
 ... @@ -41,7 +41,7 @@ Plate bending stiffness :math:D=\dfrac{Eh^3}{12(1-\nu^2)} and shear stiffness ... @@ -41,7 +41,7 @@ Plate bending stiffness :math:D=\dfrac{Eh^3}{12(1-\nu^2)} and shear stiffness with a shear correction factor :math:\kappa = 5/6 for a homogeneous plate with a shear correction factor :math:\kappa = 5/6 for a homogeneous plate of thickness :math:h:: of thickness :math:h:: thick = Constant(1e-2) thick = Constant(1e-3) D = E*thick**3/(1-nu**2)/12. D = E*thick**3/(1-nu**2)/12. F = E/2/(1+nu)*thick*5./6. F = E/2/(1+nu)*thick*5./6. ... @@ -52,7 +52,7 @@ constant value in the thin plate Love-Kirchhoff limit:: ... @@ -52,7 +52,7 @@ constant value in the thin plate Love-Kirchhoff limit:: The unit square mesh is divided in :math:N\times N quadrilaterals:: The unit square mesh is divided in :math:N\times N quadrilaterals:: N = 40 N = 50 mesh = UnitSquareMesh.create(N, N, CellType.Type_quadrilateral) mesh = UnitSquareMesh.create(N, N, CellType.Type_quadrilateral) Continuous interpolation using of degree :math:d=\texttt{deg} is chosen for both deflection and rotation:: Continuous interpolation using of degree :math:d=\texttt{deg} is chosen for both deflection and rotation:: ... @@ -127,5 +127,9 @@ We then solve for the solution and export the relevant fields to XDMF files :: ... @@ -127,5 +127,9 @@ We then solve for the solution and export the relevant fields to XDMF files :: The solution is compared to the Kirchhoff analytical solution:: The solution is compared to the Kirchhoff analytical solution:: print "Kirchhoff deflection:", -1.265319087e-3*float(f/D) print "Kirchhoff deflection:", -1.265319087e-3*float(f/D) print "Reissner-Mindlin FE deflection:", min(ww.vector().get_local()) # point evaluation for quads print "Reissner-Mindlin FE deflection:", -min(ww.vector().get_local()) # point evaluation for quads # is not implemented in fenics 2017.2 # is not implemented in fenics 2017.2 \ No newline at end of file For N=50 quads per side, one finds :math:w_{FE} = 1.38182\text{e-5} for linear quads and :math:w_{FE} = 1.38176\text{e-5} for quadratic quads against :math:w_{\text{Kirchhoff}} = 1.38173\text{e-5} for the thin plate solution. \ No newline at end of file