SUPG Stabilization for Compressible Flow?

4 months ago by
Ryan W  
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$dcdt =·DcR·(vc)kc  

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=Lp·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):
    def inside(self,x,on_boundary):
        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



#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):
    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)



for n in range(num_steps):
	solve(a == L, u, inlet_BC)

Community: FEniCS Project

1 Answer

4 months ago by

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:

Thank you for the reference, and the clarification. I've been focused on the k=0 case lately just to get things up and running, so I didn't even consider that adding a reaction term might play a role in destabilizing the solution. I'll look through the article and see if there is a method which will work for this situation.
written 4 months ago by Ryan W  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »