Hi,
I am solving a solid mechanics problem - deflection of cantilever beam with a membrane embedded. The problem is as shown in figure below,
Problem definition"
The grey colored block is the solid, and the membrane is along the green plane. I impose displacement = 0 bc at the left face, and a neuman traction bc on the top face of block (load).
Governing equation:
The standard solid mechanics equation for this problem would look like, (at steady state)
+ load term
Now, I will add the membrane contribution as follows, (adding red term)
+ load term
where FS = deformation gradient*2nd Piola Kirchoff stress = First Piola Kirchoff stress. In my case I have used Neo Hookean material model to derive stresses.
What I tried:
Initially I solved the simple solid problem without membrane and found no issue. Now, when i add the membrane term (red) I am facing issues. For debugging purpose, I have made the membrane term to be ‘zero’ as following,
I = Identity(dim)
K_m = Constant(0.0)*I
res_solid = inner(F*S,grad(du))*dx - inner(f,du)*ds(24) + inner(K_m,grad(du)("+"))*dS
In the above line, since I made K_m = zero, i expect the results to be same as standard problem, as if there is no membrane. But, I get different results (I am looking at the LHS matrix of assembled form for comparison). I do not understand what is happening here.
Can anyone help me with this?
Thanks.
MWE:
from dolfin import *
import numpy as np
from scipy.sparse import csr_matrix
from mshr import *
# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(16,16,500), 2, 4, 10)
# Define subdomains
class Domain1(SubDomain):
def inside(self, x, on_boundary):
# Assuming the interface is at the halfway point in the x-direction
return x[1] <= 8
class Domain2(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 8
domain1 = Domain1()
domain2 = Domain2()
# Mark subdomains with different markers
domains = MeshFunction('size_t', mesh, mesh.topology().dim())
domains.set_all(0)
domain1.mark(domains, 1)
domain2.mark(domains, 2)
# Define and mark boundaries
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[2], 0)
class TopBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 16)
class Interface(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 8)
left_boundary = LeftBoundary()
top_boundary = TopBoundary()
interface = Interface()
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
left_boundary.mark(boundaries, 23) # Mark the left boundary as 23
top_boundary.mark(boundaries, 24) # Mark the top boundary as 24
interface.mark(boundaries, 25) # Mark the interface as 25
# Define measures for integrating over subdomains and boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dS = Measure('dS', domain=mesh, subdomain_data=boundaries) # for interfaces
nsd = mesh.geometry().dim()
I = Identity(nsd)
# define function space
V = VectorFunctionSpace(mesh, "CG", 1)
u = Function(V)
du = TestFunction(V)
# material properties
nu = Constant(0.3) # poisson ratio
mu = Constant(3.84615) # shear modulus
E = 2*mu*(1+nu) # youngs modulus
kappa = E*nu/((1+nu)*(1-2*nu)) + (2/3)*E/(2*(1 + nu))
rho = Constant(pow(10,-6))
# solid model
F = I + grad(u)
J = det(F)
C = F.T*F
E = 0.5*(C-I)
f = Constant((0,0.00005,0)) # traction force on top face
S = mu*(J**(-2/3))*(I-(tr(C) + (3-nsd))/3*inv(C)) + kappa/2*((J**2)-1)*inv(C) # neo-hookean
# Boundary conditions
bcs = [
DirichletBC(V, Constant((0,0,0)), boundaries, 23),
]
# define the variational problem
res_solid = inner(F*S,grad(du))*dx - inner(f,du)*ds(24)
K_m = Constant(0.0)*I
res_solid = res_solid + inner(K_m,grad(du)("+"))*dS(25) # line i am talking about
Dres_solid = derivative(res_solid,u) # Jacobian
# problem
problem_m = NonlinearVariationalProblem(res_solid,u,bcs,Dres_solid)
solver_m = NonlinearVariationalSolver(problem_m)
solver_m.parameters['newton_solver']\
['maximum_iterations'] = 15
solver_m.parameters['newton_solver']\
['relative_tolerance'] = 1e-6
solver_m.parameters['newton_solver']['linear_solver'] = 'mumps'
# following lines are for debugging
A, b = assemble_system(Dres_solid, -res_solid, bcs)
A_csr = csr_matrix(A.array())
np.savetxt("AK.csv", A_csr.toarray(), delimiter=",")
b_np = b.get_local()
np.savetxt("bK.csv", b_np, delimiter=",")
# solve
solver_m.solve()
# output u
ufile = File("u.pvd")
ufile << u
Note: Paraview files might visually look same for both cases, but when i look into my matrices, I can see significant differences.