Mixed linear and nonlinear elasticity problem

variational formulation
1

Hi
I recently installed the mixed-dimensional branch and I am trying to make sense out of it. I want to solve a problem on a simple 2D geometry. The geometry is a rectangle fixed at the left and right edges under distributed pressure load at the top edge. This rectangle is divided equally into 2 rectangles (left and right subdomains). The left subdomain is a linear elastic material but the right subdomain is a hyperelastic material . This figure shows the problem that I am looking at:

BN
BN904Ă—358 6.99 KB

In other words, I want to solve a linear and a nonlinear problem simultaneously. I have already defined the geometry mesh, marked domains and boundaries and also defined the variational problem on each subdomain but in the last step, I am not sure how to solve it (Maybe because there are new syntax that I am not familiar with). Does anybody know how I can solve it?

Here is the minimal work:

from dolfin import *

mesh = RectangleMesh(Point(0., 0.), Point(10, 2), 20, 5)

n = FacetNormal(mesh)
# Ditributed load at the top
p = Constant(10)

# Define boundaries
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], 10.0)

class Top_Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 2.0) and x[0] >= 0.0 and x[0] <= 5.0

class Top_Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 2.0) and x[0] >= 5.0 and x[0] <= 10.0

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
subdomains = MeshFunction('size_t', mesh, mesh.topology().dim())

left = Left()
right = Right()
top_left = Top_Left()
top_right = Top_Right()

boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
top_left.mark(boundaries, 3)
top_right.mark(boundaries, 4)

# Define subdomains
class LEFT_HALF(SubDomain):
  def inside(self, x, on_boundary):
      return True if x[0] >= 0.0 and x[0] <= 5.0 and x[1] >= 0.0 and x[1] <= 2.0 else False

class RIGHT_HALF(SubDomain):
  def inside(self, x, on_boundary):
      return True if x[0] >= 5.0 and x[0] <= 10.0 and x[1] >= 0.0 and x[1] <= 2.0 else False

left_half = LEFT_HALF()
right_half = RIGHT_HALF()

# Marking left and right subdomains
subdomains.set_all(0)
left_half.mark(subdomains, 1)
right_half.mark(subdomains, 2)

# Linear elasticity on the left domain
def eps(v):
    return sym(grad(v))

E = Constant(1e5)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

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

# Function spaces for the left and right subdomains
V_left = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
V_right = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
W = MixedFunctionSpace(V_left, V_right)

# test and trial functions for the left and right subdomains
du_left, du_right = TrialFunctions(W)
v_left, v_right = TestFunctions(W)

# BCs
bc_left = DirichletBC(W.sub_space(0), Constant((0.,0.)), boundaries , 1)
bc_right = DirichletBC(W.sub_space(1), Constant((0.,0.)), boundaries , 2)
bcs = [bc_left,bc_right]

# Variational form for the left subdomain
a = inner(sigma(du_left), eps(v_left))*dx(1) + inner(sigma(du_left), eps(v_left))*dx(1) # dx(1) corresponds to the left subdomain
l = inner(-p*n, v_left) * ds(3)  # ds(3) corresponds to the top edge of the left subdomain

u_mixed = Function(W)
#solve(a == l, u_mixed, bcs)

# Hyperelasticity on the right subdomain
I = Identity(2)             # Identity tensor
F = I + grad(u_mixed)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

# Total potential energy  ### dx(2) corresponds to the right subdomain
Pi = psi*dx(2) - inner(-p*n, v_right) * ds(4)   # ds(4) corresponds to the top edge of the right subdomain

# directional derivative
F = derivative(Pi, u_mixed, v_right)

# Compute Jacobian of F
J = derivative(F, u_mixed, du_right)

# Solve variational problem
#solve(F == 0, u_mixed, bcs, J=J)

# How to solve it?!
2

You shouldn’t need the mixed dimensional branch for this. You can just solve it with a single displacement field, but assembling different integrands on the subdomains. Also, note that, if you use a monolithic solution scheme (which is my interpretation of “simultaneously”), the coupled problem will itself be a nonlinear one, even if one of the subproblems happens to be linear, so you can essentially ignore the fact that the left subproblem is linear when setting up the coupled problem.

# ...

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

####### Everything above here is the same #######

# Function spaces for the left and right subdomains
#V_left = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
#V_right = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
#W = MixedFunctionSpace(V_left, V_right)

# test and trial functions for the left and right subdomains
#du_left, du_right = TrialFunctions(W)
#v_left, v_right = TestFunctions(W)

# Solve as a nonlinear problem for a single displacement field:
V = VectorFunctionSpace(mesh,"Lagrange",degree=2)
u = Function(V)
v = TestFunction(V)

# BCs
#bc_left = DirichletBC(W.sub_space(0), Constant((0.,0.)), boundaries , 1)
#bc_right = DirichletBC(W.sub_space(1), Constant((0.,0.)), boundaries , 2)
bc_left = DirichletBC(V, Constant((0.,0.)), boundaries , 1)
bc_right = DirichletBC(V, Constant((0.,0.)), boundaries , 2)
bcs = [bc_left,bc_right]

# Variational form for the left subdomain
#a = inner(sigma(du_left), eps(v_left))*dx(1) + inner(sigma(du_left), eps(v_left))*dx(1) # dx(1) corresponds to the left subdomain
#l = inner(-p*n, v_left) * ds(3)  # ds(3) corresponds to the top edge of the left subdomain

# Formulate linear subproblem as contribution to
# nonlinear residual for full problem.  (Why was the
# interior term added twice in the original code?)
F_lin = inner(sigma(u), eps(v))*dx(1) - inner(-p*n, v)*ds(3)

#u_mixed = Function(W)
#solve(a == l, u_mixed, bcs)

# Hyperelasticity on the right subdomain
I = Identity(2)             # Identity tensor
#F = I + grad(u_mixed)             # Deformation gradient
F = I + grad(u)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

# Stored strain energy density (compressible neo-Hookean model)
#psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

# Make the right half softer to show the effect visually:
mu_hyp = Constant(1e-3)*mu
lmbda_hyp = Constant(1e-2)*lmbda
psi = (mu_hyp/2)*(Ic - 3) - mu_hyp*ln(J) + (lmbda_hyp/2)*(ln(J))**2

# Total potential energy  ### dx(2) corresponds to the right subdomain
#Pi = psi*dx(2) - inner(-p*n, v_right) * ds(4)   # ds(4) corresponds to the top edge of the right subdomain

# Note that to get the external force via derivative(),
# you need the displacement, not the test function to
# mulitply the force.
Pi = psi*dx(2) - inner(-p*n, u)*ds(4)

# directional derivative
#F = derivative(Pi, u_mixed, v_right)
F_hyp = derivative(Pi, u, v)

# Combine nonlinear residuals for subproblems:
F = F_lin + F_hyp

# Compute Jacobian of F
#J = derivative(F, u_mixed, du_right)
J = derivative(F, u)

# Solve variational problem
#solve(F == 0, u_mixed, bcs, J=J)
solve(F==0, u, bcs, J=J)

from matplotlib import pyplot as plt
plot(u,mode="displacement")
plt.show()
2 Likes
3

Thanks for your solution and comments. Sorry if I had some typos in my code (Putting the interior term twice in the bilinear part or multiplying the test function instead of function in total potential energy).
That was really helpful.