# Postprocessing of CFD results in Mathematica

I was trying to perform computational fluid dynamics (CFD) post-processing in Mathematica after solving primitive variables from Navier-Stokes equation, mainly $(u,v)$, $(Px,Py)$. Since I am not an expert with the finite difference and NDSolve functionality of Mathematica, I used FORTRAN to perform numerical simulation.

The objective is to generate streamlines and vorticity from $u$ and $v$ values in Mathematica.

I looked online and found this link : http://bugman123.com/FluidMotion/index.html where vorticity and streamlines calculations were demonstrated.

I tried cloning the same code using

u = Import["F:\\Lid_Driven_Cavit\\Run_Max_Interations\\u_final.dat", "CSV"];

v = Import["F:\\Lid_Driven_Cavit\\Run_Max_Interations\\v_final.dat","CSV"];

(* u and v represents velocity components in x and y directions
respectively *)

n = Dimensions[u, 1];

dx = 1/(n - 1);

dy = dx;

\[Psi]= \[Integral]u dx= \[Integral]v dy
\[Omega]= \[PartialD]v/dx-\[PartialD]u/dy

psi = Table[0, {n - 1}, {n - 1}]; Do[
psi[[i, j]] =psi[[i, j - 1]] + dx (u[[i + 1, j - 1]] + u[[i + 1, j]])/2,{j, 2, n - 1}, {i, 1, n - 1}];

psi = ListInterpolation[psi, {{0, 1}, {0, 1}}];

ContourPlot[psi[x, y], {x, 0, 1}, {y, 0, 1}, PlotPoints -> 50,
PlotRange -> All, ContourShading -> False,
Contours -> {-0.08, -0.077,-0.07, -0.06, -0.045, -0.025, -0.01,-0.0025, 0,-8*^-6, 2*^-6, 3*^-5, 8*^-5, 1*^-4, 3*^-4, 6*^-4, 9*^-4},
ContourStyle -> Table[{Hue[2 (1 - x)/3]}, {x, 0, 1, 1/16}]]


But this is not working. I am not sure of the error but I am trying to generate the same code as mentioned in the website.

Also for the omega

omega = ListInterpolation[Table[(v[[i, j]] - v[[i - 1, j]])/dx - (u[[i, j]]-  u[[i, j - 1]])/ dx,
{i, 2, n - 1}, {j, 2, n - 1}], {{0, 1}, {0, 1}}];
ContourPlot[omega[x, y], {x, 0, 1}, {y, 0, 1}, PlotPoints -> 50,
PlotRange -> All, ContourShading -> False, Contours -> Range[-14, 7],
ContourStyle -> Table[{Hue[2 (1 - x)/3]}, {x, 0, 1, 1/21}]]


But this is not working as well.

Here is the initial data:

u.dat & v.dat

To be honest, the code you found isn't a good example of coding in Mathematica. I think the following 3 lines are enough for you:

ListStreamPlot[Transpose[{u, v}, {3, 1, 2}], DataRange -> {{0, 1}, {0, 1}}]
curl = Most /@ Differences@v - Most@Differences[u, {0, 1}];
ListContourPlot@LowpassFilter[curl, 1]


You may also want to try ListDensityPlot:

ListDensityPlot[LowpassFilter[curl, 1], ColorFunction -> "DeepSeaColors"]


OK, I forgot about ListLineIntegralConvolutionPlot:

ListLineIntegralConvolutionPlot[Transpose[{u, v}, {3, 1, 2}],
ColorFunction -> "RoseColors", DataRange -> {{0, 1}, {0, 1}}]
(* The following piece of code produces the same graph. *)
(*
{fu, fv} = ListInterpolation[#, {{0, 1}, {0, 1}}] & /@ {u, v}
LineIntegralConvolutionPlot[{fu[x, y], fv[x, y]}, {x, 0, 1}, {y, 0, 1},
ColorFunction -> "RoseColors"]
*)


• Hi, is it possible to use LineIntegralConvolutionplot. The u and v are discrete values and using the function Interpolation, the discrete solution can be converted into a black-box continuous function. Something like : blog.wolfram.com/2013/07/09/… I will be trying to use this method and will post the code if I am able to run it. Commented Aug 20, 2015 at 16:10
• @user11948 Code added, have a look. Commented Aug 21, 2015 at 3:24
• If I wish to evaluate stream-function using trapezoidal rule and perform integration u*dy over the complete domain, is there a mathematica inbuilt function to evaluate it or I need to run a loop. Thanks Commented Sep 9, 2015 at 1:17
• @user11948 Simply use Mean@(Most@# + Rest@#)/2 &@u. You can also try the combination of Interpolation and NIntegrate if you're not insist on trapezoidal rule. Well… I really suggest you to set more effort on learning the basic of Mathematica. Commented Sep 15, 2015 at 10:42

The following part is incorrect, and redundant anyway:

\[Psi]= \[Integral]u dx= \[Integral]v dy
\[Omega]= \[PartialD]v/dx-\[PartialD]u/dy


Otherwise, the code you copied has only one major problem, in the definition of n. If you try to execute it, you get Do::iterb: Iterator {j,2,n-1} does not have appropriate bounds. >> which suggests where the problem might lie. In fact, n is defined as the output of Dimensions, so it is a list and not a number as defined. You can fix that by using:

n = First@Dimensions[u, 1];