Hi everyone!
I’m trying to simulate the post-buckling behavior under compresion of a rectangular shell. First of all, I verified my model using a linear solver getting good results, but when I try to use a nonlinear solver, adding some coordinates coupling terms which I found in a tutorial, instability never occures.
I’ve been trying to find the error for a month, but I can’t find what I’m doing wrong.
Here is my example. I tried to reduce it as much as I could:
import matplotlib.pyplot as plt
import numpy as np
import os
from dolfin import *
from ufl import replace
print("-> Geometry and Mesh:")
B = 10.
H = 10.
seeds = 30
Nx = seeds
Ny = seeds
mesh = RectangleMesh(Point(0., 0.), Point(B, H), Nx, Ny, "right")
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], B)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], H)
# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
# ------------------------------------------------------------------------------------------------------------------------------------------
#-> Model Parameters:
b = assemble(1*ds(2))
h = Constant(0.46)
E = Constant(70e3)
nu = Constant(0.3)
lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2 / (1 + nu)
lmbda_ps = 2 * lmbda * mu / (lmbda + 2 * mu)
# Neumann Load:
p0 = 5000. # MPa
p = Expression("p0*t", t=1., p0=p0, degree=0)
# ------------------------------------------------------------------------------------------------------------------------------------------
#-> Element formulation:
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), degree = 1)
P2 = FiniteElement("Lagrange", mesh.ufl_cell(), degree = 2)
bubble = FiniteElement("B", mesh.ufl_cell(), degree = 3)
enriched = P1 + bubble
element = MixedElement([VectorElement(P2, dim=2), P2, VectorElement(enriched, dim=2)])
V = FunctionSpace(mesh, element)
v = Function(V)
u, w, theta = split(v)
v_ = TestFunction(V)
u_, w_, theta_ = split(v_)
dv = TrialFunction(V)
# ------------------------------------------------------------------------------------------------------------------------------------------
#-> Shell Formulation:
eps = sym(grad(u)) + 0.5*outer(grad(w), grad(w))
kappa = sym(grad(theta))
gamma = grad(w) - theta
def plane_stress_elasticity(e):
return lmbda_ps * tr(e) * Identity(2) + 2 * mu * e
N = h * plane_stress_elasticity(eps)
M = h ** 3 / 12 * plane_stress_elasticity(kappa)
Q = mu * h * gamma
# ------------------------------------------------------------------------------------------------------------------------------------------
print("-> Boundary Conditions:")
bi_u = DirichletBC(V.sub(0), Constant((0.,0.)), boundaries, 4) # Borde inferior -> ux = uy = 0
bi_w = DirichletBC(V.sub(1), Constant(0.), boundaries, 4) # Borde inferior -> uz = 0
bs_w = DirichletBC(V.sub(1), Constant(0.), boundaries, 2) # Borde superior -> uz = 0
bli_w = DirichletBC(V.sub(1), Constant(0.), boundaries, 1) # Lat Izquierdo -> uz = 0
bld_w = DirichletBC(V.sub(1), Constant(0.), boundaries, 3) # Lat Derecho -> uz = 0
bcs = [bi_u, bi_w, bs_w, bli_w, bld_w]
# ------------------------------------------------------------------------------------------------------------------------------------------
#-> Variational Formulation and solution:
Wdef = (
inner(N, eps)
+ inner(M, kappa)
+ dot(Q, gamma)
)*dx # Strain energy
Wext = dot(p, u[1])*ds(2) # External work
form = derivative(Wdef - Wext, v, v_) # in-> Formulation, Function, TestFunction
J = derivative(form, v, dv) # in-> Jacobian, Function, TrialFunction
init = Function(V)
v.assign(init)
problem = NonlinearVariationalProblem(form, v, bcs, J = J)
solver = NonlinearVariationalSolver(problem)
solver.parameters["newton_solver"]["absolute_tolerance"] = 1E-8
solver.solve()
u_h, w_h, theta_h = v.split(deepcopy=True)
# ------------------------------------------------------------------------------------------------------------------------------------------
# Post-Process
im = plot(u_h[1], mode="color")
plt.colorbar(im)
plt.show()
Any comment will be very helpful.