Commit e9ee2373 by Jeremy BLEYER

parent 53f479f5
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
 # # # .. # gedit: set fileencoding=utf8 : # .. # gedit: set fileencoding=utf8 : # .. raw:: html # #

# # # .. _LinearElasticity2D: # # # .. _LinearElasticity2D: # # # ========================= # ========================= # 2D linear elasticity # 2D linear elasticity ... ...
 # raw:: html

# .. _ModalAnalysis: # .. _ModalAnalysis: # # # ========================================== # ========================================== ... ...
 # # # .. # gedit: set fileencoding=utf8 : # .. # gedit: set fileencoding=utf8 : # .. raw:: html # #

# # # .. _OrthotropicElasticity: # .. _OrthotropicElasticity: # # ... @@ -11,10 +14,10 @@ ... @@ -11,10 +14,10 @@ # Introduction # Introduction # ------------ # ------------ # # # In this numerical tour, we will show how to tackle the case of orthotropic elasticity (in a 2D setting). The corresponding file can be obtained from # In this numerical tour, we will show how to tackle the case of orthotropic elasticity (in a 2D setting). The corresponding file can be obtained from # :download:orthotropic_elasticity.py. # :download:orthotropic_elasticity.py. # # # We consider here the case of a square plate perforated by a circular hole of # We consider here the case of a square plate perforated by a circular hole of # radius :math:R, the plate dimension is :math:2L\times 2L with :math:L \gg R # radius :math:R, the plate dimension is :math:2L\times 2L with :math:L \gg R # Only the top-right quarter of the plate will be considered. Loading will consist # Only the top-right quarter of the plate will be considered. Loading will consist # of a uniform traction on the top/bottom boundaries, symmetry conditions will also # of a uniform traction on the top/bottom boundaries, symmetry conditions will also ... @@ -27,7 +30,7 @@ from mshr import * ... @@ -27,7 +30,7 @@ from mshr import * L, R = 1., 0.1 L, R = 1., 0.1 N = 50 # mesh density N = 50 # mesh density domain = Rectangle(Point(0.,0.), Point(L, L)) - Circle(Point(0., 0.), R) domain = Rectangle(Point(0.,0.), Point(L, L)) - Circle(Point(0., 0.), R) mesh = generate_mesh(domain, N) mesh = generate_mesh(domain, N) ... @@ -35,7 +38,7 @@ mesh = generate_mesh(domain, N) ... @@ -35,7 +38,7 @@ mesh = generate_mesh(domain, N) # Constitutive relation # Constitutive relation # --------------------- # --------------------- # # # Constitutive relations will be defined using an engineering (or Voigt) notation (i.e. # Constitutive relations will be defined using an engineering (or Voigt) notation (i.e. # second order tensors will be written as a vector of their components) contrary # second order tensors will be written as a vector of their components) contrary # to the :ref:LinearElasticity2D example which used an intrinsic notation. In # to the :ref:LinearElasticity2D example which used an intrinsic notation. In # the material frame, which is assumed to coincide here with the global :math:(Oxy) # the material frame, which is assumed to coincide here with the global :math:(Oxy) ... @@ -46,13 +49,13 @@ mesh = generate_mesh(domain, N) ... @@ -46,13 +49,13 @@ mesh = generate_mesh(domain, N) # .. math:: # .. math:: # \begin{Bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ 2\varepsilon_{xy} # \begin{Bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ 2\varepsilon_{xy} # \end{Bmatrix} = \begin{bmatrix} 1/E_x & -\nu_{xy}/E_x & 0\\ # \end{Bmatrix} = \begin{bmatrix} 1/E_x & -\nu_{xy}/E_x & 0\\ # -\nu_{yx}/E_y & 1/E_y & 0 \\ 0 & 0 & 1/G_{xy} \end{bmatrix}\begin{Bmatrix} # -\nu_{yx}/E_y & 1/E_y & 0 \\ 0 & 0 & 1/G_{xy} \end{bmatrix}\begin{Bmatrix} # \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} # \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} # \end{Bmatrix} # \end{Bmatrix} # # # with :math:E_x, E_y the two Young's moduli in the orthotropy directions, :math:\nu_{xy} # with :math:E_x, E_y the two Young's moduli in the orthotropy directions, :math:\nu_{xy} # the in-plane Poisson ration (with the following relation ensuring the constitutive # the in-plane Poisson ration (with the following relation ensuring the constitutive # relation symmetry :math:\nu_{yx}=\nu_{xy}E_y/E_x) and :math:G_{xy} being the # relation symmetry :math:\nu_{yx}=\nu_{xy}E_y/E_x) and :math:G_{xy} being the # shear modulus. This relation needs to be inverted to obtain the stress components as a function # shear modulus. This relation needs to be inverted to obtain the stress components as a function # of the strain components :math:\boldsymbol{\sigma}=\mathbf{C}\boldsymbol{\varepsilon} with # of the strain components :math:\boldsymbol{\sigma}=\mathbf{C}\boldsymbol{\varepsilon} with # :math:\mathbf{C}=\mathbf{S}^{-1}:: # :math:\mathbf{C}=\mathbf{S}^{-1}:: ... @@ -89,7 +92,7 @@ def sigma(v): ... @@ -89,7 +92,7 @@ def sigma(v): # Problem position and resolution # Problem position and resolution # -------------------------------- # -------------------------------- # # # Different parts of the quarter plate boundaries are now defined as well as the # Different parts of the quarter plate boundaries are now defined as well as the # exterior integration measure ds:: # exterior integration measure ds:: class Top(SubDomain): class Top(SubDomain): ... @@ -127,7 +130,7 @@ a = inner(sigma(du), eps(u_))*dx ... @@ -127,7 +130,7 @@ a = inner(sigma(du), eps(u_))*dx T = Constant((0, 1e-3)) T = Constant((0, 1e-3)) l = dot(T, u_)*ds(1) l = dot(T, u_)*ds(1) # Symmetric boundary conditions are applied on the Top and Left boundaries # Symmetric boundary conditions are applied on the Top and Left boundaries # and the problem is solved:: # and the problem is solved:: # symmetry boundary conditions # symmetry boundary conditions ... ...
 # raw:: html

# .. _ReissnerMindlinDG: # .. _ReissnerMindlinDG: # # # ============================================================== # ============================================================== ... ...
 # raw:: html

# .. _ReissnerMindlinQuads: # .. _ReissnerMindlinQuads: # # # ========================================== # ========================================== ... ...
 .. _vonMisesPlasticity: .. _vonMisesPlasticity: .. raw:: html

================================================== ================================================== Elasto-plastic analysis of a 2D von Mises material Elasto-plastic analysis of a 2D von Mises material ================================================== ================================================== ... ...
 .. # gedit: set fileencoding=utf8 : .. # gedit: set fileencoding=utf8 : .. raw:: html

.. _LinearElasticity2D: .. _LinearElasticity2D: ========================= ========================= 2D linear elasticity 2D linear elasticity ... ...
 .. # gedit: set fileencoding=utf8 : .. # gedit: set fileencoding=utf8 : .. raw:: html

.. _OrthotropicElasticity: .. _OrthotropicElasticity: ... @@ -11,10 +14,10 @@ ... @@ -11,10 +14,10 @@ Introduction Introduction ------------ ------------ In this numerical tour, we will show how to tackle the case of orthotropic elasticity (in a 2D setting). The corresponding file can be obtained from In this numerical tour, we will show how to tackle the case of orthotropic elasticity (in a 2D setting). The corresponding file can be obtained from :download:orthotropic_elasticity.py. :download:orthotropic_elasticity.py. We consider here the case of a square plate perforated by a circular hole of We consider here the case of a square plate perforated by a circular hole of radius :math:R, the plate dimension is :math:2L\times 2L with :math:L \gg R radius :math:R, the plate dimension is :math:2L\times 2L with :math:L \gg R Only the top-right quarter of the plate will be considered. Loading will consist Only the top-right quarter of the plate will be considered. Loading will consist of a uniform traction on the top/bottom boundaries, symmetry conditions will also of a uniform traction on the top/bottom boundaries, symmetry conditions will also ... @@ -27,7 +30,7 @@ between a rectangle and a circle:: ... @@ -27,7 +30,7 @@ between a rectangle and a circle:: L, R = 1., 0.1 L, R = 1., 0.1 N = 50 # mesh density N = 50 # mesh density domain = Rectangle(Point(0.,0.), Point(L, L)) - Circle(Point(0., 0.), R) domain = Rectangle(Point(0.,0.), Point(L, L)) - Circle(Point(0., 0.), R) mesh = generate_mesh(domain, N) mesh = generate_mesh(domain, N) ... @@ -35,7 +38,7 @@ between a rectangle and a circle:: ... @@ -35,7 +38,7 @@ between a rectangle and a circle:: Constitutive relation Constitutive relation --------------------- --------------------- Constitutive relations will be defined using an engineering (or Voigt) notation (i.e. Constitutive relations will be defined using an engineering (or Voigt) notation (i.e. second order tensors will be written as a vector of their components) contrary second order tensors will be written as a vector of their components) contrary to the :ref:LinearElasticity2D example which used an intrinsic notation. In to the :ref:LinearElasticity2D example which used an intrinsic notation. In the material frame, which is assumed to coincide here with the global :math:(Oxy) the material frame, which is assumed to coincide here with the global :math:(Oxy) ... @@ -46,13 +49,13 @@ frame, the orthotropic constitutive law writes :math:\boldsymbol{\varepsilon}=\ ... @@ -46,13 +49,13 @@ frame, the orthotropic constitutive law writes :math:\boldsymbol{\varepsilon}=\ .. math:: .. math:: \begin{Bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ 2\varepsilon_{xy} \begin{Bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ 2\varepsilon_{xy} \end{Bmatrix} = \begin{bmatrix} 1/E_x & -\nu_{xy}/E_x & 0\\ \end{Bmatrix} = \begin{bmatrix} 1/E_x & -\nu_{xy}/E_x & 0\\ -\nu_{yx}/E_y & 1/E_y & 0 \\ 0 & 0 & 1/G_{xy} \end{bmatrix}\begin{Bmatrix} -\nu_{yx}/E_y & 1/E_y & 0 \\ 0 & 0 & 1/G_{xy} \end{bmatrix}\begin{Bmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} \end{Bmatrix} \end{Bmatrix} with :math:E_x, E_y the two Young's moduli in the orthotropy directions, :math:\nu_{xy} with :math:E_x, E_y the two Young's moduli in the orthotropy directions, :math:\nu_{xy} the in-plane Poisson ration (with the following relation ensuring the constitutive the in-plane Poisson ration (with the following relation ensuring the constitutive relation symmetry :math:\nu_{yx}=\nu_{xy}E_y/E_x) and :math:G_{xy} being the relation symmetry :math:\nu_{yx}=\nu_{xy}E_y/E_x) and :math:G_{xy} being the shear modulus. This relation needs to be inverted to obtain the stress components as a function shear modulus. This relation needs to be inverted to obtain the stress components as a function of the strain components :math:\boldsymbol{\sigma}=\mathbf{C}\boldsymbol{\varepsilon} with of the strain components :math:\boldsymbol{\sigma}=\mathbf{C}\boldsymbol{\varepsilon} with :math:\mathbf{C}=\mathbf{S}^{-1}:: :math:\mathbf{C}=\mathbf{S}^{-1}:: ... @@ -89,7 +92,7 @@ second-order tensor or using the Voigt engineering notation:: ... @@ -89,7 +92,7 @@ second-order tensor or using the Voigt engineering notation:: Problem position and resolution Problem position and resolution -------------------------------- -------------------------------- Different parts of the quarter plate boundaries are now defined as well as the Different parts of the quarter plate boundaries are now defined as well as the exterior integration measure ds:: exterior integration measure ds:: class Top(SubDomain): class Top(SubDomain): ... @@ -101,7 +104,7 @@ exterior integration measure ds:: ... @@ -101,7 +104,7 @@ exterior integration measure ds:: class Bottom(SubDomain): class Bottom(SubDomain): def inside(self, x, on_boundary): def inside(self, x, on_boundary): return near(x[1],0) and on_boundary return near(x[1],0) and on_boundary # exterior facets MeshFunction # exterior facets MeshFunction facets = MeshFunction("size_t", mesh, 1) facets = MeshFunction("size_t", mesh, 1) facets.set_all(0) facets.set_all(0) ... @@ -122,12 +125,12 @@ on the top boundary:: ... @@ -122,12 +125,12 @@ on the top boundary:: u_ = TestFunction(V) u_ = TestFunction(V) u = Function(V, name='Displacement') u = Function(V, name='Displacement') a = inner(sigma(du), eps(u_))*dx a = inner(sigma(du), eps(u_))*dx # uniform traction on top boundary # uniform traction on top boundary T = Constant((0, 1e-3)) T = Constant((0, 1e-3)) l = dot(T, u_)*ds(1) l = dot(T, u_)*ds(1) Symmetric boundary conditions are applied on the Top and Left boundaries Symmetric boundary conditions are applied on the Top and Left boundaries and the problem is solved:: and the problem is solved:: # symmetry boundary conditions # symmetry boundary conditions ... @@ -135,12 +138,12 @@ and the problem is solved:: ... @@ -135,12 +138,12 @@ and the problem is solved:: DirichletBC(V.sub(1), Constant(0.), facets, 3)] DirichletBC(V.sub(1), Constant(0.), facets, 3)] solve(a == l, u, bc) solve(a == l, u, bc) import matplotlib.pyplot as plt import matplotlib.pyplot as plt p = plot(sigma(u)[1,1]/T[1], mode='color') p = plot(sigma(u)[1,1]/T[1], mode='color') plt.colorbar(p) plt.colorbar(p) plt.title(r"$\sigma_{yy}$",fontsize=26) plt.title(r"$\sigma_{yy}$",fontsize=26) The :math:\sigma_{xx} and :math:\sigma_{yy} components should look like The :math:\sigma_{xx} and :math:\sigma_{yy} components should look like that: that: ... ...
 .. raw:: html

.. _ModalAnalysis: .. _ModalAnalysis: ... ...
 .. raw:: html

.. _ReissnerMindlinDG: .. _ReissnerMindlinDG: ... ...
 .. raw:: html

 ... @@ -11,10 +11,10 @@ What is it about ? ... @@ -11,10 +11,10 @@ What is it about ? ------------------------ ------------------------ These numerical tours will introduce you to a wide variety of topics in These numerical tours will introduce you to a wide variety of topics in computational continuum and structural mechanics using the finite element software FEniCS. computational continuum and structural mechanics using the finite element software FEniCS. Many covered topics can be considered as standard and will help the reader in Many covered topics can be considered as standard and will help the reader in getting started with FEniCS using solid mechanics examples. getting started with FEniCS using solid mechanics examples. Other topics will also be more exploratory and will reflect currently investigated research topics, Other topics will also be more exploratory and will reflect currently investigated research topics, illustrating the versatility of FEniCS. illustrating the versatility of FEniCS. ... @@ -22,33 +22,40 @@ illustrating the versatility of FEniCS. ... @@ -22,33 +22,40 @@ illustrating the versatility of FEniCS. The full set of demos can be obtained from the *COmputational MEchanics Toolbox* (COMET) available at The full set of demos can be obtained from the *COmputational MEchanics Toolbox* (COMET) available at https://gitlab.enpc.fr/jeremy.bleyer/comet-fenics. https://gitlab.enpc.fr/jeremy.bleyer/comet-fenics. -------------------- Citing and license -------------------- If you find these demos useful for your research work, please consider citing them using the following If you find these demos useful for your research work, please consider citing them using the following Zenodo DOI https://doi.org/10.5281/zenodo.1287832 Zenodo DOI https://doi.org/10.5281/zenodo.1287832 .. code-block:: none .. code-block:: none @article{bleyer2018numericaltours, @article{bleyer2018numericaltours, title={Numerical Tours of Computational Mechanics with FEniCS}, title={Numerical Tours of Computational Mechanics with FEniCS}, DOI={10.5281/zenodo.1287832}, DOI={10.5281/zenodo.1287832}, publisher={Zenodo}, publisher={Zenodo}, author={Jeremy Bleyer}, author={Jeremy Bleyer}, year={2018}} year={2018}} All this work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License _ |license|. .. |license| image:: https://i.creativecommons.org/l/by-sa/4.0/88x31.png ----------------------- ----------------------- How do I get started ? How do I get started ? ----------------------- ----------------------- You can find instructions on how to install FEniCS on the FEniCS project website You can find instructions on how to install FEniCS on the FEniCS project website http://fenicsproject.org. In the following numerical tours, we will use the http://fenicsproject.org. In the following numerical tours, we will use the Python interface for the different FEniCS scripts. These demos have been written Python interface for the different FEniCS scripts. These demos have been written using FEniCS 2017.2.0 but many should work with older versions. using FEniCS 2017.2.0 but many should work with older versions. FEniCS is also distributed along with an important number of documented or FEniCS is also distributed along with an important number of documented or undocumented examples, some of them will be revisited in these tours but do not undocumented examples, some of them will be revisited in these tours but do not hesitate over looking at other interesting examples. hesitate over looking at other interesting examples. In the following, we will assume that readers possess basic knowledge of FEniCS commands. In the following, we will assume that readers possess basic knowledge of FEniCS commands. In particular, we advise you to go first through the documentation and tutorials https://fenicsproject.org/tutorial/ In particular, we advise you to go first through the documentation and tutorials https://fenicsproject.org/tutorial/ if this is not the case. if this is not the case. ... @@ -56,9 +63,9 @@ if this is not the case. ... @@ -56,9 +63,9 @@ if this is not the case. About the author About the author ---------------------- ---------------------- Jeremy Bleyer _ is a researcher Jeremy Bleyer _ is a researcher in Solid and Structural Mechanics at Laboratoire Navier _, in Solid and Structural Mechanics at Laboratoire Navier _, a joint research unit of Ecole Nationale des Ponts et Chaussées _, a joint research unit of Ecole Nationale des Ponts et Chaussées _, IFSTTAR _ and CNRS _ (UMR 8205). IFSTTAR _ and CNRS _ (UMR 8205). email: jeremy.bleyer@enpc.fr email: jeremy.bleyer@enpc.fr ... ...
 ... @@ -51,8 +51,8 @@ ... @@ -51,8 +51,8 @@

Elasto-plastic analysis of a 2D von Mises material

Elasto-plastic analysis of a 2D von Mises material

Introduction

Introduction

This example is concerned with the incremental analysis of an elasto-plastic

This example is concerned with the incremental analysis of an elasto-plastic ... @@ -410,6 +410,7 @@ when considering a zero hardening modulus.

... @@ -410,6 +410,7 @@ when considering a zero hardening modulus.

 ... @@ -51,7 +51,7 @@ ... @@ -51,7 +51,7 @@

2D linear elasticity

2D linear elasticity

Introduction

Introduction

... @@ -220,6 +220,7 @@ writing/reading. Prefered output format is now .xdmf:

... @@ -220,6 +220,7 @@ writing/reading. Prefered output format is now .xdmf:

 ... @@ -217,7 +217,7 @@ div.nboutput div.output_area.stderr { ... @@ -217,7 +217,7 @@ div.nboutput div.output_area.stderr { .ansi-bold { font-weight: bold; } .ansi-bold { font-weight: bold; } .ansi-underline { text-decoration: underline; } .ansi-underline { text-decoration: underline; }