(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 = ( \
	inner(grad(vr), mu*(grad(ur) + grad(ur).T)) 			\
	- 2.0/3.0*mu*div(vr)*div(ur)					\
	- inner(vr, rom2*ur) - div(vr)*pr 					\
 	- qr*div(ur) 							\
 	- tau1*inner(-grad(qr) + rom2*vr, -grad(pr) + rom2*ur)					\
	+ 0.5*tau2*rom2*rom2*inner(grad(vr) - grad(vr).T, grad(ur) - grad(ur).T )   )*dx \
	+ ks*mu*(2.0*cn*vr[0]*ui[0] + (cn*cn-sn*sn)/cn*vr[1]*ui[1])*ds(1)

ai = ( \
	inner(grad(vi), mu*(grad(ui) + grad(ui).T)) 			\
	- 2.0/3.0*mu*div(vi)*div(ui)					\
	- inner(vi, rom2*ui) - div(vi)*pi 					\
 	- qi*div(ui) 							\
 	- tau1*inner(-grad(qi) + rom2*vi, -grad(pi) + rom2*ui)					\
	+ 0.5*tau2*rom2*rom2*inner(grad(vi) - grad(vi).T, grad(ui) - grad(ui).T )   )*dx \
	- 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

1 Answer


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)​
Please login to add an answer/comment or follow this question.
The thread is closed. No new answer/comment may be added.

Similar posts:
Search »
  • Nothing matches yet.