This example was originally provided by fred2 <at> qnet.com in http://thread.gmane.org/gmane.comp.python.fipy/911

The model problem of a (compressible) gas confined between two infinite plates is a useful test case for a Navier-Stokes solver. For stationary walls maintained at fixed (but different) temperatures, and assuming continuum no-slip boundary conditions, the the Navier-Stokes equations can be solved analytically.

The energy equation simplifies to the (1D) heat equation (i.e., a DiffusionTerm) to be solved for the temperature profile.

from fipy import *

nx = 401
dx = 1. / nx
mesh = Grid1D(nx=nx, dx=dx)
x, = mesh.getCellCenters()
x = CellVariable(mesh=mesh, value=x)

temperature = CellVariable(name="$T / T_L$", mesh=mesh) is the exponent of the "variable hard sphere (VHS)" viscosity/temperature relationship (i.e., )

omega = 0.657
kappa = temperature.getArithmeticFaceValue()**omega

energy_eqn = DiffusionTerm(coeff=kappa)


which is solved with no-slip boundary conditions (strictly applicable only for Knudsen number ) such that

TracMath macro processor has detected an error. Please fix the problem before continuing.

The command:

'/usr/bin/pdflatex -interaction=nonstopmode 39af16190f9cca54e7b2929b39ea6bed9bb0ba6a.tex'
failed with the following output:
"This is pdfTeX, Version 3.1415926-1.40.10 (TeX Live 2009/Debian)\nentering extended mode\n(./39af16190f9cca54e7b2929b39ea6bed9bb0ba6a.tex\nLaTeX2e <2009/09/24>\nBabel <v3.8l> and hyphenation patterns for english, usenglishmax, dumylang, noh\nyphenation, loaded.\n(/usr/share/texmf-texlive/tex/latex/base/article.cls\nDocument Class: article 2007/10/19 v1.4h Standard LaTeX document class\n(/usr/share/texmf-texlive/tex/latex/base/size10.clo))\n(/usr/share/texmf-texlive/tex/latex/base/inputenc.sty\n(/usr/share/texmf-texlive/tex/latex/base/utf8.def\n(/usr/share/texmf-texlive/tex/latex/base/t1enc.dfu)\n(/usr/share/texmf-texlive/tex/latex/base/ot1enc.dfu)\n(/usr/share/texmf-texlive/tex/latex/base/omsenc.dfu)))\n(/usr/share/texmf-texlive/tex/latex/cmap/cmap.sty)\n(/usr/share/texmf/tex/latex/cm-super/type1ec.sty\n(/usr/share/texmf-texlive/tex/latex/base/t1cmr.fd))\n(/usr/share/texmf-texlive/tex/latex/base/fontenc.sty\n(/usr/share/texmf-texlive/tex/latex/base/t1enc.def)<<t1.cmap>>)\n(/usr/share/texmf-texlive/tex/latex/amsmath/amsmath.sty\nFor additional information on amsmath, use the ?' option.\n(/usr/share/texmf-texlive/tex/latex/amsmath/amstext.sty\n(/usr/share/texmf-texlive/tex/latex/amsmath/amsgen.sty))\n(/usr/share/texmf-texlive/tex/latex/amsmath/amsbsy.sty)\n(/usr/share/texmf-texlive/tex/latex/amsmath/amsopn.sty))\n(/usr/share/texmf-texlive/tex/latex/amscls/amsthm.sty)\n(/usr/share/texmf-texlive/tex/latex/amsfonts/amssymb.sty\n(/usr/share/texmf-texlive/tex/latex/amsfonts/amsfonts.sty))\n(/usr/share/texmf-texlive/tex/latex/tools/bm.sty)\n(/usr/share/texmf/tex/latex/preview/preview.sty\n(/usr/share/texmf/tex/latex/preview/prtightpage.def))\nNo file 39af16190f9cca54e7b2929b39ea6bed9bb0ba6a.aux.\nPreview: Fontsize 10pt\nPreview: PDFoutput 1\n<<ot1.cmap>> (/usr/share/texmf-texlive/tex/latex/amsfonts/umsa.fd)\n(/usr/share/texmf-texlive/tex/latex/amsfonts/umsb.fd)\nPreview: Tightpage -32891 -32891 32891 32891\n[1{/var/lib/texmf/fonts/map/pdftex/updmap/pdftex.map}]\n(./39af16190f9cca54e7b2929b39ea6bed9bb0ba6a.aux) )\n!pdfTeX error: /usr/bin/pdflatex (file ecrm0700): Font ecrm0700 at 600 not foun\nd\n ==> Fatal error occurred, no output PDF file produced!\n"
"\nkpathsea: Running mktexpk --mfmode / --bdpi 600 --mag 1+0/600 --dpi 600 ecrm0700\nmkdir: cannot create directory ././.texmf-var': Permission denied\nmktexpk: /usr/share/texmf/web2c/mktexdir /.texmf-var/fonts/pk/ljfour/jknappen/ec failed.\nkpathsea: Appending font creation commands to missfont.log.\n"

with a wall temperature ratio .

chi = 3.72
BCs = (FixedValue(faces=mesh.getFacesLeft(), value=1.),
FixedValue(faces=mesh.getFacesRight(), value=chi))


The normal momentum equation simplifies to from which the density profile can be obtained.

rho_dimensional = 1. / temperature
rho = rho_dimensional / rho_dimensional[nx / 2]
rho.name = r"$\rho / \rho_0$"


The exact steady-state temperature profile is where and is the plate separation distance.

temperature_analytical = ((chi**(omega+1) - 1.)*x + 1.)**(1./(omega+1))
temperature_analytical.name = temperature.name + " analytical"

rho_dimensional_analytical = 1. / temperature_analytical
rho_analytical = rho_dimensional_analytical / rho_dimensional_analytical[nx / 2]
rho_analytical.name = rho.name + " analytical"

viewer = Viewer(vars=(temperature, temperature_analytical, rho, rho_analytical),
title=r'''gas between heated plates (Navier-Stokes, no-slip)
$T_{left}/T_{right} = %g$, $\omega = %4.3f$''' %  (chi, omega))


Starting from a linear profile, it takes only a few sweeps to reach convergence.

temperature.setValue(1 + (chi - 1.)*x)

res = 1e10
while res > 1.e-6:
res = energy_eqn.sweep(var=temperature,
boundaryConditions=BCs)

viewer.plot()

print temperature.allclose(temperature_analytical, atol=1.e-4)
print rho.allclose(rho_analytical, atol=1.e-4)


Reference Wadsworth:1993 compares finite difference Navier-Stokes and particle simulation method predictions in the rarefied flow (i.e., slip) regime for which some experimental data are available.

@article{Wadsworth:1993,