Postprocessing of CFD results in Mathematica

Question

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

Answer 1

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"]
*)

Answer 2

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];

instead of your definition.

Once you do that, all the rest works as intended: