Commit bee702af authored by Jeremy BLEYER's avatar Jeremy BLEYER

Some updates of Reissner DG

parent 411bbced
...@@ -10,16 +10,18 @@ Introduction ...@@ -10,16 +10,18 @@ Introduction
------------- -------------
This program solves the Reissner-Mindlin plate equations on the unit This program solves the Reissner-Mindlin plate equations on the unit
square with uniform transverse loading and simply supported boundary conditions. square with uniform transverse loading and clamped boundary conditions.
The corresponding file can be obtained from :download:`reissner_mindlin_dg.py`. The corresponding file can be obtained from :download:`reissner_mindlin_dg.py`.
It uses a Discontinuous Galerkin interpolation for the rotation field to It uses a Discontinuous Galerkin interpolation for the rotation field to
remove shear-locking issues in the thin plate limit. Details of the formulation remove shear-locking issues in the thin plate limit. Details of the formulation
can be found in *P. Hansbo et al., Comput. Methods Appl. Mech. Engrg. 200 (2011) 638–648*, can be found in [HAN2011]_.
https://doi.org/10.1016/j.cma.2010.09.009
The solution for :math:`w` in this demo will look as follows: The solution for :math:`\theta_x` on the middle line of equation :math:`y=0.5`
will look as follows for 10 elements and a stabilization parameter :math:`s=1`:
.. image:: dg_rotation_N10_s1.png
:scale: 15%
...@@ -30,6 +32,7 @@ Implementation ...@@ -30,6 +32,7 @@ Implementation
Material properties and loading are the same as in :ref:`ReissnerMindlinQuads`:: Material properties and loading are the same as in :ref:`ReissnerMindlinQuads`::
from __future__ import print_function
from fenics import * from fenics import *
E = Constant(1e3) E = Constant(1e3)
...@@ -45,37 +48,37 @@ The unit square mesh is here divided in triangles and we get the facet MeshFunct ...@@ -45,37 +48,37 @@ The unit square mesh is here divided in triangles and we get the facet MeshFunct
mesh = UnitSquareMesh(N, N) mesh = UnitSquareMesh(N, N)
facets = MeshFunction("size_t", mesh, 1) facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0) facets.set_all(0)
ds = Measure("ds")[facets] ds = Measure("ds", subdomain_data=facets)
Continuous interpolation using of degree 2 is chosen for the deflection :math:`w` Continuous interpolation using of degree 2 is chosen for the deflection :math:`w`
whereas the rotation field :math:`\underline{\theta}` is discretized using discontinuous linear polynomials:: whereas the rotation field :math:`\underline{\theta}` is discretized using discontinuous linear polynomials::
We = FiniteElement("Lagrange", mesh.ufl_cell(), 2) We = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
Te = VectorElement("DG", mesh.ufl_cell(), 1) Te = VectorElement("DG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh,MixedElement([We,Te])) V = FunctionSpace(mesh, MixedElement([We, Te]))
Clamped boundary conditions on the lateral boundary are defined as:: Clamped boundary conditions on the lateral boundary are defined as::
def border(x, on_boundary): def border(x, on_boundary):
return on_boundary return on_boundary
bc = [DirichletBC(V.sub(0),Constant(0.), border)] bc = [DirichletBC(V.sub(0), Constant(0.), border)]
Standard part of the variational form is the same (without full integration):: Standard part of the variational form is the same (without full integration)::
def strain2voigt(eps): def strain2voigt(eps):
return as_vector([eps[0,0],eps[1,1],2*eps[0,1]]) return as_vector([eps[0, 0], eps[1, 1], 2*eps[0, 1]])
def voigt2stress(S): def voigt2stress(S):
return as_tensor([[S[0],S[2]],[S[2],S[1]]]) return as_tensor([[S[0], S[2]], [S[2], S[1]]])
def curv(u): def curv(u):
(w,theta) = split(u) (w, theta) = split(u)
return sym(grad(theta)) return sym(grad(theta))
def shear_strain(u): def shear_strain(u):
(w,theta) = split(u) (w, theta) = split(u)
return theta-grad(w) return theta-grad(w)
def bending_moment(u): def bending_moment(u):
DD = as_tensor([[D,nu*D,0],[nu*D,D,0],[0,0,D*(1-nu)/2.]]) DD = as_tensor([[D, nu*D, 0], [nu*D, D, 0],[0, 0, D*(1-nu)/2.]])
return voigt2stress(dot(DD,strain2voigt(curv(u)))) return voigt2stress(dot(DD,strain2voigt(curv(u))))
def shear_force(u): def shear_force(u):
return F*shear_strain(u) return F*shear_strain(u)
...@@ -84,10 +87,10 @@ Standard part of the variational form is the same (without full integration):: ...@@ -84,10 +87,10 @@ Standard part of the variational form is the same (without full integration)::
u_ = TestFunction(V) u_ = TestFunction(V)
du = TrialFunction(V) du = TrialFunction(V)
L = f*u_[0]*dx L = f*u_[0]*dx
a = inner(bending_moment(u_),curv(du))*dx + dot(shear_force(u_),shear_strain(du))*dx a = inner(bending_moment(u_), curv(du))*dx + dot(shear_force(u_), shear_strain(du))*dx
We then add the contribution of jumps in rotation across all internal facets plus We then add the contribution of jumps in rotation across all internal facets plus
a stabilization term involing a user-defined parameter :math:`s`:: a stabilization term involing a user-defined parameter :math:`s`::
...@@ -96,40 +99,48 @@ a stabilization term involing a user-defined parameter :math:`s`:: ...@@ -96,40 +99,48 @@ a stabilization term involing a user-defined parameter :math:`s`::
h = CellVolume(mesh) h = CellVolume(mesh)
h_avg = (h('+')+h('-'))/2 h_avg = (h('+')+h('-'))/2
stabilization = Constant(10.) stabilization = Constant(10.)
(dw,dtheta) = split(du) (dw, dtheta) = split(du)
(w_,theta_) = split(u_) (w_, theta_) = split(u_)
a -= dot(avg(dot(bending_moment(u_),n)),jump(dtheta))*dS + dot(avg(dot(bending_moment(du),n)),jump(theta_))*dS \ a -= dot(avg(dot(bending_moment(u_), n)), jump(dtheta))*dS + dot(avg(dot(bending_moment(du), n)), jump(theta_))*dS \
- stabilization*D/h_avg*dot(jump(theta_),jump(dtheta))*dS - stabilization*D/h_avg*dot(jump(theta_), jump(dtheta))*dS
Because of the clamped boundary conditions, we also need to add the corresponding Because of the clamped boundary conditions, we also need to add the corresponding
contributions of the external facets (the imposed rotation is zero on the boundary contributions of the external facets (the imposed rotation is zero on the boundary
so that no term arise in the linear functional):: so that no term arise in the linear functional)::
a -= dot(dot(bending_moment(u_),n),dtheta)*ds + dot(dot(bending_moment(du),n),theta_)*ds \ a -= dot(dot(bending_moment(u_), n), dtheta)*ds + dot(dot(bending_moment(du), n), theta_)*ds \
- 2*stabilization*D/h*dot(theta_,dtheta)*ds - 2*stabilization*D/h*dot(theta_, dtheta)*ds
We then solve for the solution and export the relevant fields to XDMF files :: We then solve for the solution and export the relevant fields to XDMF files ::
solve(a == L, u, bc) solve(a == L, u, bc)
(w,theta) = split(u) (w, theta) = split(u)
Vw = FunctionSpace(mesh,We) Vw = FunctionSpace(mesh, We)
Vt = FunctionSpace(mesh,Te) Vt = FunctionSpace(mesh, Te)
ww = Function(Vw, name="Deflection") ww = Function(Vw, name="Deflection")
tt = Function(Vt, name="Rotation") tt = Function(Vt, name="Rotation")
ww.assign(project(w, Vw)) ww.assign(project(w, Vw))
tt.assign(project(theta, Vt)) tt.assign(project(theta, Vt))
file_results = XDMFFile("RM_DG_results.xdmf") file_results = XDMFFile("RM_DG_results.xdmf")
file_results.parameters["flush_output"] = True file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True file_results.parameters["functions_share_mesh"] = True
file_results.write(ww, 0.) file_results.write(ww, 0.)
file_results.write(tt, 0.) file_results.write(tt, 0.)
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:", ww(0.5,0.5) print("Reissner-Mindlin FE deflection:", -ww(0.5, 0.5))
\ No newline at end of file
For :math:`h=0.001` and 50 elements per side, one finds :math:`w_{FE} = 1.38322\text{e-5}` against :math:`w_{\text{Kirchhoff}} = 1.38173\text{e-5}` for the thin plate solution.
-----------
References
-----------
.. [HAN2011] Peter Hansbo, David Heintz, Mats G. Larson, A finite element method with discontinuous rotations for the Mindlin-Reissner plate model, *Computer Methods in Applied Mechanics and Engineering*, 200, 5-8, 2011, pp. 638-648, https://doi.org/10.1016/j.cma.2010.09.009.
\ No newline at end of file
...@@ -10,13 +10,13 @@ Introduction ...@@ -10,13 +10,13 @@ Introduction
------------- -------------
This program solves the Reissner-Mindlin plate equations on the unit This program solves the Reissner-Mindlin plate equations on the unit
square with uniform transverse loading and fully clamped boundary conditions. square with uniform transverse loading and fully clamped boundary conditions.
The corresponding file can be obtained from :download:`reissner_mindlin_quads.py`. The corresponding file can be obtained from :download:`reissner_mindlin_quads.py`.
It uses quadrilateral cells and selective reduced integration (SRI) to It uses quadrilateral cells and selective reduced integration (SRI) to
remove shear-locking issues in the thin plate limit. Both linear and remove shear-locking issues in the thin plate limit. Both linear and
quadratic interpolation are considered for the transverse deflection quadratic interpolation are considered for the transverse deflection
:math:`w` and rotation :math:`\underline{\theta}`. :math:`w` and rotation :math:`\underline{\theta}`.
The solution for :math:`w` in this demo will look as follows: The solution for :math:`w` in this demo will look as follows:
...@@ -32,11 +32,12 @@ Implementation ...@@ -32,11 +32,12 @@ Implementation
Material parameters for isotropic linear elastic behavior are first defined:: Material parameters for isotropic linear elastic behavior are first defined::
from __future__ import print_function
from fenics import * from fenics import *
E = Constant(1e3) E = Constant(1e3)
nu = Constant(0.3) nu = Constant(0.3)
Plate bending stiffness :math:`D=\dfrac{Eh^3}{12(1-\nu^2)}` and shear stiffness :math:`F = \kappa Gh` Plate bending stiffness :math:`D=\dfrac{Eh^3}{12(1-\nu^2)}` and shear stiffness :math:`F = \kappa Gh`
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`::
...@@ -60,38 +61,38 @@ Continuous interpolation using of degree :math:`d=\texttt{deg}` is chosen for bo ...@@ -60,38 +61,38 @@ Continuous interpolation using of degree :math:`d=\texttt{deg}` is chosen for bo
deg = 1 deg = 1
We = FiniteElement("Lagrange", mesh.ufl_cell(), deg) We = FiniteElement("Lagrange", mesh.ufl_cell(), deg)
Te = VectorElement("Lagrange", mesh.ufl_cell(), deg) Te = VectorElement("Lagrange", mesh.ufl_cell(), deg)
V = FunctionSpace(mesh,MixedElement([We,Te])) V = FunctionSpace(mesh, MixedElement([We, Te]))
Clamped boundary conditions on the lateral boundary are defined as:: Clamped boundary conditions on the lateral boundary are defined as::
def border(x, on_boundary): def border(x, on_boundary):
return on_boundary return on_boundary
bc = [DirichletBC(V,Constant((0.,0.,0.)), border)] bc = [DirichletBC(V, Constant((0., 0., 0.)), border)]
Some useful functions for implementing generalized constitutive relations are now Some useful functions for implementing generalized constitutive relations are now
defined:: defined::
def strain2voigt(eps): def strain2voigt(eps):
return as_vector([eps[0,0],eps[1,1],2*eps[0,1]]) return as_vector([eps[0, 0], eps[1, 1], 2*eps[0, 1]])
def voigt2stress(S): def voigt2stress(S):
return as_tensor([[S[0],S[2]],[S[2],S[1]]]) return as_tensor([[S[0], S[2]], [S[2], S[1]]])
def curv(u): def curv(u):
(w,theta) = split(u) (w, theta) = split(u)
return sym(grad(theta)) return sym(grad(theta))
def shear_strain(u): def shear_strain(u):
(w,theta) = split(u) (w, theta) = split(u)
return theta-grad(w) return theta-grad(w)
def bending_moment(u): def bending_moment(u):
DD = as_tensor([[D,nu*D,0],[nu*D,D,0],[0,0,D*(1-nu)/2.]]) DD = as_tensor([[D, nu*D, 0], [nu*D, D, 0],[0, 0, D*(1-nu)/2.]])
return voigt2stress(dot(DD,strain2voigt(curv(u)))) return voigt2stress(dot(DD,strain2voigt(curv(u))))
def shear_force(u): def shear_force(u):
return F*shear_strain(u) return F*shear_strain(u)
The contribution of shear forces to the total energy is under-integrated using The contribution of shear forces to the total energy is under-integrated using
a custom quadrature rule of degree :math:`2d-2` i.e. for linear (:math:`d=1`) a custom quadrature rule of degree :math:`2d-2` i.e. for linear (:math:`d=1`)
quadrilaterals, the shear energy is integrated as if it were constant (1 Gauss point instead of 2x2) quadrilaterals, the shear energy is integrated as if it were constant (1 Gauss point instead of 2x2)
and for quadratic (:math:`d=2`) quadrilaterals, as if it were quadratic (2x2 Gauss points instead of 3x3):: and for quadratic (:math:`d=2`) quadrilaterals, as if it were quadratic (2x2 Gauss points instead of 3x3)::
...@@ -99,37 +100,37 @@ and for quadratic (:math:`d=2`) quadrilaterals, as if it were quadratic (2x2 Gau ...@@ -99,37 +100,37 @@ and for quadratic (:math:`d=2`) quadrilaterals, as if it were quadratic (2x2 Gau
u_ = TestFunction(V) u_ = TestFunction(V)
du = TrialFunction(V) du = TrialFunction(V)
dx_shear = dx(metadata={"quadrature_degree":2*deg-2}) dx_shear = dx(metadata={"quadrature_degree": 2*deg-2})
L = f*u_[0]*dx L = f*u_[0]*dx
a = inner(bending_moment(u_),curv(du))*dx + dot(shear_force(u_),shear_strain(du))*dx_shear a = inner(bending_moment(u_), curv(du))*dx + dot(shear_force(u_), shear_strain(du))*dx_shear
We then solve for the solution and export the relevant fields to XDMF files :: We then solve for the solution and export the relevant fields to XDMF files ::
solve(a == L, u, bc) solve(a == L, u, bc)
(w,theta) = split(u) (w, theta) = split(u)
Vw = FunctionSpace(mesh,We) Vw = FunctionSpace(mesh, We)
Vt = FunctionSpace(mesh,Te) Vt = FunctionSpace(mesh, Te)
ww = Function(Vw, name="Deflection") ww = Function(Vw, name="Deflection")
tt = Function(Vt, name="Rotation") tt = Function(Vt, name="Rotation")
ww.assign(project(w, Vw)) ww.assign(project(w, Vw))
tt.assign(project(theta, Vt)) tt.assign(project(theta, Vt))
file_results = XDMFFile("RM_results.xdmf") file_results = XDMFFile("RM_results.xdmf")
file_results.parameters["flush_output"] = True file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True file_results.parameters["functions_share_mesh"] = True
file_results.write(ww, 0.) file_results.write(ww, 0.)
file_results.write(tt, 0.) file_results.write(tt, 0.)
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
For ``N=50`` quads per side, one finds :math:`w_{FE} = 1.38182\text{e-5}` for linear quads For :math:`h=0.001` and 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 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. the thin plate solution.
\ No newline at end of file \ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment