2D Kelvin-Helmholtz instability - velocity diverges

Uncategorized
1

Hi! I’m trying to learn FEniCS and I thought I’d try applying it to the initial conditions from this arxiv article I found on simulating the 2D Kelvin-Helmholtz instability.

I’ve managed to collect a bunch of ad hoc solutions to various problems I’ve run into while setting this up, and I think I’ve got something now. I’m mostly adapting the example of the FEniCS tutorial’s Navier-Stokes solver. Here’s my code:

from sympy import *

x, y = symbols('x[0] x[1]')
U, C, delta = symbols('U_oo c_n delta')
phi = U * exp(-(y-0.5)**2/delta**2) * (cos(8*pi*x) + cos(20 * pi * x))
ux = U * tanh((2 * y - 1)/delta) + C * phi.diff(y)
uy = - C * phi.diff(x)

ux_expression = ccode(ux).replace("M_PI", "pi")
uy_expression = ccode(uy).replace("M_PI", "pi")

from fenics import *
from mshr import *
import numpy as np
from matplotlib import pyplot as plt

T = 0.05            # final time
num_steps = 50   # number of time steps
dt = T / num_steps # time step size
mu = 0.0         # dynamic viscosity
rho = 1            # density
U_oo = 1
c_n = 1e-3
delta = 1/28

# Create mesh
channel = Rectangle(Point(0, 0), Point(1, 1))
domain = channel
mesh = generate_mesh(domain, 64)

class PeriodicBoundary(SubDomain):
    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        return bool(x[0] < DOLFIN_EPS and x[0] > -DOLFIN_EPS and on_boundary)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
        y[0] = x[0] - 1.0
        y[1] = x[1]

pbc = PeriodicBoundary()        
        
# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2, constrained_domain=pbc)
Q = FunctionSpace(mesh, 'P', 1, constrained_domain=pbc)

# Define boundaries
inflow   = 'near(x[0], 0)'
outflow  = 'near(x[0], 1)'
walls    = 'near(x[1], 0) || near(x[1], 1)'
# Define inflow profile
inflow_profile = ('4.0*1.5*x[1]*(1 - x[1]) / pow(1, 2)', '0')
initial_condition_profile = (
    ux_expression, uy_expression,
)

initial_condition = Expression(initial_condition_profile,
                               degree=2,
                               U_oo = U_oo,
                               c_n = c_n,
                               delta = delta,
                              )

# Define boundary conditions
bcu_free_slip = DirichletBC(V.sub(1), 0.0, walls)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [
       bcu_free_slip,
]
bcp = [bcp_outflow]

# 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)
u_n.interpolate(initial_condition)
u_.interpolate(initial_condition)
p_n = Function(Q)
p_  = Function(Q)
u_n

# Define expressions used in variational forms
U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
f  = Constant((0, 0))                # is this gravity?
k  = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
    # time stepping
    # convective term
    # stress tensor - diffusion and pressure
    # pressure on boundary
    # diffusion
    # gravity????
F1 = rho*dot((u - u_n) / k, v)*dx \
   + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds \
   - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f, 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 - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Create XDMF files for visualization output
# xdmffile_u = XDMFFile('navier_stokes_kh/velocity.xdmf')
# xdmffile_p = XDMFFile('navier_stokes_kh/pressure.xdmf')

# Time-stepping
t = 0
from tqdm import auto
with auto.tqdm(total=num_steps) as pbar:
    for n in range(num_steps):

        # Update current time
        t += dt

        # Step 1: Tentative velocity step
        b1 = assemble(L1)
        [bc.apply(b1) for bc in bcu]
        solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

        # Step 2: Pressure correction step
        b2 = assemble(L2)
        [bc.apply(b2) for bc in bcp]
        solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

        # Step 3: Velocity correction step
        b3 = assemble(L3)
        solve(A3, u_.vector(), b3, 'cg', 'sor')

        # Plot solution
        plt.figure()
        mappable = plot(u_**2, title='Velocity squared')
        plt.colorbar(mappable)
#         plt.figure()
#         plot(p_, title='Pressure')

#         # Save solution to file (XDMF/HDF5)
#         xdmffile_u.write(u_, t)
#         xdmffile_p.write(p_, t)

        # Update previous solution
        u_n.assign(u_)
        p_n.assign(p_)

        pbar.set_postfix({'u max': u_.vector().max(), "t": t})
        pbar.update(1)

    
# xdmffile_u.close()
# xdmffile_p.close()

However, on running this within a jupyter notebook, my velocity blows up to infinity somewhere on the mesh - the plot doesn’t display it. I have tried bumping up the number of steps to 500, but that doesn’t seem to have helped.

Is there anything obvious that I’m missing here?

2

I’m having the same problem.

3

Hi, I run it with a uniform mesh and some dissipation (mu=0.0001) and it is working ok.
But the boundary conditions for pressure are giving something non-periodic.
If you run for a long time (T=2) you will see 2 instabilities forming one on the center and the other on the boundary, that one getting destroyed very quickly because the boundary.
I guess one should give bc for the pressure on the top-botton part.