Inaccurate solution in linear elasticity simulation with volume forces and quadrature elements(FEniCS 2019.1.0)

Legacy DOLFIN
,
1

Hi everyone,

I started working with fenics 2019.1.0 a few weeks ago. For my programming project I need to use quadrature elements. To illustrate my problem, I have created a small example that solves simple linear elasticity on a unit cube. The last part of the code checks the solution of the displacement field for correctness. To estimate the error, the modulus of elasticity and Poisson’s ratio are calculated from the solution of the displacement field and the difference to the initial value is calculated.

In my initial code I created a mesh with gmsh added the surfaces and volume as physical groups. If a surface force is applied as a boundary condition or if a displacement on the top surface is prescribed, the error is in the range of 10^{-13}, which seems correct to me.

The problem is that if I apply a volume force (e.g. gravity) then the error is in the range of 0.4-1.2\%, which seems a little too high to me. When I change the mesh and generate it in fenics the error is even higher. Since the code works for the first two types of boundary conditions, I could imagine that it is a problem with the definition of the integration measure or the mesh.

Attached is a code example where the mesh is created in fenics.

Does anyone have any idea how to solve this problem?
Thank you very much for your time!

import dolfin as fe
fe.parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)

import numpy as np

n = 1
mesh = fe.UnitCubeMesh(n, n, n)

V = fe.VectorFunctionSpace(mesh, "CG", 2)

# Test- und Trial-functions
u = fe.TrialFunction(V)
v = fe.TestFunction(V)

# solution
u_sol = fe.Function(V,name="Displacement")

# quadrature space
deg_stress = 2
dim = 6

We = fe.VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=dim,quad_scheme='default')
W = fe.FunctionSpace(mesh, We)
sig = fe.Function(W, name="Stress")

# Dirichlet boundary conditions
tol = fe.DOLFIN_EPS

def boundary_Dx(x, on_boundary):
    return on_boundary and (fe.near(x[0], 0, tol))

def boundary_Dy(x, on_boundary):
    return on_boundary and (fe.near(x[1], 0, tol))

def boundary_Dz(x, on_boundary):
    return on_boundary and (fe.near(x[2], 0, tol))

def top(x, on_boundary):
    return on_boundary and (fe.near(x[1], 1, tol))

bc0 = fe.DirichletBC(V.sub(0), fe.Constant(0.), boundary_Dx)
bc1 = fe.DirichletBC(V.sub(1), fe.Constant(0.), boundary_Dy)
bc2 = fe.DirichletBC(V.sub(2), fe.Constant(0.), boundary_Dz)
bc3 = fe.DirichletBC(V.sub(1), fe.Constant(0.2), top)
# bcs = [bc0, bc1, bc2, bc3]
bcs = [bc0, bc1, bc2]

# Linear elasticity material law
E = fe.Constant(210e9)
nu = fe.Constant(0.3)
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))

b = 1/(1+nu)
a = b*(nu/(1-2*nu))
# Modulus = E*C_voigt
C_voigt = fe.as_tensor([
    [a+b, a, a, 0, 0, 0],
    [a, a+b, a, 0, 0, 0],
    [a, a, a+b, 0, 0, 0],
    [0, 0, 0, (1/2)*b, 0, 0],
    [0, 0, 0, 0, (1/2)*b, 0],
    [0, 0, 0, 0, 0, (1/2)*b]
])

# strain tensor
def epsilon(u_):
    return fe.sym(fe.grad(u_))

#stress tensor in voigt-notation
def sigma(u_):
    e = epsilon(u_)
    e_voigt = fe.as_vector([e[0, 0], e[1, 1], e[2, 2], 2 * e[1, 2], 2 * e[0, 2], 2 * e[0, 1]])
    s_voigt = fe.dot(E*C_voigt, e_voigt)
    return fe.as_vector([s_voigt[0], s_voigt[1], s_voigt[2], 2 * s_voigt[3], 2 * s_voigt[4], 2 * s_voigt[5]])

def as_3D_tensor(X):
    return fe.as_tensor([[X[0], X[5]/2, X[4]/2],
                         [X[5]/2, X[1],   X[3]/2],
                         [X[4]/2, X[3]/2, X[2]]])

# integration measure
metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = fe.dx(metadata=metadata)

# weak equation
sig_ = sigma(u)
a_lhs = fe.inner(epsilon(v),as_3D_tensor(sig_) ) * dxm

rho_g = fe.Constant((0.,1e3, 0.))
# rho_g = fe.Constant((0.,0., 0.))
F_Gravitation = fe.dot(rho_g,v) * dxm

f_rhs = F_Gravitation

# solve equation
fe.solve(a_lhs == f_rhs , u_sol, bcs)

    
# check for correctness
check_eps = epsilon(u_sol)
check_sig = as_3D_tensor(sigma(u_sol))

Veps = fe.TensorFunctionSpace(mesh, "DG", degree=0)
eps = fe.Function(Veps, name="Strain")
eps.assign(fe.project(check_eps, Veps))

Vsig = fe.TensorFunctionSpace(mesh, "DG", degree=0)
sig = fe.Function(Vsig, name="Stress")
sig.assign(fe.project(check_sig, Vsig))

e_xx, e_xy, e_xz, e_yx, e_yy, e_yz, e_zx, e_zy, e_zz = eps.split()
s_xx, s_xy, s_xz, s_yx, s_yy, s_yz, s_zx, s_zy, s_zz = sig.split()

e_yy_vec = e_yy.compute_vertex_values()
num_nodes = len(e_yy_vec)
e_yy_bar = np.sum(e_yy_vec) / num_nodes

s_yy_vec = s_yy.compute_vertex_values()
num_nodes = len(s_yy_vec)
s_yy_bar = np.sum(s_yy_vec) / num_nodes

print("-----------")
print("Check E = ", float(E))
E_bar = s_yy_bar/e_yy_bar
print("Mean E calculated: ", E_bar)
print("Error for E in %: ", 100* float(abs(E - E_bar)/E))

print("------")
print("Check nu = ", float(nu))

e_xx_vec = e_xx.compute_vertex_values()
num_nodes = len(e_xx_vec)
e_xx_bar = np.sum(e_xx_vec) / num_nodes

nu_bar = - e_xx_bar/ e_yy_bar

print("Mean calculated nu = ", float(nu_bar))
print("Error for nu in %: ", 100* float(abs(nu - nu_bar)/nu))