### Linear Elasticity Using Discontinuous Galerkin

237
views
2
6 months ago by
Hello,

I am trying to simulate an elastic bar deforming under its own weight and a load. I am using the discontinuous Galerkin method. Here are the relevant equations:

The expression for stress:
$\boldsymbol{\sigma} = 2G\triangledown^s \textbf{u} + 2G \frac{\nu}{1-2\nu} (\bigtriangledown \cdot \textbf{u})\textbf{I}$

The differential equation to be solved:
$\bigtriangledown \cdot \boldsymbol{\sigma} + \rho \textbf{g} = \textbf{0}$

The discontinuous Galerkin weak form of the above differential equation:

$\Sigma_{E \epsilon \varepsilon_h} \int_{\Omega_E} \boldsymbol{\sigma}:\bigtriangledown \boldsymbol{v} d\Omega_E - \Sigma_{e \epsilon \Gamma_h \cup \Gamma_D } \int_{\Gamma_E}\{\boldsymbol{\sigma n}\} \cdot [\boldsymbol{v}] d\Gamma_E - \\ \Sigma_{e \epsilon \Gamma_N} \int_{\Gamma_E} \boldsymbol{t}\cdot \boldsymbol{v} d\Gamma_E -\Sigma_{E \epsilon \varepsilon_h} \int (\rho \boldsymbol{g}) \cdot \boldsymbol{v} d\Omega_E + \Sigma_{e \epsilon \Gamma_h \cup \Gamma_D} \frac{\sigma_0}{|e|^{\beta_0}}\int_{\Gamma_E} [\boldsymbol{u}] \cdot [\boldsymbol{v}]d\Gamma_E \\ - \Sigma_{e \epsilon \Gamma_D} \frac{\sigma_0}{|e|^{\beta_0}}\int_{\Gamma_E} \boldsymbol{u_D} \cdot \boldsymbol{v} d\Gamma_E = \boldsymbol{0}$

The boundary conditions are:
• Load on the top (y = height)
• Fixed displacement on the left ($u_x = 0.0$ ).
• Fixed displacement on the bottom ($u_y = 0.0$ ).
• Traction-free boundary on the right.

I am currently getting incorrect results. Is there a problem below with the way the boundary conditions have been specified in FEniCS? Please note that I am trying to impose all the boundary conditions weakly.

Also, does FEniCS support strongly imposed Dirichlet boundary conditions for DG methods?

Here is my code. Thank you for helping:

from fenics import *
import numpy as np
import ufl
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pdb

g_const = -9.81
g = Constant((0.0,g_const))
density_const = 3000
density = Constant(density_const)
G = 10.0E9
nu = 0.0
width = 1.0
height = 3.0

num_x_elements = 70
num_y_elements = 10
num_x_nodes = num_x_elements + 1
num_y_nodes = num_y_elements + 1
total_nodes = num_x_nodes*num_y_nodes

mesh = RectangleMesh(Point(0.0,0.0),Point(width,height),num_x_elements,num_y_elements)

element = VectorElement('DG',triangle,1)
V = FunctionSpace(mesh,element)

top = CompiledSubDomain('near(x[1],{})'.format(height))
bottom = CompiledSubDomain('near(x[1],0.0)')
right = CompiledSubDomain('near(x[0],{})'.format(width))
left = CompiledSubDomain('near(x[0],0.0)')

boundary_markers = FacetFunction('size_t',mesh)
top.mark(boundary_markers,0)
right.mark(boundary_markers, 1)
bottom.mark(boundary_markers,2)
left.mark(boundary_markers, 3)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

ux_left = Constant(0.0)
uy_bottom = Constant(0.0)
traction = Expression(('0.0','near(x[1],height) ? load : 0.0'), degree=1, height=height, load=load)

v = TestFunction(V)
(v_1, v_2) = split(v)

u_ = Function(V)
(u_x, u_y) = split(u_)

d = 2

sigma_0 = Constant(100000.0)
beta_0_1 = Constant(3.0/(d-1.0))
e = FacetArea(mesh)
n = FacetNormal(mesh)

def sigma(u):
return 2.0*G*strain + 2.0*G*nu/(1.0-2.0*nu)*div(u_)*Identity(d)

Eqn = inner(sigma(u_),grad(v))*dx \
- dot(avg(sigma(u_)*n),jump(v))*dS \
- dot(sigma(u_)*n,v)*(ds(2)+ds(3)) \
- dot(traction,v)*(ds(0)+ds(1)) \
- dot(density*g, v)*dx \
+ sigma_0/pow(avg(e),beta_0_1)*dot(jump(u_),jump(v))*dS \
+ sigma_0/pow(e,beta_0_1)*dot(u_,v)*(ds(2)+ds(3)) \
- sigma_0/pow(e,beta_0_1)*ux_left*v_1*ds(3) \
- sigma_0/pow(e,beta_0_1)*uy_bottom*v_2*ds(2)

J = derivative(Eqn,u_)
problem = NonlinearVariationalProblem(Eqn, u_, [], J)
solver = NonlinearVariationalSolver(problem)
newton_solver_parameters = {"nonlinear_solver":"newton",
"newton_solver":{"error_on_nonconvergence":False,
"report":True,
"maximum_iterations":10,
"linear_solver":"mumps",
"convergence_criterion":"residual",
"relaxation_parameter":0.5},
"print_matrix":False,
"print_rhs":False}
solver.parameters.update(newton_solver_parameters)
solver.solve()​
Community: FEniCS Project