I am solving a system of PDEs, where one of the unknowns, a vector field q, must satisfy the boundary condition. n is the normal vector to the surface.
q⋅n=1
How to apply boundary condition of q?Any help will be thanked.
As i said in the previous post, you need to supply the fiberinnhold PDE to anyone to help you. This is because there are different methods for different PDEs (Lagrange multipliers or Nitsches method). It is also dependent of the shape of the domain, as you can restrict it with Dirichlet conditions if your grid is aligned with the coordinate system.
Here is my program. The domain is a square with crack inside,I want to simulate hydraulic fracture.The third PDE is about this question. h is a fluid mass flux vector,n is normal vector.dot(h,n)=0.02.I do not know how to apply this kind of boundary.Thank you for your patience.
Preliminaries and mesh
from dolfin import *
mesh = Mesh(‘mesh.xml’)
Define Space
V = FunctionSpace(mesh, ‘CG’, 1)
W = VectorFunctionSpace(mesh, ‘CG’, 1)
WW = FunctionSpace(mesh, ‘DG’, 0)
WWW = TensorFunctionSpace(mesh, ‘CG’, 1)
d,q = TrialFunction(V),TestFunction(V)
u,v = TrialFunction(W),TestFunction(W)
h,k = TrialFunction(W),TestFunction(W)
Introduce manually the material parameters
mu = 98000000.0
g = Constant((0,-9.8))
l = 1
lmbda = 147000000.0
rouf = 1000.0
M = 100000000.0
m0 = 700.0
b = 1.0
k1 = 2.11e-13
yita_ = 0.001
yita = 0.0
yipuxilong = 50.0
Gc = 2700
t = 0
deltaT = 0.1
tol = 1e-3
Constituive functions
def epsilon(u):
return sym(grad(u))
def sigma(u):
return 2.0muepsilon(u)+lmbdatr(epsilon(u))Identity(len(u))
def psi(u):
return 0.5(lmbda+mu)(0.5*(tr(epsilon(u))+abs(tr(epsilon(u)))))2+mu*inner(dev(epsilon(u)),dev(epsilon(u)))
def H(uold,unew,Hold):
return conditional(lt(psi(uold),psi(unew)),psi(unew),Hold)
def norm1(d):
return sqrt(dot(grad(d),grad(d)))
def n1(d):
return grad(d)/norm1(d)
def w(d):
return dot(n1(d),dot(epsilon(uold),n1(d)))*1e-2
def permeability(d):
return rouf2k1Identity(2)+pow(d,yipuxilong)rouf**2(w(d)**2/(12yita_)-k1)(Identity(2)-outer(n1(d),n1(d)))
def pressure(u):
return -M * (b * tr(epsilon(u))-mold/rouf)
def pressurenew(m):
return -M * (b * tr(epsilon(unew))-m/rouf)
Boundary conditions
top = CompiledSubDomain(“near(x[1], 0.5) && on_boundary”)
bot = CompiledSubDomain(“near(x[1], -0.5) && on_boundary”)
left = CompiledSubDomain(“near(x[0], -0.5) && on_boundary”)
right = CompiledSubDomain(“near(x[0], 0.5) && on_boundary”)
cracktop = CompiledSubDomain(“x[1]<=0.00051 and x[1]>=0.00049 and abs(x[0]) <= 0.1 && on_boundary”)
crackbot = CompiledSubDomain(“x[1]>=-0.00051 and x[1]<=-0.00049 and abs(x[0]) <= 0.1 && on_boundary”)
crackleft = CompiledSubDomain(“x[0]>=-0.11 and x[0]<=-0.09 and abs(x[1]) <= 0.0005 && on_boundary”)
crackright = CompiledSubDomain(“x[0]<=0.11 and x[0]>=0.09 and abs(x[1]) <= 0.0005 && on_boundary”)
def Crack(x):
return abs(x[1]) <= 1e-3 and abs(x[0]) <= 0.1
bcbot = DirichletBC(W, Constant((0.0,0.0)), bot)
bctop = DirichletBC(W, Constant((0.0,0.0)), top)
bcleft = DirichletBC(W, Constant((0.0,0.0)), left)
bcright = DirichletBC(W, Constant((0.0,0.0)), right)
bc_u = [bcbot, bctop,bcleft,bcright]
bc_d = [DirichletBC(V, Constant(1.0), Crack)]
bccracktop = DirichletBC(W.sub(1), 0.002, cracktop)
bccrackbot = DirichletBC(W.sub(1), -0.002, crackbot)
bccrackleft = DirichletBC(W.sub(0), -0.002, crackleft)
bccrackright = DirichletBC(W.sub(0), 0.002, crackright)
bc_h = [bccracktop, bccrackbot,bccrackleft,bccrackright]
m_d = Expression(‘0.0’,degree = 0)
mold = project(m_d,V)
unew,uold = Function(W),Function(W)
dnew,dold,Hold = Function(V),Function(V),Function(WW)
hnew,hold = Function(W),Function(W)
Variational form
#E_u = ((1.0-dold)**2)*inner(grad(v),sigma(u))*dx
E_u = ((1.0-dold)**2)inner(grad(v),sigma(u))dx+inner((Mb(tr(epsilon(u))-mold/rouf)Identity(2)),grad(v))dx-(m0+mold)inner(g,v)dx
E_d = Gclinner(grad(d),grad(q))dx+(Gcd/l-(1-d)H(uold,unew,Hold))qdx
#E_d = (Gclinner(grad(d),grad(q))+((Gc/l)+2.0H(uold,unew,Hold))inner(d,q)-2.0H(uold,unew,Hold)q)dx
E_h = dot((grad(-M/rouf(btr(epsilon(uold))-mold/rouf))-grad(dot(((SpatialCoordinate(mesh))+ uold),g))+dot(inv(permeability(dold)),h)),k)*ds
p_u = LinearVariationalProblem(lhs(E_u), rhs(E_u), unew, bc_u)
p_d = LinearVariationalProblem(lhs(E_d), rhs(E_d), dnew, bc_d)
p_h = LinearVariationalProblem(lhs(E_h), rhs(E_h), hnew, bc_h)
solver_u = LinearVariationalSolver(p_u)
solver_d = LinearVariationalSolver(p_d)
solver_h = LinearVariationalSolver(p_h)
Initialization of the iterative procedure and output requests
conc_d = File ("./ResultsDirnew11/d.pvd")
conc_h = File ("./ResultsDirnew21/h.pvd")
conc_p = File ("./ResultsDirnew31/p.pvd")
conc_m = File ("./ResultsDirnew41/m.pvd")
conc_perm = File ("./ResultsDirnew51/perm.pvd")
conc_u = File ("./ResultsDirnew61/u.pvd")
fname = open(‘ForcevsDisp.txt’, ‘w’)
Staggered scheme
while t<=150.0 :
t += deltaT
iter = 0
err = 1
while err > tol:
iter += 1
solver_u.solve()
solver_d.solve()
solver_h.solve()
err_disp = errornorm(unew,uold,norm_type = 'l2',mesh = None)
err_phi = errornorm(dnew,dold,norm_type = 'l2',mesh = None)
err_flux = errornorm(hnew,hold,norm_type = 'l2',mesh = None)
err = max(err_disp,err_phi,err_flux)
uold.assign(unew)
dold.assign(dnew)
hold.assign(hnew)
mold -= deltaT*div(hnew)
Hold.assign(project(psi(unew),WW))
pnew = project(pressurenew(mold),V)
mnew = project(mold,V)
perm = project(permeability(dnew),WWW)
if err < tol:
print ('Iterations:', iter, ', Total time', t)
if round(t*1e4) % 10 == 0:
conc_d << dnew
conc_h << hnew
conc_p << pnew
conc_m << mnew
conc_perm << perm
conc_u << unew
fname.close()
print (‘Simulation completed’)
First of all, please format your code using ```.
Secondly, you should make the code as minimal as possible, as it makes it easier for people to help you, and is more benifical for the community. I would suggest you only add a working code solving the third PDE.
You also need to specify which boundary/boundaries you would like to apply inner(q,n)=1.
Sorry about that.I want apply inner(h,n)=1 into crack defined in the program.
# Preliminaries and mesh
from dolfin import *
mesh = Mesh('mesh.xml')
# Define Space
W = VectorFunctionSpace(mesh, 'CG', 1)
h,k = TrialFunction(W),TestFunction(W)
# Introduce manually the material parameters
#...
# Constituive functions
#...
# Boundary conditions
def Crack(x):
return abs(x[1]) <= 1e-3 and abs(x[0]) <= 0.1
hnew,hold = Function(W)，Function(W)
E_h = dot(h,k)*dx
p_h = LinearVariationalProblem(lhs(E_h), rhs(E_h), hnew, bc_h)
solver_h = LinearVariationalSolver(p_h)
# Initialization of the iterative procedure and output requests
#...
fname = open('ForcevsDisp.txt', 'w')
# Solving scheme
while t<=150.0 :
t += deltaT
iter = 0
err = 1
while err > tol:
iter += 1
solver_h.solve()
err= errornorm(hnew,hold,norm_type = 'l2',mesh = None)
hold.assign(hnew)
if err < tol:
print ('Iterations:', iter, ', Total time', t)
if round(t*1e4) % 10 == 0:
conc_h << hnew
print ('Simulation completed')
As far as I can tell, the crack here is an actual 2D domain (depends on the mesh size), not a 1D domain. As this will be an internal boundary (aligned with the mesh coordinates if your mesh is structured), I guess you can use a DirichletCondition, combined with a MeshFunction and SubDomain marker, see this tutorial. To let it mark the specific internal boundaries create your crack as
class Crack(SubDomain):
def inside(self, x, on_boundary):
return abs(x[1]) <= 1e-3 and abs(x[0]) <= 0.1
# the mark mesh function
Thanks a lot.I will study the tutorial first and have a try.
After mark the specific internal boundary, I still don’t know how to apply inner(h,n)=1 on the crack,like this:bcbot = DirichletBC(W, Constant((0.0,0.0)), bot)? But I don’t know the specific value of vector h.
As your crack goes along the x-axis, n=(0,1). This means that you can use W.sub(1), and a Constant(0,1) in the bd
What should I do if the crack is a circle? Then the n is not specific.
And I don’t know why I need to define the subdomain.
If the crack is a circle you need to use Nitsches method, Lagrange multipliers or the mixed dimensional finite element method.
In any case you need to define the subdomain such that you can mark a MeshFunction, which in turn can be used either in
a) A DirichletBC as described above
b) dS_crack =Measure(“dS”, domain=mesh, subdomain_data=mf, subdomain_id=crack_marker) which can be used for integration over this facet for Nitsche or Lagrange multipliers.
c) Create a submesh for the mixed dimensional FEM
I just started using this software.Do you have similar boundary condition example to share?
I applied boundary condition in this way.But it seems wrong.I don’t know why.And do you know when to use mixed element?
You need to supply a minimal example for anyone to help you. There are plenty of demos for FEniCS at bitbucket, for instance this mixed poisson demo and poisson with pure Neumann conditions
Please follow Read before posting: How do I get my question answered?
which describes how to make a minimal working example.
It is a coupled problem,so it is a little difficult to make a minimal working example.This nearly the simplest program.test program
There are always ways for making minimal examples. You cannot expect that other people have time to debug your whole code for you.
Anyhow, your link to google drive is not publicly available.
The first step is to see if you can solve each sub problem of your code and see if you can apply your boundary conditions on a simpler problem with no coupling. Then, if that works, you can try to couple the two problems.
OK.Thank you very much