Hi all,
I am solving a homogenization problem on a composite consisting of 2 different materials. I have defined the composite with 2 subdomains on one mesh. I just wanted to know how I can imply continuity of stress and displacement at the interface? does it count as a boundary condition in fenics?
1 Like
If the subdomains are topologically-connected in the mesh, then solving for displacement in a CG function space defined on this mesh will strongly enforce displacement continuity at the interface. If you just vary the material properties on the two subdomains within a standard Galerkin weak form, that will also weakly enforce traction compatibility.
To see this, suppose the subdomains are \Omega_1 and \Omega_2, with \Omega = \Omega_1\cup\Omega_2. Galerkinâs method poses the weak problem
\int_{\Omega_1}\pmb{\sigma}_1:\nabla\mathbf{v} + \int_{\Omega_2}\pmb{\sigma}_2:\nabla\mathbf{v} = \int_\Omega\mathbf{f}\cdot\mathbf{v}~~~~~\forall\mathbf{v}\text{ ,}
where \pmb{\sigma}_i is the Cauchy stress in \Omega_i. Integrating by parts on each \Omega_i, one gets
-\int_{\Omega_1}(\nabla\cdot\pmb{\sigma}_1)\cdot\mathbf{v} -\int_{\Omega_2}(\nabla\cdot\pmb{\sigma}_2)\cdot\mathbf{v} + \int_\Gamma(\pmb{\sigma}_1\cdot\mathbf{n}_1 + \pmb{\sigma}_2\cdot\mathbf{n}_2)\cdot\mathbf{v} = \int_\Omega\mathbf{f}\cdot\mathbf{v}\text{ ,}
where \Gamma = \overline{\Omega}_1\cap\overline{\Omega}_2, i.e., the shared boundary or interface (and Iâve assumed that strongly-enforced Dirichlet and/or periodic BCs are enforced on any external boundaries), and \mathbf{n}_i is the outward-facing normal to \Omega_i. Because \mathbf{n}_1 = -\mathbf{n}_2 on \Gamma, the EulerâLagrange equation corresponding to the \Gamma integral is
(\pmb{\sigma}_1 - \pmb{\sigma}_2)\cdot\mathbf{n}_1 = \mathbf{0}\text{ ,}
which is the stress continuity that you want physically, for equilibrium, while the volume integrals correspond to the strong form of -\nabla\cdot\pmb{\sigma}_i = \mathbf{f} on each subdomain.
4 Likes
Thanks a lot for your reply. It was very helpful.
So, as you said, I just used CG function space but I get the following result which is shown in the picture below. I can tell that stress contour is discontinuous. Am I right?
It depends on what scalar quantity youâre plotting as âstressâ. Equilibrium doesnât imply that all components of the stress tensor are continuous; you only expect to have (\pmb{\sigma}_1 - \pmb{\sigma}_2)\cdot\mathbf{n} = \mathbf{0}.
Thatâs right. But this is the contour for y-axis normal stress.
I just simplified my code to check out if I made a mistake in my variational formulation, but I got the same discontinuity in my result. Here is the simplified code that I just wrote and the final result that I got.
Thanks in advance.
for name in dir():
if not name.startswith(â_â):
del globals()[name]
from dolfin import *
import numpy as np
import matplotlib.pyplot as pl
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)
top = Top()
class Left1(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and between(x[1], (0.5, 1))
leftup=Left1()
class Right1(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0) and between(x[1], (0.5, 1))
rightup = Right1()
class Interface(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.5)
interface = Interface()
class Left2(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and between(x[1], (0.0, 0.5))
leftdown = Left2()
class Right2(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0) and between(x[1], (0.0, 0.5))
rightdown = Right2()
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
bottom = Bottom()
#Create meshes
n = 50
mesh = UnitSquareMesh(n, n)
class Solidup(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.5
solidup= Solidup()
class Soliddown(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.5
soliddown = Soliddown()
regions = MeshFunction(âsize_tâ, mesh, mesh.topology().dim())
regions.set_all(0)
solidup.mark(regions, 1)
soliddown.mark(regions, 2)
boundaries = MeshFunction(âsize_tâ, mesh, 1)
boundaries.set_all(0)
top.mark(boundaries, 1)
leftup.mark(boundaries, 2)
rightup.mark(boundaries, 3)
interface.mark(boundaries, 4)
leftdown.mark(boundaries, 5)
rightdown.mark(boundaries, 6)
bottom.mark(boundaries, 7)
#define a submesh, composed of region 0 and 1
solidup = MeshFunction(âsize_tâ, mesh, mesh.topology().dim())
solidup.set_all(0)
solidup.array()[regions.array() == 1] = 1
submeshsolidup = SubMesh(mesh, solidup, 1)
soliddown = MeshFunction(âsize_tâ, mesh, mesh.topology().dim())
soliddown.set_all(0)
soliddown.array()[regions.array() == 2] = 2
submeshsoliddown = SubMesh(mesh, soliddown, 2)
#define new meshfunctions on this submesh with the same values as their original mesh
solidup_regions = MeshFunction(âsize_tâ, submeshsolidup, submeshsolidup.topology().dim())
solidup_regions.set_all(0)
solidup_boundaries = MeshFunction(âsize_tâ, submeshsolidup, 1)
solidup_boundaries.set_all(0)
soliddown_regions = MeshFunction(âsize_tâ, submeshsoliddown, submeshsoliddown.topology().dim())
soliddown_regions.set_all(0)
soliddown_boundaries = MeshFunction(âsize_tâ, submeshsoliddown, 1)
soliddown_boundaries.set_all(0)
W = VectorFunctionSpace(mesh,âCGâ, 2) # Electroelasticity
w = Function(W)
uz = TestFunction(W)
dx = Measure(âdxâ, domain=mesh, subdomain_data=regions)
bcs = [DirichletBC(W, Constant((0.,0.)), bottom)
,DirichletBC(W, Constant((0.,0.4)), top)]
E = Constant(1e5)
nu = Constant(0.3)
model = âplane_stressâ
Landa1=100
Landa2=150
mu = E/2/(1+nu)
lmbda = Enu/(1+nu)/(1-2nu)
if model == âplane_stressâ:
lmbda = 2mulmbda/(lmbda+2*mu)
def sigma(v):
return lmbdatr(grad(v))Identity(2) + 2.0mugrad(v)
Wint =inner(Landa1*sigma(uz), grad(w))dx(1)+ inner(Landa2sigma(uz), grad(w))*dx(2)
solve(Wint==0, w, bcs)
#Plot sigma and u
plt.figure()
p3 = plot(w[1],title=âksiâ)
plt.colorbar(p3)
plt.show()
plt.figure()
ps = plot(grad(w)[1,1],title=ây-axis Normal stress [MPa]â)
plt.colorbar(ps)
plt.show()
First, I believe thereâs a typo in your formulation, since, in the line Wint = ..., it would be more standard for uz to be the Function and w to be the TestFunction, and for the nonlinear solve to have uz as its unknown solution. Second, if w is interpreted as a displacement, grad(w)[1,1] is not the normal stress in the y-direction, it is the extensional strain. (In general, strain would be the symmetrized gradient of displacement, but symmetrization doesnât change diagonal components.) You do not expect strain to be continuous. The stress tensor would be Landa1*sigma(uz) or Landa2*sigma(uz), depending on which subdomain youâre in (assuming the code is changed so that uz is the displacement, rather than the test function).
1 Like
Just to add to what David has already said, if your function is in CG2 (piecewise quadratic approximation), then itâs gradient should be in DG1 (piecewise linear discontinuous). This is what you should see in your results. I donât know right off the bat, but plot may/may not handle DG1 functions. The best way to visualize is by writing them using write_checkpoint to a xdmf file and open in Paraview.
Also, others would be able to help you better, if you enclose your code withing triple fences ``` such that it is reproducible.
1 Like
Youâre totally right. I had some typos in that code, sorry for that. Now, I totally get it. just a further question if you donât mind. As I am new to fenics, I have no idea that how I will be able to plot sigma(uz) in each subdomain individually. I mean how to separate sigma(uz) in each subdomain to multiply to different Landa.
Thanks for your time.
Thanks for your reply. I got the xdmf file and the result which is shown in paraview is same as the one that I got in by using plot. As david said I have to plot uz in in each subdomain to use the right coefficient for that. I just wanted to ask you if you know how to plot sigma(uz) in each subdomain individually. I mean how to separate sigma(uz) in each subdomain to multiply to different Landa.
So, the way the existing plot function works is that most UFL objects you pass as the first argument are projected (in the sense of L^2) onto a CG1 (linear, continuous) finite element space for visualization. (If you pass a CG1 function, it is plotted directly, or, if you pass a DG0 function, that can also be plotted directly.) You can mimic this behavior for a piecewise function that is f1 on subdomain 1 and f2 on subdomain 2 by writing a custom projection function:
def split_project(f1,f2,V):
u = TrialFunction(V)
v = TestFunction(V)
uh = Function(V)
solve(u*v*dx == f1*v*dx(1) + f2*v*dx(2),uh)
return uh
The return value of this will be a Function in the space V that is the L^2 projection of the piecewise function defined by f1 and f2. If V is a CG1 or DG0 space, that can be plotted with plot directly. (Otherwise, it will be projected onto CG1 to plot.) As @bhaveshshrimali pointed out, you can get a higher-fidelity visualization by choosing V to be a higher-degree DG space, but you would need to write it out to a file and plot with Paraview.
An extra note: L^2-projecting a discontinuous function onto a CG space will give oscillatory results; you see this in your original plot of grad(w)[1,1] (which was projected onto CG1, as explained above), where there is some over- and under-shoot near the discontinuity. One way around such oscillatory behavior is to use a lumped-mass L^2 projection, as demonstrated in my answer here; that is also more computationally-efficient than solving a linear system for the consistent L^2 projection.
2 Likes
How about using a DG0 function landa for that takes values Landa1 and Landa2 as you wish; then using landa*sigma(uz) as your function to be (local-) projected onto a DG1 space and finally written to Paraview using write_checkpoint., which in your case should roughly translate to
from dolfin import *
metadata={"quadrature_degree": 2}
def local_project(v,V):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv,v_)*dx(metadata=metadata)
b_proj = inner(v,v_)*dx(metadata=metadata)
solver = LocalSolver(a_proj,b_proj)
solver.factorize()
u = Function(V)
solver.solve_local_rhs(u)
return u
V1 = TensorFunctionSpace(mesh, "DG", 1) #for sigma
Vdg = FunctionSpace(mesh, "DG", 0) # for landa
landa = Function(Vdg)
x = Vdg.tabulate_dof_coordinates()
for i in range(x.shape[0]):
if regions.array()[i] == 1:
landa.vector().vec().setValueLocal(i, 100) # `Landa1`
else:
landa.vector().vec().setValueLocal(i, 150) # `Landa2`
sigma_vis = local_project(landa*sigma(uz), V1)
with XDMFFile(MPI.comm_world, "sigma.xdmf") as fil:
fil.write_checkpoint(sigma_vis, "sigma", 0)
4 Likes
Thanks a lot.That works perfectly.
Great.That works perfectly.
Please, you could post the complete code as bhaveshshrimali commented, all of it encapsulated with ``` !
And have you remade considering the application of a load on the upper domain?
Dear kamensky,
Could you point me to a book or bibliographic reference where you got this approach from the two-domain garlekin method?
The basic idea is probably as old as the structural analysis that modern finite element methods historically grew out of, since structural analysts are happy to use adjacent elements of different stiffness. I based the particular notation I used in my post on this paper about fluidâstructure interaction.
2 Likes