# SUPG Stabilization for Compressible Flow?

288

views

1

I am trying to implement the advection-diffusion equation:

$\frac{dc}{dt}=\nabla\cdot D\nabla c-R\nabla\cdot\left(\vec{v}c\right)-kc$

with velocity specified as follows:

$\vec{v}=-K\cdot\nabla p$→

$\nabla K\cdot\nabla p=-L_p\cdot p$∇

My domain is a square with a circular hole. The hole serves as a source of pressure (Neumann) and metabolyte (Dirichlet). I have implemented the solution with SUPG stabilization and implicit time-stepping. However, my solutions are never stable, and always have oscilations. The solution should be symmetric and smooth. I think the problem is that SUPG is designed for stabilizing incompressible flow, whereas the pressure in this particular case is compressible. Are there stabilization methods for compressible advection diffusion, and are there examples of how to implement it? I have also attached a working example, in case others would like to further investigate.

$\frac{dc}{dt}=\nabla\cdot D\nabla c-R\nabla\cdot\left(\vec{v}c\right)-kc$

`d``c``d``t`=∇·`D`∇`c`−`R`∇·(→`v``c`)−`k``c`with velocity specified as follows:

$\vec{v}=-K\cdot\nabla p$→

`v`=−`K`·∇`p`$\nabla K\cdot\nabla p=-L_p\cdot p$∇

`K`·∇`p`=−`L`_{p}·`p`My domain is a square with a circular hole. The hole serves as a source of pressure (Neumann) and metabolyte (Dirichlet). I have implemented the solution with SUPG stabilization and implicit time-stepping. However, my solutions are never stable, and always have oscilations. The solution should be symmetric and smooth. I think the problem is that SUPG is designed for stabilizing incompressible flow, whereas the pressure in this particular case is compressible. Are there stabilization methods for compressible advection diffusion, and are there examples of how to implement it? I have also attached a working example, in case others would like to further investigate.

```
from fenics import *
from mshr import *
T = 563.0
num_steps = 100
dt = T / num_steps
B = Rectangle(Point(-30/4.5,-30/4.5),Point(30/4.5,30/4.5))
C = Circle(Point(0,0),0.5/4.5,segments=100)
domain = B - C
mesh = generate_mesh(domain, 64)
#Define Constants
D0 = Constant(1.531E-6*60/(4.5*4.5)) #Baseline Diffusivity mm2 min-1
Kappa0 = Constant(0.01*17.0*60.0/(4.5*4.5*4.5*4.5)) #Baseline Permeability, mm4 N-1 min-1
Lp0 = Constant(0.022*60.0/(4.5*4.5)) #Lp*S/V, Permability of vessels to water, mm3 N-1 min-1
R = Constant(1.0) #Retardation factor for velocity, R [0,1]
K_clear0 = Constant(0) #Scaling factor for clearance, liposome/vessel permaility
CathPressure = Constant(3.121E-2*(4.5*4.5*4.5)) #Pressure gradient at the catheter inlet (calculate from velocity), N mm-3
CathConc = Constant(1.0)
mats = FacetFunction('size_t', mesh) #create the object to store BC flags
mats.set_all(0) #set all BCs to intially have 0-flux
class boundary_CathInlet_SD(SubDomain):
def __init__(self):
SubDomain.__init__(self)
def inside(self,x,on_boundary):
tol=1.1
return on_boundary and near(x[0],0,tol) and near(x[1],0,tol)
bcSD = boundary_CathInlet_SD() #Instantiate the boundary subdomain
bcSD.mark(mats,1) #Flag the catheter tip subdomain
ds_p=Measure('ds',domain=mesh,subdomain_data=mats)
n=FacetNormal(mesh)
g=CathPressure
#Pressure Model
Q = FunctionSpace(mesh, 'P', 2)
p = TrialFunction(Q)
q = TestFunction(Q)
PressureVar = -Kappa0*dot(grad(p), grad(q))*dx - Lp0*p*q*dx + g*q*ds_p(1)
a_p, L_p = lhs(PressureVar), rhs(PressureVar)
Ap = assemble(a_p)
Bp= assemble(L_p)
p = Function(Q)
solve(Ap, p.vector(), Bp, 'bicgstab', 'hypre_amg')
vel = -Kappa0*grad(p)
V = FunctionSpace(mesh, 'P', 2)
u = TrialFunction(V)
v = TestFunction(V)
#SUPG Stabilization
h = CellSize(mesh)
nb = sqrt(inner(vel,vel))
Pe = nb*4.5/D0
tau = 0.5*h*pow(4.0/(Pe*h)+2.0*nb,-1.0)
v+=tau*inner(vel,grad(v)) #add the stabilization terms to the test function
IC_map = Constant(0.0)
u_n = project(IC_map, V)
u_n.rename('u', 'initial value')
def boundary_CathInlet(x,on_boundary):
tol=1.1
return on_boundary and near(x[0],0,tol) and near(x[1],0,tol)
inlet_BC = DirichletBC(V, CathConc, boundary_CathInlet)
F = -dt*K_clear0*u*v*dx - D0*dt*dot(grad(u),grad(v))*dx - R*dt*div(vel*u)*v*dx + v*(u_n-u)*dx
a,L = lhs(F),rhs(F)
u=Function(V)
u.rename('u','solution')
t=0.0
for n in range(num_steps):
t+=dt
solve(a == L, u, inlet_BC)
u_n.assign(u)
```

### 1 Answer

2

You've changed the test function $v$ for all terms of the variational form; SUPG specifies that only the advective term should be modified. In addition, the equation you posted is an advection-diffusion-*reaction* equation ( the $kc$ term); you require a different form of stabilization for this type of problem. I have found the following reference useful for studying these types of problems:

https://www.sciencedirect.com/science/article/pii/S0045782597002065

written
6 months ago by
Ryan W

Please login to add an answer/comment or follow this question.