5.3.3. Rayleigh-Taylor instability
Now, we can make use of the power of pyoomph to easily couple the Navier-Stokes equations with an advection-diffusion equation. Here, we will solve a denser fluid residing atop of a less dense fluid. We apply the Boussinesq approximation, i.e. the mass density \(\rho\) is assumed to be constant in the continuity equation and also in the inertia term, but the bulk force will depend on the varying density. We set the bulk force to \(\vec{f}=-100 c \vec{e}_y\) and solve for the composition field \(c\) that ranges from \(0\) (bottom light fluid) to \(1\) (heavier top fluid). We choose a small diffusivity so that advection can dominate. For the equations for the advection-diffusion of \(c\) and the velocity-pressure pair, we use the predefined equations from pyoomph. The corresponding classes are NavierStokesEquations and AdvectionDiffusionEquations. However, we could also use the equations we have developed in the previous sections instead.
from pyoomph import *
# Use the pre-defined equations for Navier-Stokes and advection-diffusion
from pyoomph.equations.navier_stokes import *
from pyoomph.equations.advection_diffusion import *
class RayleighTaylorProblem(Problem):
def __init__(self):
super(RayleighTaylorProblem,self).__init__()
self.W,self.H=0.25, 1 # Size of the box
self.rho,self.mu=0.01, 1 # density and viscosity
self.Nx=4 # elmenents in x-direction
self.max_refinement_level=4 # max. 4 times refining
def define_problem(self):
# add the mesh
self.add_mesh(RectangularQuadMesh(size=[self.W,self.H],N=[self.Nx,int(self.Nx*self.H/self.W)]))
eqs=MeshFileOutput() # output
bulkforce=100*var("c")*vector(0,-1) # bulkforce: depends on the composition c
# Advection diffusion equation: Advected by the velocity
eqs+=AdvectionDiffusionEquations(fieldnames="c",wind=var("velocity"),diffusivity=0.0001,space="C1")
# Navier-Stokes with the c-dependent bulk formce
eqs+=NavierStokesEquations(bulkforce=bulkforce,dynamic_viscosity=self.mu,mass_density=self.rho)
# Initial condition
xrel,yrel=var("coordinate_x")/self.W,var("coordinate_y")/self.H-0.5
eqs+=InitialCondition(c=tanh(100*(yrel-0.0125*cos(2*pi*xrel))))
# Refinements based on c and velocity
eqs+=SpatialErrorEstimator(c=1,velocity=1)
# Adding no-slip conditions
for wall in ["left","right","top","bottom"]:
eqs+=DirichletBC(velocity_x=0,velocity_y=0) @ wall
# Fix one pressure degree
eqs+=DirichletBC(pressure=0)@"bottom/left"
self.add_equations(eqs@"domain")
if __name__=="__main__":
with RayleighTaylorProblem() as problem:
problem.run(10,numouts=50,spatial_adapt=1)
If you have read the tutorial up to here, you should understand all steps. The idea is to create a NavierStokesEquations object with a body force depending on the variable var("c"). This variable is defined and solved in the AdvectionDiffusionEquations, which in turn get the field var("velocity") for the advection. Thereby, both parts are coupled. Since we allow no normal outflow, we have to fix a single pressure degree, which we do in the lower left corner. The results are depicted in Fig. 5.8.
Fig. 5.8 Miscible Rayleigh-Taylor instability by coupling the Navier-Stokes equations with an advection-diffusion equation.