Nonlinear convection diffusion equation

Question

Hi, I am fairly new to fenics, and apologize if my question is too basic. I would like to solve a nonlinear convection-diffusion equation with robin b.c.s. Here is the governing equation

$\frac{\partial c_i}{\partial t} + \tilde{u}\frac{\partial c_i}{\partial x} + \frac{\partial}{\partial z}(c_i(1-c_i)) = \frac{1}{Pe}\frac{\partial^2c_i}{\partial z^2}$ with robin b.c.s at $z=0$ and $z=1$ as follows: $\frac{\partial}{\partial z}(c_i(1-c_i)) = \frac{1}{Pe}\frac{\partial^2c_i}{\partial z^2}$. I have re-written the robin b.c.s as follows: $c_i(1-c_i) = \frac{1}{Pe}\frac{\partial c_i}{\partial z}$. I have attached the complete and minimal code below.

Basically, this equation has two convection term. The first one is the material moves with velocity (linear velocity presented below), the second convection term makes the material goes upward. So material enters with concentration with 0.5, should end up with something greater than 0.5 in above (because they move upward) and less than 0.5 at below. The solution I get from following does not show concentration less than 0.5 at below. But it does have the "correct" the shape of concentration contour.

I think maybe my problem is that

1. the robin boundary condition is only on z direction, which I don't think my code correctly reflect that. ( corresponding code: u(1.0-u)vdsTop - u(1.0-u)vdsBottom)

2. This might not be related to the problem I see at here, but more toward my future work. The second convection term can be a bit more complicated than $\frac{\partial}{\partial z}(c_i(1-c_i))$, and can actually have $\frac{\partial}{\partial z}(g(z)c_i(1-c_i))$. In this case, I can not take $g(z)$ out of derivative, and do dot product. In this case, should I manually take divergence on this expression in fenics?

   from dolfin import *
   # Create mesh and build function space
   mesh = UnitSquareMesh(200, 200, 'crossed')
   V = FunctionSpace(mesh, "Lagrange", 1)
   # Create boundary markers
   boundary_parts = FacetFunction('size_t', mesh)
   left   = AutoSubDomain(lambda x: near(x[0], 0.0))
   top  = AutoSubDomain(lambda x: near(x[1], 1.0))
   bottom = AutoSubDomain(lambda x: near(x[1], 0.0))
   left  .mark(boundary_parts, 1)
   top   .mark(boundary_parts, 2)
   bottom.mark(boundary_parts, 3)
   
   #  initial condition
   ic = Expression("0.5", degree=3)
   
   # Equation coefficients
   Pe = Constant(100.) # Pe number
   velocity = Expression(("x[1]", "0"), degree=3) # convecting velocity
   local_shear_rate = Expression(("0", "1"), degree=3)  #local shear rate
   
   # Define boundary measure on Neumann part of boundary
   dsTop = Measure("ds", subdomain_id=2, subdomain_data=boundary_parts)
   dsBottom = Measure("ds", subdomain_id=3, subdomain_data=boundary_parts)
   
   def separation_v(u):
      return local_shear_rate*(1.0 - u)
   
   # Define steady part of the equation
   def operator(u, v):
    return  1.0/Pe*inner(u.dx(1), v.dx(1))*dx  + dot(velocity, grad(u))*v*dx + dot(local_shear_rate, grad(u*(1.0-u)))*v*dx - u*(1.0-u)*v*dsTop - u*(1.0-u)*v*dsBottom
   
   # Define trial and test function and solution at previous time-step
     u = Function(V) # set as nonlinear term
     v = TestFunction(V)
     u0 = Function(V)
   
   # Time-stepping parameters
      t_end = 1
      dt = 0.05
      theta = Constant(0.5) # Crank-Nicolson scheme
   
    # Define time discretized equation
       F = (1.0/dt)*inner(u-u0, v)*dx + theta*operator(u, v) + (1.0-theta)*operator(u0, v)
   
     # Define boundary condition
        bc = DirichletBC(V, Constant(0.5), boundary_parts, 1)
   
     # Prepare solution function and solver
         solve(F == 0, u, bc)
   
     # Prepare initial condition
          u0.interpolate(ic)
   
     # Create file for storing results
          vtkfile = File('mass_transfer/mass_transfer.pvd')
   
     # Time-stepping
          t = 0.0
          u.rename("u", "concentraiton")
          u.interpolate(ic)
   
     # save initial solution
          vtkfile << u
   
     while t < t_end:
   
    # Solve the problem
      solve(F == 0, u, bc)
   
     # Store solution to file and plot
       vtkfile << (u, t)
       plot(u, title='Solution at t = %g' % t)
   
     # Move to next time step
        u0.assign(u)
         t += dt
   

Answer 1

Dear dengzhekai,

few comments:

1. It is very skeptical, to use degree=3 for constant Expression :)
2. I do not understand your strange use of initial conditions.
3. Your intuition is right, condition - u*(1.0-u)*v*dsTop - u*(1.0-u)*v*dsBottom is not good. The normals are oriented in opposite way on top and bottom boundary! You should change it to - u*(1.0-u)*v*dsTop + u*(1.0-u)*v*dsBottom.

I have changed your code with the above adjustments to:

from dolfin import *
# Create mesh and build function space
mesh = UnitSquareMesh(50, 50, 'crossed')
V = FunctionSpace(mesh, "Lagrange", 1)
# Create boundary markers
boundary_parts = FacetFunction('size_t', mesh)
left   = AutoSubDomain(lambda x: near(x[0], 0.0))
top  = AutoSubDomain(lambda x: near(x[1], 1.0))
bottom = AutoSubDomain(lambda x: near(x[1], 0.0))
left  .mark(boundary_parts, 1)
top   .mark(boundary_parts, 2)
bottom.mark(boundary_parts, 3)

#  initial condition
ic = Expression("0.5", degree=0)

# Equation coefficients
Pe = Constant(10.) # Pe number
velocity = Expression(("x[1]", "0"), degree=3) # convecting velocity
local_shear_rate = Expression(("0", "1"), degree=3)  #local shear rate

# Define boundary measure on Neumann part of boundary
dsTop = Measure("ds", subdomain_id=2, subdomain_data=boundary_parts)
dsBottom = Measure("ds", subdomain_id=3, subdomain_data=boundary_parts)

def separation_v(u):
    return local_shear_rate*(1.0 - u)

# Define steady part of the equation
def operator(u, v):
    return  1.0/Pe*inner(u.dx(1), v.dx(1))*dx  + dot(velocity, grad(u))*v*dx + dot(local_shear_rate, grad(u*(1.0-u)))*v*dx - u*(1.0-u)*v*dsTop + u*(1.0-u)*v*dsBottom

# Define trial and test function and solution at previous time-step
u = Function(V) # set as nonlinear term
v = TestFunction(V)
u0 = Function(V)

# Time-stepping parameters
t_end = 1
dt = 0.05
theta = Constant(0.5) # Crank-Nicolson scheme

# Define time discretized equation
F = (1.0/dt)*inner(u-u0, v)*dx + theta*operator(u, v) + (1.0-theta)*operator(u0, v)

# Define boundary condition
bc = DirichletBC(V, Constant(0.5), boundary_parts, 1)

# Prepare initial condition
u0.interpolate(ic)

# Create file for storing results
vtkfile = File('mass_transfer/mass_transfer.pvd')

# Time-stepping
t = 0.0

while t < t_end:

    # Solve the problem
    solve(F == 0, u, bc)

    # Store solution to file and plot
    u.rename("u", "concentration")
    vtkfile << u

    # Move to next time step
    u0.assign(u)
    t += dt

Do you expect (after some time) pictures like this?

Comments

Wow. Thanks for the answer. I did not realize the direction of normal vector is different from top to the bottom. Thanks!

---

The normal vector to the exterior boundary is defined as pointing outwards.

---

Hi mhabera and nate, thank you very much for answering my question. I have a new question about evaluating facetnormal on the internal boundary line, see here. I basically see the same behavior as I posted in this thread. I wonder would it be possible to provide some insight in my new question ?