Hello everyone,
for my master thesis I am trying to solve the incompressible Navier-Stokes equations. My problem involves simulating an oscillating fluid using ultra sound. The Navier-Stokes equations are scaled with the parameters beta and epsilon which are derived from the dimensions of an ultra sound wave. You can find the pde I am trying to solve in the beginning of my code.
For spatial discretization I am using the IPCS I have found in the “solving pde in python tutorial”.
I am passing an initial pressure and for my boundary conditions I am using Dirichlet BC at the Inlet and the Outlet for both the velocity and the pressure to create the oscillation in the fluid. For the horizontal walls I am passing slip BC.
My problem compiles but explodes after a couple of time steps. Here is my code:
"""
Incompressible Navier-Stokes equations using the Incremental Pressure Correction
Scheme (IPCS).
beta * u' + beta*epsilon * u . nabla(u)) - div(sigma(u, p)) = 0
div(u) = 0
"""
from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
#####-ultrasound wave parameters
f=1.e6 # frequency of ultrasound taken from JV [Hz]
omega= 2. * np.pi * f # angular frequency of ultrasound calculated [Hz]
c_f=1484. # speed of sound in unperturbated fluid taken from JV[m/s]
tau=1./f # time period of sound taken from JV [1.e-6s]
wavelength=c_f/f # wavelength of sound [1.484e3/1e6=1.484e-3m]
#####-the fluid parameters are for water at room temperature
rho=998. # unpertubated mass density of the fluid taken from JV[kg/m^3]
p_0=101325. # unpertubated pressure of fluid taken from JV [Pa]
nu=mu/rho # kinematic shear viscocity [m^2/s]
R = 1e-6 # characteristic length of colloidal particle taken from JV[m]
#####-ultra sound amplitudes
deltap=10e3 #pressure amplitude of ultrasound taken from JV [Pa]
deltav=deltap/(rho*c_f) #velocity amplitude of ultrasound wave traveling in x direction taken from JV [m/s] and Landau Lifschitz Fluid Mechanics pg 252 eq. (64.11)
a = deltav/omega # characteristic oszillation amplitude of colloidal particle [m] << R
#####-dimensionless parameters
beta = R**2*omega/nu # acoustic Reynolds number =6.258
epsilon = a/R # smallness parameter =U/(w*R)=0.001
p_scale = a*omega*mu/R
v_scale = a*omega
#initial scaled pressure
p_0_ent = p_0/p_scale #=14976
#scaled inlet/outlet
deltav_ent=deltav/v_scale #=1
deltap_ent=deltap/p_scale #=1478
#####-computational parameters
l1, l2= wavelength/R, 200e-6/R #domain size taken from JV =1484, 200
t_max=40*tau #simulation duration
num_steps=4000
dt=t_max/(num_steps)
tol=1E-14
t=0
nx=2*1484
ny=20
deltax=2*l1/nx
# Create XDMF files for visualization output
xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity.xdmf')
xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure.xdmf')
# Create mesh
mesh = RectangleMesh(Point(0,0), Point(2*l1, l2), nx, ny)
#define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)
# Define boundaries
inflow = 'near(x[0], 0)'
outflow = 'near(x[0], 2*1.484e-3/1e-6)'
walls = 'near(x[1], 0) || near(x[1], 200e-6/1e-6)'
# Define inital pressure and initial velocity
u_ini = Expression(('v_x_0', 'v_y_0'), v_x_0=Constant(0.0), v_y_0=Constant(0.0), degree=1)
u_n = interpolate(u_ini, V)
p_ini = Expression('p_0_ent', p_0_ent=Constant(p_0_ent), degree=1)
p_n = interpolate(p_ini, Q)
# Define inflow profile
inflow_profile_p=Expression('p_0+deltap_ent*sin(2*pi*f*t)', p_0=Constant(p_0_ent), pi=Constant(np.pi), f=Constant(f), t=t, deltap_ent=Constant(deltap_ent), degree=1)
inflow_profile_u=Expression(('vx_0+deltav_ent*sin(2*pi*f*t)', 'vy_0'), vx_0=Constant(0), vy_0=Constant(0), pi=Constant(np.pi), f=Constant(f), t=t, deltav_ent=Constant(deltav_ent), degree=1)
# Define boundary conditions
bcu_inflow = DirichletBC(V, inflow_profile_u, inflow)
bcp_inflow=DirichletBC(Q, inflow_profile_p, inflow)
bcu_outflow = DirichletBC(V, inflow_profile_u, outflow)
bcp_outflow=DirichletBC(Q, inflow_profile_p, outflow)
bcu = [bcu_inflow,bcu_outflow ]
bcp = [bcp_inflow, bcp_outflow]
xdmffile_u.write(u_n, t)
xdmffile_p.write(p_n, t)
# Define expressions used in variational forms
U = 0.5*(u_n + u)
n = FacetNormal(mesh)
k = Constant(dt)
beta = Constant(beta)
epsilon = Constant(epsilon)
param = Constant(beta*epsilon)
# Define strain-rate tensor
def epsilon(u):
return sym(nabla_grad(u))
# Define stress tensor
def sigma(u, p):
return 2*epsilon(u) - p*Identity(len(u))
# Define variational problem for step 1
F1 = beta* dot((u - u_n) / k, v)*dx + param * dot(dot(u_n, nabla_grad(u_n)), v)*dx + inner(sigma(U, p_n), epsilon(v))*dx
a1 = lhs(F1)
L1 = rhs(F1)
# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (beta/k)*div(u_)*q*dx
# Define variational problem for step 3
a3 = beta*dot(u, v)*dx
L3 = beta*dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx
# create progress bar
progress = Progress('Time-stepping')
set_log_level(16)
t = 0
#####---time loop and solving of the variational problem
for i in range(num_steps):
#####-update current time
t+= dt
#####-update boundary condition
inflow_profile_u.t=t
inflow_profile_p.t=t
bcu_inflow = DirichletBC(V, inflow_profile_u, inflow)
bcp_inflow=DirichletBC(Q, inflow_profile_p, inflow)
bcu_outflow = DirichletBC(V, inflow_profile_u, outflow)
bcp_outflow=DirichletBC(Q, inflow_profile_p, outflow)
bcu = [bcu_inflow,bcu_outflow ]
bcp = [bcp_inflow, bcp_outflow]
#solve the first step
solve(a1==L1, u_, bcu, solver_parameters={'linear_solver' : 'bicgstab', 'preconditioner':'hypre_amg'})
solve(a2==L2, p_, bcp, solver_parameters={'linear_solver' : 'bicgstab', 'preconditioner':'hypre_amg'})
solve(a3==L3, u_, bcu, solver_parameters={'linear_solver':'cg', 'preconditioner':'sor'})
#save solution to vtk file every 100th time step
#if (i % 100) ==0:
xdmffile_u.write(u_, t)
xdmffile_p.write(p_, t)
#update previous solution
u_n.assign(u_)
p_n.assign(p_)
print(t)
I can’t really get to the bottom of this and I am wondering if the problem lies in the scaling as in the tutorial the non scaled incompressible NSE are solved or somewhere else my implementation.
Any help would be greatly appreciated! Thanks!