### (Deleted) infinite displacement encountered in elasticity problem

142
views
-1
7 months ago by

I am solving a time-harmonic elasticity problem where I can find the answer (displacements) analytically. It is defined within the code below as "uexact*" with real and imaginary parts. For boundary conditions, I am enforcing the uexact on the left face, and I have Dirichlet BC for u[1] on top/bottom face, and u[2] on front and back. For the right face, I have this Robin BC:

The results (displacements) that I get by solving the equation with Fenics 2016.2.0 are inf when I printed a few of them. I was wondering what is going wrong here. Any help would be appreciated. Here is the full code:

rom dolfin import *

# Material properties:
mu = 1.0				# shear modulus
rom2 = 36*pi*pi				# rho * omega^2 = 36pi^2 gives wavelength = 1/3.
ks = sqrt(rom2/mu)

# Create functions for boundary conditions
top = 0.5
moden = 1.0
sn = moden*pi/(ks*top)
cn = sqrt(1.0-sn*sn)
zero = Constant(0)

# Exact Solution:
uexactr = Expression(("(-1.0*s/c*cos(m*pi*x[1]/height)*sin(c*k*x[0]))", "(sin(m*pi*x[1]/height)*cos(c*k*x[0]))", "0.0"),
m = moden, s = sn, c = cn, k = ks, height = top, degree=1)

uexacti = Expression( ( "(s/c*cos(m*pi*x[1]/height)*cos(c*k*x[0]))", "(sin(m*pi*x[1]/height)*sin(c*k*x[0]))", "0.0"),
m = moden, s = sn, c = cn, k = ks, height = top, degree=1)

# Geometry and boundaries:
top = 0.50
length = 2.0
width  = 0.5
resolution = 8			# Number of elements along x axis per (real or total) wavelength = 2pi/ks.

# Mesh and mesh parameters:
# Rectangle geometry
nel = int(round(resolution*length*ks/(2*pi)))	# Number of elements along x axis.
h = length/nel					# Element size.
nelx = nel					# Number of elements along x axis.
nely = int(round(top/h))			# Number of elements along y axis.

mesh  = BoxMesh(Point(0, 0, 0), Point(length, top, width), nelx, nely, nely)

# Right boundary for boundary condition:
class Right(SubDomain):
def inside(self, x, on_boundary):
return x[0] > length - DOLFIN_EPS

right = Right()

# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
right.mark(boundaries, 1)

ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

# Computational Parameters:
c = cos(ks*h)
alpha = ks*ks*h*h/6.0
tau1 = h*h/mu/2.0
tau2 = ((1+tau1*rom2)*alpha*(2.0+c)/(1.0-c) - 1.0)/(rom2*ks*ks)

# ========================================================================================================
# Define function spaces
p = 1
V_ele = VectorElement("CG", mesh.ufl_cell(), p)
F_ele = FiniteElement("CG", mesh.ufl_cell(), p)
V = FunctionSpace(mesh, MixedElement([V_ele, V_ele, F_ele, F_ele]))

# Define Dirichlet boundaries: top, bottom, right, left:
def lboundary(x):
return (x[0] < DOLFIN_EPS)

def rboundary(x):
return x[0] > length - DOLFIN_EPS

def tbboundary(x):
return x[1] < DOLFIN_EPS or x[1] > top - DOLFIN_EPS

def frontbackboundary(x):
return x[2] > width - DOLFIN_EPS or x[2] < DOLFIN_EPS

def pboundary(x):
return bool(x[1] < DOLFIN_EPS or x[1] > top - DOLFIN_EPS or x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS)

# Enforce Dirichlet condition on y=x1 component of the displacement on top/bottom boundaries:
bctbr = DirichletBC(V.sub(0).sub(1), zero, tbboundary)     # x1 or y component
bctbi = DirichletBC(V.sub(1).sub(1), zero, tbboundary)     # x1 or y component

# Enforce exact solution as Dirichlet condition on left boundary:
bclr = DirichletBC(V.sub(0), uexactr, lboundary)
bcli = DirichletBC(V.sub(1), uexacti, lboundary)

# Enforce zero normal displacement as Dirichlet condition on front+back boundaries:
bcfbr = DirichletBC(V.sub(0).sub(2), zero, frontbackboundary)
bcfbi = DirichletBC(V.sub(1).sub(2), zero, frontbackboundary)

# Collect boundary conditions
# bcs = [bctbr, bctbi, bclr, bcli]
bcs = [bctbr, bctbi, bclr, bcli, bcfbr, bcfbi]

# Define variational problem
(ur,ui,pr,pi) = TrialFunctions(V)
(vr,vi,qr,qi) =  TestFunctions(V)
f = Constant((0, 0, 0))

ar = ( \
- 2.0/3.0*mu*div(vr)*div(ur)					\
- inner(vr, rom2*ur) - div(vr)*pr 					\
- qr*div(ur) 							\
+ ks*mu*(2.0*cn*vr[0]*ui[0] + (cn*cn-sn*sn)/cn*vr[1]*ui[1])*ds(1)

ai = ( \
- 2.0/3.0*mu*div(vi)*div(ui)					\
- inner(vi, rom2*ui) - div(vi)*pi 					\
- qi*div(ui) 							\
- ks*mu*(2.0*cn*vi[0]*ur[0] + (cn*cn-sn*sn)/cn*vi[1]*ur[1])*ds(1)

a = ar + ai

L = inner(vr + tau1*grad(qr), f)*dx

# Compute solution
w = Function(V)
solve(a == L, w, bcs, solver_parameters={"linear_solver": "mumps"})
Community: FEniCS Project

0
6 months ago by
Hi, consider this formulation. I didn't try it on your problem but it worked for me.

V_ele = VectorElement("CG", mesh.ufl_cell(), p)
F_ele = FiniteElement("CG", mesh.ufl_cell(), p)

T  = V_ele*F_ele
V  = FunctionSpace(mesh, T)​