Hello everyone,
I am working on modelling arterial blood flow using 1D Navier-Stokes for flow compliant tubes. My full equations are:
\frac{ \partial A }{ \partial t } + \frac{ \partial Au }{ \partial x } = 0
\frac{ \partial u }{ \partial t } + u \frac{ \partial u }{ \partial x } + \frac{ 1 }{ \rho }\frac{ \partial p }{ \partial x } + K_Ru = 0
In my code, they implemented as follows (with some in-between steps) with a \theta time scheme:
# Continuity Equation
AF = A*psi_a*dx \
- An*psi_a*dx \
+ dt * (
theta * ( grad(A * u )[0] * psi_a * dx )
+ (1-theta) * ( grad(An * un)[0] * psi_a * dx )
)
# Momentum Equation
uF = u*psi_u*dx \
- un*psi_u*dx \
+ dt * (
# Advection
theta * alpha * u * grad(u )[0] * psi_u * dx
+ (1-theta) * alpha * un * grad(un)[0] * psi_u * dx
# Area / Pressure
+ theta * ( ( beta/2/rho/Ai/pow(A +DOLFIN_EPS,0.5) ) * grad(A )[0] * psi_u * dx )
+ (1-theta) * ( ( beta/2/rho/Ai/pow(An+DOLFIN_EPS,0.5) ) * grad(An)[0] * psi_u * dx )
# Viscous Dissipation
+ theta * ( ( 2 * (gamma+2) * np.pi * mu * u ) / ( rho * (A +DOLFIN_EPS) ) ) * psi_u * dx
+ (1-theta) * ( ( 2 * (gamma+2) * np.pi * mu * un ) / ( rho * (An+DOLFIN_EPS) ) ) * psi_u * dx
)
I have the above equations working with a single Dirichlet BC for a time-varying inlet (the rest being Neumann conditions); however, when trying to apply any sort of second Dirichlet condition at the outlet my solution diverges. The error is as following:
Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
From what I have read, I should be able to apply one condition at the inlet and one at the outlet for sub-critical problems. Currently, I am trying to apply a half-sin wave for the inlet,
def sinusoidal_wave(t):
return abs( amp * np.sin(n * np.pi * t / T) * np.heaviside(T/n-t, 1) )
and an outlet condition controlling the backward propagating wave at the outlet:
u(x=L) = W_f(1-R_t) - u(x=L-dx)
where,
W_f = u(x=L-dx) - u_0 + 4 ( c - c_0 )
is the forward propagating wave. R_t is the reflection coefficient. c is the wave speed.
To apply them, I create the BCs as
bcs = [
DirichletBC(V.sub(0), inlet, inlet_bdry),
DirichletBC(V.sub(1), outlet, outlet_bdry)
]
and I edit the inlet/outlet variables within the time loop using assign:
inlet.assign(sinusoidal_wave(t))
outlet.assign(Ur) # Ur == u(x=L) from the above eqn
Am I missing something that causes my solution to diverge with the additional Dirichlet BC? Does anyone have any experience with this type of system?
A full example is below with an inlet condition applied area and the above outlet condition applied to the velocity. It runs if the second BC is removed. Apologies for the long example code.
from fenics import *
from tqdm import tqdm
import numpy as np
set_log_level(LogLevel.ERROR)
#%% #* PARAMETERS
#? time parameters
t = 0 # initial time
T = 1 # time period / final time
dt = 3e-05 # time step
Nt = round(T / dt) # number of time steps
#? fluid properties
rho = 1060 # density, kg / m3
mu = 0.004 # viscosity
pext = 0 # external pressure
gamma = 9 # flow profile parameter
alpha = 1 # momentum correction coefficient
#? geometry parameters
L = 1 # domain length, m
Nx = 100 # number of cells
#? initial conditions
R0 = 0.01 # inital radius, m
A0 = np.pi*R0**2 # initial area, m
u0 = 0.0 # initial velocity, m/s
#? vessel properties
v = 0.5 # poisson's ratio, -
E = 400000 # young's modulus, kPa
h = 0.0015 # vessel thickness, m
beta = E * h * np.sqrt(np.pi) / ( 1 - v**2)
#? Flow Properties
c0 = pow(beta/2/rho,0.5) / pow(A0,0.25)
#? inlet parameters
amp = A0*0.1
n = 5
#? outlet parameters
Rt = 0
#? scheme parameters
theta = 0.55
#%% #* FE DOMAIN
mesh = IntervalMesh(Nx, 0, L)
elV = FiniteElement('CG', mesh.ufl_cell(), 2)
V = FunctionSpace(mesh, elV*elV) # mixed space: Area, Velocity
Ai = Expression('A0', degree=2, A0=A0) # initial area
def inlet_bdry(x, on_boundary):
return on_boundary and near(x[0], 0, 1e-14)
def outlet_bdry(x, on_boundary):
return on_boundary and near(x[0], 0, 1e-14)
#%% #* VARIATIONAL PROBLEM
#? New Solution
U = Function(V)
A, u = split(U)
#? Old Solution
Un = Function(V)
An, un = split(Un)
Un.assign(Expression(('A0', 'u0'), degree=2, A0=Ai, u0=u0)) # initial condition
#? Test Functions
psi_a, psi_u = TestFunction(V)
#? Variational Equations : Crank-Nicholson
#* Area Equations
AF = A*psi_a*dx \
- An*psi_a*dx \
+ dt * (
theta * ( grad(A *u )[0]*psi_a*dx )
+ (1-theta) * ( grad(An*un)[0]*psi_a*dx )
)
#* Momentum Equation
uF = u*psi_u*dx \
- un*psi_u*dx \
+ dt * (
#? Momentum Advection
theta * alpha * u * grad(u)[0] * psi_u * dx
+ (1-theta) * alpha * un * grad(un)[0] * psi_u * dx
#? Area Effect Through Pressure
+ theta * ( ( beta/2/rho/Ai/pow(A +DOLFIN_EPS,0.5) ) * grad(A )[0] * psi_u * dx )
+ (1-theta) * ( ( beta/2/rho/Ai/pow(An+DOLFIN_EPS,0.5) ) * grad(An)[0] * psi_u * dx )
#? Viscous Dissipation
+ theta * ( ( 2 * (gamma+2) * np.pi * mu * u ) / ( rho * (A +DOLFIN_EPS) ) ) * psi_u * dx
+ (1-theta) * ( ( 2 * (gamma+2) * np.pi * mu * un ) / ( rho * (An+DOLFIN_EPS) ) ) * psi_u * dx
)
#* Full Variational Form
F = AF + uF
#%% #* BOUNDARY CONDITIONS
def sinusoidal_wave(t):
return abs( amp * np.sin(n * np.pi * t / T) * np.heaviside(T/n-t, 1) + A0 )
inlet = Constant(sinusoidal_wave(t))
outlet = Constant(0)
bcs = [
DirichletBC(V.sub(0), inlet, inlet_bdry),
DirichletBC(V.sub(1), outlet, outlet_bdry) # removing this allows the code to run
]
#%% #* SOLVER
solver_params = {
'newton_solver': {
'linear_solver': 'mumps',
'maximum_iterations': 100,
}
}
for step in tqdm(range(Nt), desc='Processing'):
#? Boundary Conditions
inlet.assign(sinusoidal_wave(t))
Al = Un(L,)[0]
Ul = Un(L,)[1]
c = pow(beta/2/rho/A0,0.5) * pow(Al,0.25)
Wf = Ul - u0 + 4 * ( c - c0 )
Ur = Wf * ( 1 - Rt ) - Ul
outlet.assign(Ur)
#? Solution
J = derivative(F, U)
problem = NonlinearVariationalProblem(F, U, bcs, J=J)
solution = NonlinearVariationalSolver(problem)
solution.parameters.update(solver_params)
solution.solve()
#? Step foward in time
t += dt
Un.assign(U)
Thank you in advance for any insights!