Taylor-Green vortex sheet, reduced units

The Taylor-Green vortex sheet is a solution to the 2D Navier-Stokes equations for an incompressible Newtonian fluid:

$latex frac{d mathbf{u}}{d t}= nu nabla^2 mathbf{u} – nabla p/rho ,$

where $latex mathbf{u}$ is the velocity field, $latex p$ is the pressure, $latex nu=mu/rho$ is the kinematic viscosity, and $latex rho$ is the fixed density of the fluid. The time derivative is a total derivative:

$latex frac{d mathbf{u}}{d t} =  frac{partial mathbf{u}}{partial t} + (mathbf{u} cdot nabla) mathbf{u}$

It is common to choose parameters that simplify the equations, but that can obscure the role of the different parameters. In the following, I provide expressions with all relevant parameters included, with their physical dimensions. I later pass to dimensionless, or reduced, units, in terms of the Reynolds and Courant numbers.


Sigue leyendo

Quick and dirty Rayleigh-Taylor instability

The Rayleigh-Taylor instability is a well-known benchmarck for CFD codes. The idea is to start with two phases, on on top of the other, the lighter one being underneath. The interface is slightly perturbed, and this plume appears. I describe a quick and dirty way of getting this instability.

Sigue leyendo

Sound waves with attenuation

Just a simple derivation of the role of attenuation in the standard sound wave equation.

Starting with the Navier-Stokes momentum equation

$latex frac{partial }{partial t} mathbf{u} + mathbf{u} nabla mathbf{u} = – frac{1}{rho} nabla p + mu nabla^2 mathbf{u} $

Plus mass continuity
$latex nabla cdot mathbf{rho u} + frac{partial rho}{partial t} =0 $

Now, let us suppose small pressure and density fluctiations about a constant background. I.e:
$latex mathbf{u} approx mathbf{u} qquad papprox p_0+delta p qquad rhoapprox rho_0+delta rho .$

The first equation means the medium onto which sounds travel is not moving (should be modified for e.g. wind). In this case, the fluctuations in pressure and density, being small, should be related by a constant which, with some foresight, we will call $latex c^2$:
$latex delta p approx c^2 delta rho .$

Then, neglecting second and higher order perturbation terms we obtain linearized equations involving only the velocity and the pressure (not the density):

$latex frac{partial }{partial t} mathbf{u} = – frac{1}{rho_0} nabla p + nu nabla^2 mathbf{u} $

$latex rho_0 c^2 nabla cdot mathbf{rho u} + frac{partial p}{partial t} =0 $

Here, we have dropped the “$latex delta$s” (it’s confusing, but the expressions are so much cleaner). The kinematic viscosity is defined as $latex nu=mu/rho_0$.

If we differentiate with respect to time in the first, and with respect to space (applying $latex nabla$) in the second, we can eliminate the pressure term $latex nabla(partial p/partial t)$ (I know, exchanging both derivatives, which makes mathematicians nervous.) The following equation for the velocity results:

$latex frac{partial^2 }{partial t^2 } mathbf{u} =c^2 nabla^2 mathbf{u} + nu frac{partial }{partial t} nabla^2 mathbf{u} .$

Clearly, a wave equation with no viscosity. Let’s try a wave solution of the form

$latex mathbf{u} = mathbf{u}_0 e^{i(alpha t – k x)} $


$latex -alpha^2 = -c^2 k^2 – i nu alpha k^2 $


$latex alpha^2 – 2 i beta alpha – c^2 k^2 =0, $

where we introduce $latex beta=(nu k^2)/2$.

The solution to the equation is

$latex alpha = i beta pm k c sqrt{ 1 – (beta/c k)^2}$

Therefore our solution looks like

$latex mathbf{u} = mathbf{u}_0 e^{ – beta t } e^{i(omega t – k x)} ,$

with an exponential dampening governed by $latex beta$: the relaxation time is $latex tau=1/beta=2/(nu k^2)$, so that higher viscosity dampens waves faster, and high $latex k$ waves (short wavelengths) are dampened faster than low $latex k$ (long wavelength) ones.

The angular frequency is $latex omega= k csqrt{1 – (beta/c k)^2}$, meaning that this waves have dispersion. Long waves have the expected $latex omega= k c $ non-dispersive relationship, and all of them have the same phase velocity regardless of their wavelength, $latex v_p = omega /k = c $, the speed of sound. Shorter waves, on the other hand, will have that phase velocity reduced:

$latex v_p = omega /k = c sqrt{1 – (nu k/ (2c) ) ^2} ,$

until $latex k = 2 c / nu $, beyond which sound waves do not oscilate and the overdamped region begins.

The group velocity, on the other hand, is

$latex v_g = frac{d omega }{ d k} = c frac{1 – frac12 (frac{nu k }{c})^2 ) }{ sqrt{1 – (frac{nu k }{2 c})^2} } ,$

which vanishes before the phase velocity, at $latex k = sqrt{2} c / nu $, and then goes negative (!?).

Sparse Poisson problem in eigen

OwlgenBack to scientific computing. Lately, I have been using the Eigen libraries for linear algebra. They were recommended by the CGAL project, and they indeed share a common philosophy. Thanks to the rather high C++ level they can accomplish this sort of magic:

[code language=”cpp”] int n = 100;

VectorXd x(n), b(n);

SpMat A(n,n); fill_A(A); fill_b(b); // solve Ax = b

Eigen::BiCGSTAB solver;

//ConjugateGradient solver; solver.compute(A); x = solver.solveWithGuess(b,x0); [/code]


Notice that A is a sparse matrix! I am next describing how to use this in order to solve the 1D Poisson equation.

Sigue leyendo

OpenFOAM cheatlist

A quick cheatsheet for OpenFOAM. In italics, things that are useful but not part of OpenFOAM proper. Interesting read: The OpenFOAM Technology Primer


Shortcuts to directories

(type “alias” to reveal these)

  • run (go to own’s running directory)
  • foam
  • foamfv
  • foam3rdParty (hit <tab> for these longish commands!)
  • tut
  • app
  • sol
  • util
  • lib
  • src

Environment variables

  • echo $FOAM_ <tab>  (directories)
  • echo $WM_ <tab>  (building, aka compiling, settings)


  • Generation
  • Import / export
    • foamMeshToFluent
    • fluentMeshToFoam, etc …
  • Operations
    • refineHexMesh
    • transformPoints
    • makeAxialMesh
    • collapseEdges
    • autoPatch
    • mirorMesh
  • Properties
    • checkMesh


  • setFields
  • topoSet
  • patchAverage
  • patchIntegrate
  • vorticity
  • yPlusRAS
  • yPlusLES
  • boxTurb
  • applyBoundaryLayer
  • R
  • wallShearStress


  • foamHelp (e.g. foamHelp boudary -field U)


  • sample
  • paraview
  • probeLocations


  • icoFoam
  • interFoam
  • many, many others


  • foamJob
  • decomposePar
  • reconstructPar
  • mpirun
  • nohup


  • foamNew source App …
  • doxygen doxyfile
  • gdb
  • valgrind
  • wmake
  • wclean
  • aliases for settings: wm32, wm64, wmSP, wmDP, wmSET, wmUNSET

Sorting eivenvectors in octave

A brief note. Doing ee= eig(m); produces a vector with all eigenvalues of matrix m. These are unsorted, so ee=sort(eig(m)); produces a vector with a sorted list of all eigenvalues of matrix m.

If we want the eigenvectors, [vv,eed]=eig(m); produces a diagonal matrix eed whose elements are the same as ee, and are unsorted. Who to produce a sorted vector with the eigenvalues, and re-order the eigenvectors accordingly? Thus:


[ee,perm] = sort(diag(eed));


from gnuplot to python

Part of my moving to python for science.

I have lately been using gnuplot to monitor the progress and end results of my simulations.  These are 2D hydrodynamic simulations, which involve 2D scalar and vector fields.

2D plots

With gnuplot, I used to run stuff like

plot ‘8000/points.dat’ u 1:2

for a 2D plot of a scalar field. On column 1 I have the x values of the coordinates and on column 2 the y values, so this just plots the coordinates. The “8000” is the directory for time step 8000, where the data file is saved.

OK, now running ipython –pylab I have to run

dt = pylab.loadtxt(‘8000/points.dat’)




2D scatter plot with colors


Important things to notice:

  • python starts numbering c-style: at 0! so, despite pylab having been designed in order to mimic matlab syntax, there is a clear departure here.
  • it’s easy to assign a color to the points, in this case column 12 (11 for python) is a density field. The default color scale looks ok (it goes from blue to red).
  • the size is set to 50 for a nice full-screen view.
  • New plots will appear on top of this one. If this is not desired, close the window, or run pylab.clf()

Vector plots

This is quite easy in gnuplot:

plot ‘8000/points.dat’ u 1:2:($5/10):($6/10) w vec

On columns 5 and 6 I store the x and y components of the velocity. In this example they are reduced ten-fold in order to hace a decent plot (and this has to set by hand…).

Now we may run:


which does the same. A bonus is that the vector length is automatically calculated from the average vector length. This may the right solution, if not the option “scale=xx” may be used. The higher xx is, the shorter the arrows (??). The fifth entry is a color code, which does not look to good at present because, I think, they are taken as RGB values, unlike with scatter.


Vector plot -- funny senseless colors

Vector plot — funny senseless colors

 3D plots

This one’s harder. To get a scatter 3D plot of three columns I used to run

splot ‘8000/points.dat’ u 1:2:12

Now, it’s more like

from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()

ax = fig.add_subplot(111, projection=’3d’)



That produces a quite nice scatter plot with tick marks on the three axes.


3D scatter plot



Getting rid of capitalized ribbons in office 2013

A_(capital_and_small)Do you hate capital letters? I do, since FORTRAN. Fortunately, people are learning they should not be used. Specially editors. Nevertheless, here comes office 2013 with brand new capitalized ribbons. Straight capital letters, I mean, not small caps.

Anyway, I found somewhere a trick to make sure menus stay in lower case: append a space before, or after, the ribbon name (right-click on any ribbon to change their names). This name, by the way, is correctly typed (!).