reissner_mindlin_dg.py 5.12 KB
Newer Older
1 2 3 4
raw:: html

<a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/"><p align="center"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-sa/4.0/88x31.png"/></a><br />This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/">Creative Commons Attribution-ShareAlike 4.0 International License</a></p>

Jeremy BLEYER's avatar
Jeremy BLEYER committed
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
# .. _ReissnerMindlinDG:
#
# ==============================================================
# Reissner-Mindlin plate with a Discontinuous-Galerkin approach
# ==============================================================
#
# -------------
# Introduction
# -------------
#
# This program solves the Reissner-Mindlin plate equations on the unit
# square with uniform transverse loading and clamped boundary conditions.
# The corresponding file can be obtained from :download:`reissner_mindlin_dg.py`.
#
# It uses a Discontinuous Galerkin interpolation for the rotation field to
# remove shear-locking issues in the thin plate limit. Details of the formulation
# can be found in [HAN2011]_.
#
# 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%
#
#
#
# ---------------
# Implementation
# ---------------
#
#
# Material properties and loading are the same as in :ref:`ReissnerMindlinQuads`::

from __future__ import print_function
from fenics import *

E = Constant(1e3)
nu = Constant(0.3)
thick = Constant(1e-2)
D = E*thick**3/(1-nu**2)/12.
F = E/2/(1+nu)*thick*5./6.
f = Constant(-thick**3)

# The unit square mesh is here divided in triangles and we get the facet MeshFunction for the integration measure :math:`\text{d}s`::

N = 40
mesh = UnitSquareMesh(N, N)
facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
ds = Measure("ds", subdomain_data=facets)

# 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::

We = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
Te = VectorElement("DG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, MixedElement([We, Te]))

# Clamped boundary conditions on the lateral boundary are defined as::

def border(x, on_boundary):
    return on_boundary

bc =  [DirichletBC(V.sub(0), Constant(0.), border)]


# Standard part of the variational form is the same (without full integration)::

def strain2voigt(eps):
    return as_vector([eps[0, 0], eps[1, 1], 2*eps[0, 1]])
def voigt2stress(S):
    return as_tensor([[S[0], S[2]], [S[2], S[1]]])
def curv(u):
    (w, theta) = split(u)
    return sym(grad(theta))
def shear_strain(u):
    (w, theta) = split(u)
    return theta-grad(w)
def bending_moment(u):
    DD = as_tensor([[D, nu*D, 0], [nu*D, D, 0],[0, 0, D*(1-nu)/2.]])
    return voigt2stress(dot(DD,strain2voigt(curv(u))))
def shear_force(u):
    return F*shear_strain(u)

u = Function(V)
u_ = TestFunction(V)
du = TrialFunction(V)


L = f*u_[0]*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
# a stabilization term involing a user-defined parameter :math:`s`::

n = FacetNormal(mesh)
h = CellVolume(mesh)
h_avg = (h('+')+h('-'))/2
stabilization = Constant(10.)

(dw, dtheta) = split(du)
(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 \
   - stabilization*D/h_avg*dot(jump(theta_), jump(dtheta))*dS

# 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
# 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 \
   - 2*stabilization*D/h*dot(theta_, dtheta)*ds

# We then solve for the solution and export the relevant fields to XDMF files ::

solve(a == L, u, bc)

(w, theta) = split(u)

Vw = FunctionSpace(mesh, We)
Vt = FunctionSpace(mesh, Te)
ww = Function(Vw, name="Deflection")
tt = Function(Vt, name="Rotation")
ww.assign(project(w, Vw))
tt.assign(project(theta, Vt))

file_results = XDMFFile("RM_DG_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(ww, 0.)
file_results.write(tt, 0.)

# The solution is compared to the Kirchhoff analytical solution::

print("Kirchhoff deflection:", -1.265319087e-3*float(f/D))
print("Reissner-Mindlin FE deflection:", -ww(0.5, 0.5))

# 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.