FiPy 2D Navier Stokes implementation: problem with derivatives in one direction

Question

I am trying to implement the example from https://nbviewer.jupyter.org/github/barbagroup/CFDPython/blob/master/lessons/14_Step_11.ipynb, but I am facing some problems. I think that the main problem is that I am having some issues with the boundary conditions, as well as defining the terms in the equations.

The PDEs are:

And the initial & boundary conditions are:

I understand that my variables are the two components of the velocity vector (u, v), and the pressure (p). Following the example and using FiPy, I code the PDEs as follows:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import fipy
from fipy import *

# Geometry
Lx = 2   # meters
Ly = 2    # meters
nx = 41  # nodes
ny = 41

# Build the mesh:
mesh = Grid2D(Lx=Lx, Ly = Ly, nx=nx, ny=ny)

# Main variable and initial conditions
vx = CellVariable(name="x-velocity", 
                 mesh=mesh,
                 value=0.)
vy = CellVariable(name="y-velocity",
                 mesh=mesh,
                 value=-0.)
v = CellVariable(name='velocity',
                 mesh=mesh, rank = 1.,
                 value=[vx, vy])
p = CellVariable(name = 'pressure',
                 mesh=mesh,
                 value=0.0)

# Boundary conditions
vx.constrain(0, where=mesh.facesBottom)
vy.constrain(0, where=mesh.facesBottom)
vx.constrain(1, where=mesh.facesTop)
vy.constrain(0, where=mesh.facesTop)
p.grad.dot([0, 1]).constrain(0, where=mesh.facesBottom)  # <---- Can this be implemented like this?
p.constrain(0, where=mesh.facesTop)
p.grad.dot([1, 0]).constrain(0, where=mesh.facesLeft)
p.grad.dot([1, 0]).constrain(0, where=mesh.facesRight)

#Equations
nu = 0.1 #
rho = 1.
F = 0.
# Equation definition
eqvx = (TransientTerm(var = vx) == DiffusionTerm(coeff=nu, var = vx) - ConvectionTerm(coeff=v, var=vx) - ConvectionTerm(coeff= [[(1/rho)], [0]], var=p) + F)
eqvy = (TransientTerm(var = vy) == DiffusionTerm(coeff=nu, var = vy) - ConvectionTerm(coeff=v, var=vy) - ConvectionTerm(coeff= [[0], [(1/rho)]], var=p))
eqp = (DiffusionTerm(coeff=1., var = p) == -1 * rho * (vx.grad.dot([1, 0]) ** 2 + (2 * vx.grad.dot([0, 1]) * vy.grad.dot([1, 0])) + vy.grad.dot([0, 1]) ** 2))
eqv = eqvx & eqvy & eqp

# PDESolver hyperparameters

dt = 1e-3   # (s) It should be lower than 0.9 * dx ** 2 / (2 * D)
steps = 100  #
print('Total time: {} seconds'.format(dt*steps))

# Plotting initial conditions

# Solve
vxevol = [np.array(vx.value.reshape((ny, nx)))]
vyevol = [np.array(vy.value.reshape((ny, nx)))]

for step in range(steps):
  eqv.solve(dt=dt)
  v = CellVariable(name='velocity', mesh=mesh, value=[vx, vy])

  sol1 = np.array(vx.value.reshape((ny, nx)))
  sol2 = np.array(vy.value.reshape((ny, nx)))
  vxevol.append(sol1)
  vyevol.append(sol2)

My result, at time-step 100, is this one (which does not concur with the solution given in the example):

I think that the main issues are in defining the boundary conditions for one specific dimension (e.g. dp/dx = 1, dp/dy = 0), and the derivatives of the variables in one dimension in the equations (in the code, 'eqp').

Can someone enlighten me? Thank you in advance!

Answer 1

Here is what I think is an improved version. At least the result looks more reasonable. The major changes are as follows.

• Using the trick outlined in the linked CFD notebook with the divergence over the time step in the pressure equation.
• Changing the vector velocity, v, to be a face variable so that we can use .divergence directly. Certainly cleans things up, but is a different discretization. I don't know which is more valid.
• Fixing the boundary conditions. I'm not sure p.grad.dot[].constrain was doing anything sensible. Anyway, they aren't needed for a zero gradient as that's the default.
• Not solving all the equations in one matrix. That's best to do once you are confident of solving separately correctly and you have a benchmark to check against.
• The velocity vector variable was being recreated at each step which means that is was having no impact on the equations. v is now explicitly updated in the loop.
• Not using a ConvectionTerm when adding the pressure gradient to the momentum equation. The ConvectionTerm is doing weird weighting and isn't exactly a straightforward difference. In the long run it might be good to use, but not whilst debugging.

Here is the code.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from fipy import (
    CellVariable,
    ConvectionTerm,
    Grid2D,
    FaceVariable,
    TransientTerm,
    DiffusionTerm,
    CentralDifferenceConvectionTerm,
    Viewer
)

from fipy.viewers.matplotlibViewer.matplotlibStreamViewer import MatplotlibStreamViewer

# Geometry
Lx = 2   # meters
Ly = 2    # meters
nx = 41  # nodes
ny = 41

# Build the mesh:
mesh = Grid2D(Lx=Lx, Ly = Ly, nx=nx, ny=ny)

# Main variable and initial conditions
vx = CellVariable(name="x-velocity",
                  mesh=mesh,
                  value=0.,
                  hasOld=True)
vy = CellVariable(name="y-velocity",
                  mesh=mesh,
                  value=-0.,
                  hasOld=True)

v = FaceVariable(name='velocity',
                 mesh=mesh, rank = 1)

p = CellVariable(name = 'pressure',
                 mesh=mesh,
                 value=0.0,
                 hasOld=True)

# Boundary conditions

# top

vx.constrain(1, where=mesh.facesTop)
vy.constrain(0, where=mesh.facesTop)
p.constrain(0, where=mesh.facesTop)

# left

vx.constrain(0, where=mesh.facesLeft)
vy.constrain(0, where=mesh.facesLeft)

# right

vx.constrain(0, where=mesh.facesRight)
vy.constrain(0, where=mesh.facesRight)

# bottom

vx.constrain(0, where=mesh.facesBottom)
vy.constrain(0, where=mesh.facesBottom)

#Equations
nu = 0.1
rho = 1.

dt = 0.01

# Equation definition
eqvx = (TransientTerm() == DiffusionTerm(nu) - ConvectionTerm(v) - (1 / rho) * p.grad[0])
eqvy = (TransientTerm() == DiffusionTerm(nu) - ConvectionTerm(v) - (1 / rho) * p.grad[1])
eqp = (DiffusionTerm(1.) == -1 * rho * (v.divergence**2 - v.divergence / dt))


steps = 200
sweeps = 2
print('Total time: {} seconds'.format(dt*steps))
total_time = 0.0

viewer = MatplotlibStreamViewer(v)
pviewer = Viewer(p)
vxviewer = Viewer(vx)
vyviewer = Viewer(vy)

for step in range(steps):
    vx.updateOld()
    vy.updateOld()
    p.updateOld()

    for sweep in range(sweeps):
        res_p = eqp.sweep(var=p, dt=dt)
        res0 = eqvx.sweep(var=vx, dt=dt)
        res1 = eqvy.sweep(var=vy, dt=dt)

        print(f"step: {step}, sweep: {sweep}, res_p: {res_p}, res0: {res0}, res1: {res1}")

        v[0, :] = vx.faceValue
        v[1, :] = vy.faceValue

    viewer.plot()
    pviewer.plot()
    vxviewer.plot()
    vyviewer.plot()

input('end')