Hello everyone, an unexpected scene occurred when I did a two-phase fluid displacement simulation in 2D space, where the injected fluid would act like a walk, with the top fluid taking a step forward and then stopping or retreating, followed by a similar movement of the bottom fluid.
the mathematical model I used is Cahn-hilliard Navier-stokes equations, I think the equation is fine,but I can’t found the reason to lead this phenomenon.
CODE：
#----------------------------------------------------------------------
# Model parameters
#----------------------------------------------------------------------
# Cells spcies and interface parameters
cequ = 1.
sigma_ow = 0.03 #is σ, surface tension of between water and oil
angle = 60.
epsilon = 1.e-6
# Diff_c a identifier sur les cylindre
#Mw = 1.e-16
Mw = 5.e-15
eta_w = 0.001 #η(water)
eta_o = 0.01 #η(oil)
#----------------------------------------------------------------------
# Geometrical parameters
#----------------------------------------------------------------------
l = sqrt(2)*2.5e-5
L = 2.e-4 #length#width 200um
R_i = 4.5e-5 #radius of pore after etching
RIC = 4.75e-5
#----------------------------------------------------------------------
# Numerical domain
#----------------------------------------------------------------------
R = Rectangle(Point(0., 0.), Point(L, L))
#----------------------------------------------------------------------
r1 = mshr.Polygon([dolfin.Point(0.5*L,-l),dolfin.Point(L+l,0.5*L),dolfin.Point(0.5*L,L+l),dolfin.Point(-l,0.5*L)])
r2 = mshr.Polygon([dolfin.Point(0.5*L,l),dolfin.Point(L-l,0.5*L),dolfin.Point(0.5*L,L-l),dolfin.Point(l,0.5*L)])
R1 = r1 - r2
#----------------------------------------------------------------------
c1 = Circle(Point(0.5*L, 0.), R_i)
c2 = Circle(Point(L, 0.5*L), R_i)
c3 = Circle(Point(0.5*L, L), R_i)
c4 = Circle(Point(0., 0.5*L), R_i)
C1 = c1 + c2 + c3 + c4
#----------------------------------------------------------------------
domain = R - (R - R1 - C1)
mesh = mshr.generate_mesh(domain, 50)
#----------------------------------------------------------------------
# Mesh writting in file
#----------------------------------------------------------------------
file_mesh = File('w1/mesh.pvd')
file_mesh << mesh
#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Time discretization of the problem
#----------------------------------------------------------------------
minute = 60.
TFIN0 = 0.001
dt0 = 0.00001
num_steps = TFIN0 / dt0
print (num_steps)
tps = np.array([]) ; radi = np.array([])
#TFIN1 = TFIN0 + (2 * minute)
TFIN1 = TFIN0 + 60.
dt1 = 0.001
num_steps = (TFIN1 - TFIN0) / dt1
print (num_steps)
#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Definition of relevant bounds
#----------------------------------------------------------------------
# This time with boundary markers to extract subdomains
#----------------------------------------------------------------------
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundary_markers.set_all(0)
tol = 1.0e-7
tolc = 1.0e-7
#----------------------------------------------------------------------
# Define boundary conditions
#----------------------------------------------------------------------
# Left bound
class Boundary_left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0, tol) and (x[1] >= (0.5*L-R_i)) and (x[1] <= (0.5*L+R_i))
bound_left = Boundary_left()
bound_left.mark(boundary_markers, 1)
#----------------------------------------------------------------------
# Top and Bottom bound
class Boundary_left(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], L, tol) and (x[0] >= (0.5*L-R_i)) and (x[0] <= (0.5*L+R_i)) \
or near(x[1], 0, tol) and (x[0] >= (0.5*L-R_i)) and (x[0] <= (0.5*L+R_i))
bound_left = Boundary_left()
bound_left.mark(boundary_markers, 2)
#----------------------------------------------------------------------
# Right bound
class Boundary_right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], L, tol) and (x[1] >= (0.5*L-R_i)) and (x[1] <= (0.5*L+R_i))
bound_right = Boundary_right()
bound_right.mark(boundary_markers, 3)
#----------------------------------------------------------------------
# Left and Right circles
xc1 = 0.
yc1 = 0.5*L
xc2 = L
yc2 = 0.5*L
class Boundary_4(SubDomain):
def inside(self, x, on_boundary):
return (abs(sqrt((x[1]-yc1)*(x[1]-yc1)+(x[0]-xc1)*(x[0]-xc1))-R_i) < tolc and on_boundary) \
or (abs(sqrt((x[1]-yc2)*(x[1]-yc2)+(x[0]-xc2)*(x[0]-xc2))-R_i) < tolc and on_boundary)
bound_4 = Boundary_4()
bound_4.mark(boundary_markers, 4)
#----------------------------------------------------------------------
# Top and Bottom circles
xc3 = 0.5*L
yc3 = 0.
xc4 = 0.5*L
yc4 = L
class Boundary_5(SubDomain):
def inside(self, x, on_boundary):
return (abs(sqrt((x[1]-yc3)*(x[1]-yc3)+(x[0]-xc3)*(x[0]-xc3))-R_i) < tolc and on_boundary) \
or (abs(sqrt((x[1]-yc4)*(x[1]-yc4)+(x[0]-xc4)*(x[0]-xc4))-R_i) < tolc and on_boundary)
bound_5 = Boundary_5()
bound_5.mark(boundary_markers, 5)
#----------------------------------------------------------------------
x1 = -l
y1 = 0.5*L
x2 = 0.5*L
y2 = L+l
x3 = L+l
y3 = 0.5*L
x4 = 0.5*L
y4 = -l
# Outer quadrilateral
class Boundary_6(SubDomain):
def inside(self, x, on_boundary):
r1 = (y2 - y1) * (x[0] - x1) / (x2 - x1) + y1 - x[1]
r2 = (y2 - y3) * (x[0] - x3) / (x2 - x3) + y3 - x[1]
r3 = (y3 - y4) * (x[0] - x4) / (x3 - x4) + y4 - x[1]
r4 = (y4 - y1) * (x[0] - x1) / (x4 - x1) + y1 - x[1]
return (on_boundary and (near(r1, 0.0, tol) and (x[0] >= x1) and (x[0] <= x2))) \
or (on_boundary and (near(r2, 0.0, tol) and (x[0] >= x2) and (x[0] <= x3))) \
or (on_boundary and (near(r3, 0.0, tol) and (x[0] >= x4) and (x[0] <= x3))) \
or (on_boundary and (near(r4, 0.0, tol) and (x[0] >= x1) and (x[0] <= x4)))
bound_6 = Boundary_6()
bound_6.mark(boundary_markers, 6)
#----------------------------------------------------------------------
x5 = l
y5 = 0.5*L
x6 = 0.5*L
y6 = L-l
x7 = L-l
y7 = 0.5*L
x8 = 0.5*L
y8 = l
# Inner quadrilateral
class Boundary_7(SubDomain):
def inside(self, x, on_boundary):
r5 = (y6 - y5) * (x[0] - x5) / (x6 - x5) + y5 - x[1]
r6 = (y6 - y7) * (x[0] - x7) / (x6 - x7) + y7 - x[1]
r7 = (y7 - y8) * (x[0] - x8) / (x7 - x8) + y8 - x[1]
r8 = (y8 - y5) * (x[0] - x5) / (x8 - x5) + y5 - x[1]
return (on_boundary and (near(r5, 0.0, tol) and (x[0] >= x5) and (x[0] <= x6))) \
or (on_boundary and (near(r6, 0.0, tol) and (x[0] >= x6) and (x[0] <= x7))) \
or (on_boundary and (near(r7, 0.0, tol) and (x[0] >= x8) and (x[0] <= x7))) \
or (on_boundary and (near(r8, 0.0, tol) and (x[0] >= x5) and (x[0] <= x8)))
bound_7 = Boundary_7()
bound_7.mark(boundary_markers, 7)
#----------------------------------------------------------------------
# Boundaries writting in file
file_boundary = File('w1/boundary.pvd')
file_boundary << boundary_markers
# Redefiniction dx and ds
dx = Measure("dx", domain = mesh)
ds = Measure("ds", domain = mesh, subdomain_data = boundary_markers)
#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Definition of functional space, trial & test functions
#----------------------------------------------------------------------
# Define function space
CP1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
CV2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
CHS = FunctionSpace(mesh, MixedElement(CP1, CP1, CP1, CV2))
#----------------------------------------------------------------------
# Define function & parameters
dsol = TrialFunction(CHS)
v_m, v_c, q, w = TestFunctions(CHS)
# Solution vector
sol = Function(CHS) # Current solution
sol0 = Function(CHS) # Previous solution
m, c, p, u = split(sol)
m0, c0, p0, u0 = split(sol0)
# Definition of the normal vector
n = FacetNormal(mesh)
x = SpatialCoordinate(mesh)
#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Definition of Initial condition
#----------------------------------------------------------------------
u_init = SphericalInitialConditions(0., 0.5*L, RIC, cequ, degree = 1)
sol.interpolate(u_init)
sol0.interpolate(u_init)
#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Definition of dirichlet boundary conditions
#----------------------------------------------------------------------
v_null = Constant((0.0, 0.0))
zero = Constant(0.)
# Top and Bottom bound vy = 0
bc0 = DirichletBC(CHS.sub(3).sub(1), zero, boundary_markers, 2)
# No-slip boundary condition for velocity
bc1 = DirichletBC(CHS.sub(3), v_null, boundary_markers, 4)
bc2 = DirichletBC(CHS.sub(3), v_null, boundary_markers, 5)
bc3 = DirichletBC(CHS.sub(3), v_null, boundary_markers, 6)
bc4 = DirichletBC(CHS.sub(3), v_null, boundary_markers, 7)
# Inflow boundary condition for velocity
#v_prof = Expression(("3.5e-6*sin((x[1] + 4.5e-5)*(pi/9.e-5))"), degree=2)
v_prof = Expression(("3.5e-5*sin((x[1] - 5.5e-5)*(pi/9.e-05))"), degree=2)
bc5 = DirichletBC(CHS.sub(3).sub(0), zero, boundary_markers, 1)
bc6 = DirichletBC(CHS.sub(3).sub(1), zero, boundary_markers, 1)
# Outflow vy = 0
bc7 = DirichletBC(CHS.sub(3).sub(1), zero, boundary_markers, 3)
# Collect boundary conditions
bcs = [bc0, bc1, bc2, bc3, bc4, bc5, bc6, bc7]
#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Define variational problem
import dolfin as df
I = df.Identity(2)
beta = 1./4.
dfdc = 2*beta*c*(cequ-c)**2 - 2*(cequ-c)*beta*(c**2)
#----------------------------------------------------------------------
# Weak form
#----------------------------------------------------------------------
# Defivition of divrgence and gradient in cylindrycal coordinates
# Wilson-theta
theta = 0.5
# Interface regularization parameter
alfa = (6*sqrt(2))/(cequ**3.)
# Theta-consistent quantities (to ajust no consistently implemented)
m_mid = (1.0-theta)*m0 + theta*m
u_mid = (1.0-theta)*u0 + theta*u #velocity
c_mid = (1.0-theta)*c0 + theta*c #c is w(water), local mass fraction of water 水的质量分数
eta = eta_w + c_mid*(eta_o - eta_w) #NS equations η, η=η(water)*c+η(oil)*(1-c)
# Equation (28)
L0 = c*v_m*dx - c0*v_m*dx + dt0*Mw*dot(grad(m_mid), grad(v_m))*dx + dt0*dot(u_mid, grad(c_mid))*v_m*dx
# Equation (29)
L1 = (1.0/(epsilon*epsilon))*m*(epsilon/(alfa*sigma_ow))*v_c*dx - (1.0/(epsilon*epsilon))*dfdc*v_c*dx - dot(grad(c), grad(v_c))*dx + (Constant((tan(angle*np.pi/180))**-1.0))*abs(dot(as_vector([-1*n[1], n[0]]),grad(c)))*v_c*ds(4) + (Constant((tan(angle*np.pi/180))**-1.0))*abs(dot(as_vector([-1*n[1], n[0]]),grad(c)))*v_c*ds(5) + (Constant((tan(angle*np.pi/180))**-1.0))*abs(dot(as_vector([-1*n[1], n[0]]),grad(c)))*v_c*ds(6) + (Constant((tan(angle*np.pi/180))**-1.0))*abs(dot(as_vector([-1*n[1], n[0]]),grad(c)))*v_c*ds(7)
# NS equations Divergence form
L2 = 1000.*dot(u, w)*(1./dt0)*dx - 1000.*dot(u0, w)*(1./dt0)*dx + eta*inner(grad(u), grad(w))*dx + eta*inner(grad(u).T, grad(w))*dx - p*inner(I, grad(w))*dx - m_mid*dot(grad(c_mid), w)*dx
L3 = q*div(u)*dx
# Assembling of the system of eqs (4) + (5) + (6) in weak form
L = L0 + L1 + L2 + L3
#----------------------------------------------------------------------
#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Solution of the nonlinear problem (PHASE 0 - for initial solution)
#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, sol, dsol)
# Create nonlinear problem and Newton solver
problem = NonlinearVariationalProblem(L, sol, bcs, J=a)
solver = NonlinearVariationalSolver(problem)
solver.parameters["newton_solver"]["linear_solver"] = "lu"
solver.parameters["newton_solver"]["convergence_criterion"] = "incremental"
solver.parameters["newton_solver"]["relative_tolerance"] = 1.e-10
solver.parameters["newton_solver"]['maximum_iterations'] = 50
# Output file
file0 = File("whole_chip1/potential.pvd", "compressed")
file1 = File("whole_chip1/cells.pvd", "compressed")
file2 = File("whole_chip1/pressure.pvd", "compressed")
file3 = File("whole_chip1/velocity.pvd", "compressed")
# Step in time
t = 0.
i = 1
file0 << (sol.split()[0], t)
file1 << (sol.split()[1], t)
file2 << (sol.split()[2], t)
file3 << (sol.split()[3], t)
while (t < TFIN0):
t += dt0
print('t =', t, 'dt = ', dt0, 'iteration numero', i)
sol0.vector()[:] = sol.vector()
solver.solve()
file0 << (sol.split()[0], t)
file1 << (sol.split()[1], t)
file2 << (sol.split()[2], t)
file3 << (sol.split()[3], t)
i += 1;
#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Solution of the nonlinear problem (PHASE 1 - aspiration)
#----------------------------------------------------------------------
#----------------------------------------------------------------------
v_prof = Expression(("3.5e-6*sin((x[1] - 5.5e-5)*(pi/9.e-05))"), degree=2)
bc5 = DirichletBC(CHS.sub(3).sub(0), v_prof, boundary_markers, 1)
# Collect boundary conditions
bcs = [bc0, bc1, bc2, bc3, bc4, bc5, bc6, bc7]
m0, c0, p0, u0 = sol0.split()
p_0 = Constant(0.)
# Equation (28)
L0 = c*v_m*dx - c0*v_m*dx + dt1*Mw*dot(grad(m_mid), grad(v_m))*dx + dt1*dot(u_mid, grad(c_mid))*v_m*dx
# Update of the boundary condition
L2 = 1000.*dot(u, w)*(1./dt1)*dx - 1000.*dot(u0, w)*(1./dt1)*dx + eta*inner(grad(u), grad(w))*dx + eta*inner(grad(u).T, grad(w))*dx - p*inner(I, grad(w))*dx - m_mid*dot(grad(c_mid), w)*dx + p_0*dot(n, w)*ds(3)
L = L0 + L1 + L2 + L3
# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, sol, dsol)
# Create nonlinear problem and Newton solver
problem = NonlinearVariationalProblem(L, sol, bcs, J=a)
solver = NonlinearVariationalSolver(problem)
solver.parameters["newton_solver"]["linear_solver"] = "lu"
solver.parameters["newton_solver"]["convergence_criterion"] = "incremental"
solver.parameters["newton_solver"]["relative_tolerance"] = 1.e-10
solver.parameters["newton_solver"]['maximum_iterations'] = 50
# Step in time
i = 1
file0 << (sol.split()[0], t)
file1 << (sol.split()[1], t)
file2 << (sol.split()[2], t)
file3 << (sol.split()[3], t)
while (t < TFIN1):
t += dt1
print('t =', t, 'dt = ', dt1, 'iteration numero', i)
sol0.vector()[:] = sol.vector()
solver.solve()
# File outputs each 100 time steps (each 0.001 sec.)
if i%100 == 0:
file0 << (sol.split()[0], t)
file1 << (sol.split()[1], t)
file2 << (sol.split()[2], t)
file3 << (sol.split()[3], t)
# At each time step, store time, aggregate radius and pressure
tps = np.append(tps, (t/10.)); radi = np.append(radi, (Max_height(sol, mesh)*1.E6))
i += 1;
with open('whole_chip1/aspiration_phase1.csv', 'w') as csvfile:
spamwriter = csv.writer(csvfile, delimiter=' ',
quotechar='|', quoting=csv.QUOTE_MINIMAL)
spamwriter.writerow(['tps', 'height'])
for i in range (np.size(tps)):
spamwriter.writerow([tps[i], radi[i]])