16
Nov
06

Driven Cavity

Here is the flow inside a square box where the flow across the top of the box is moving to the right at Reynolds number 1000. This program uses the finite volume method to solve the Navier-Stokes equations assuming steady, incompressible, viscous, laminar flow. This was calculated on a 300×300 non-uniform grid and took 3 hours to run on my laptop. The animation shows the motion along the streamlines/pathlines.

Driven Cavity – calculated in Fortran 90, plotted in Mathematica 4.2, 11/2/05

The following Mathematica code is not completely accurate at the boundaries, but it gives the basic idea:

(* runtime: 20 minutes *)
n = 50; dx = 1.0/(n - 1); RE = 1000; relax = 0.4; residu = residv = residp = 1;
u = Table[If[j == 1 || i == 2 || i == n, 0, 1], {i, 1, n}, {j, 1, n}];
v = p = pstar = ap = an = as = aw = ae = sc = sp = du = dv = Table[0, {n}, {n}];
lisolv[i0_, j0_, phi0_] := Module[{phi = phi0, p = q =Table[0, {n}]}, Do[q[[j0 - 1]] = phi[[i, j0 - 1]]; Do[p[[j]] = an[[i, j]]/(ap[[i, j]] -as[[i, j]]p[[j - 1]]); q[[j]] = (as[[i, j]]q[[j - 1]] + ae[[i, j]] phi[[i + 1, j]] + aw[[i, j]]phi[[i - 1, j]] + sc[[i, j]])/(ap[[i, j]] - as[[i, j]]p[[j - 1]]), {j, j0, n - 1}]; Do[phi[[i,j]] = p[[j]]phi[[i, j + 1]] + q[[j]], {j, n - 1, j0, -1}], {i, i0, n - 1}]; phi];
calc[i0_, j0_] := Module[{}, {cn, cs, ce, cw} = 0.5 {v[[i,j + 1]] + v[[i0, j0 + 1]], v[[i, j]] + v[[i0, j0]], u[[i + 1, j]] + u[[i0 + 1, j0]], u[[i, j]] + u[[i0, j0]]}dx; cp = Max[0, cn - cs + ce - cw]; {an[[i, j]], as[[i, j]], ae[[i, j]], aw[[i, j]]} = Map[Max[0, 1 - 0.5 Abs[#]] &, RE{cn, cs, ce, cw}]/RE + Map[Max[0, #] &, {-cn, cs, -ce, cw}]; sp[[i, j]] = -cp; ap[[i, j]] = an[[i, j]] + as[[i, j]] + ae[[i, j]] + aw[[i, j]] - sp[[i, j]]];
While[Max[residu, residv, residp] > 0.001, residu = residv = residp = 0; Do[calc[i - 1, j]; sc[[i, j]] = cp u[[i, j]] + dx(p[[i - 1, j]] - p[[i, j]]); du[[i, j]] = relax dx/ap[[i, j]]; residu += Abs[an[[i, j]]u[[i, j + 1]] + as[[i, j]]u[[i, j - 1]] + ae[[i, j]]u[[i + 1, j]] + aw[[i, j]]u[[i - 1, j]] - ap[[i, j]]u[[i, j]] + sc[[i, j]]], {i, 3, n - 1}, {j, 2, n - 1}]; ap /= relax; sc += (1 - relax) ap u; u = lisolv[3, 2, u]; Do[calc[i, j - 1]; sc[[i, j]] = cp v[[i, j]] + dx(p[[i, j - 1]] - p[[i, j]]); dv[[i, j]] = relax dx/ap[[i, j]];residv += Abs[an[[i, j]]v[[i, j + 1]] + as[[i, j]]v[[i, j - 1]] + ae[[i, j]]v[[i + 1, j]] + aw[[i, j]]v[[i - 1, j]] - ap[[i, j]]v[[i, j]] + sc[[i, j]]], {i, 2, n - 1}, {j, 3, n - 1}]; ap /= relax; sc += (1 - relax)ap v; v = lisolv[2, 3, v];Do[an[[i, j]] = dx dv[[i, j + 1]]; as[[i, j]] = dx dv[[i, j]]; ae[[i, j]] = dx du[[i + 1, j]]; aw[[i, j]] = dx du[[i, j]]; sp[[i, j]] = 0; sc[[i, j]] = -((v[[i, j + 1]] - v[[i, j]])dx + (u[[i + 1, j]] - u[[i, j]]) dx); residp += Abs[sc[[i, j]]]; ap[[i, j]] = an[[i, j]] + as[[i, j]] + ae[[i, j]] + aw[[i, j]] - sp[[i, j]]; pstar[[i, j]] = 0, {i, 2, n - 1}, {j, 2, n - 1}]; Do[pstar = lisolv[2, 2, pstar], {2}]; Do[If[i != 2, u[[i, j]] += du[[i, j]](pstar[[i - 1, j]] - pstar[[i, j]])]; If[j != 2, v[[i, j]] += dv[[i, j]](pstar[[i, j - 1]] - pstar[[i, j]])]; p[[i, j]] += 0.3(pstar[[i, j]] - pstar[[n - 1, n - 1]]), {i, 2, n - 1}, {j, 2, n - 1}]];

Here is some Mathematica code to plot streamlines. The blue lines are rotating clockwise and the red lines are rotating counter-clockwise. The plot on the left agrees fairly well with Ercan Erturk’s results.
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}]]


Here is some Mathematica code to plot vorticity contours. The plot on the left agrees fairly well with Ghia’s results.
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}]]

Computational Fluid Dynamics (CFD) Links

Driven CavityMathematica 4.2, 11/16/06


Here is the same driven cavity using the finite difference method. This code uses the Alternating Direction Implicit (ADI) method to solve the vorticity transport equation assuming unsteady, incompressible, viscous flow.
(* runtime: 1 minute *)
<< LinearAlgebra`Tridiagonal`; n = 20; RE = 1000; dx = 1.0/(n - 1); dt = 0.25; cxy = dt/(4dx); sxy = dt/(2RE dx^2); lambda = 1.5; free = Range[2, n - 1];
omega = v = psi = Table[0, {n}, {n}]; u = Table[If[j == n, 1, 0], {n}, {j, 1, n}];
Do[omega[[All, 1]] = 2psi[[All, 2]]/dx^2; omega[[All, n]] = 2 (psi[[All, n - 1]] + dx)/dx^2; omega[[1, All]] = 2 psi[[2, All]]/dx^2; omega[[n, All]] = 2 psi[[n - 1, All]]/dx^2; omega[[free, free]] = Transpose[Partition[TridiagonalSolve[Drop[Flatten[Table[If[i == n - 1, 0, -(cxy u[[i, j]] + sxy)], {j, 2, n - 1}, {i, 2, n - 1}]], -1], Flatten[Table[1 + 2sxy, {j, 2, n - 1}, {i, 2, n - 1}]], Drop[Flatten[Table[If[i == n - 1,0, cxy u[[i + 1, j]] - sxy], {j, 2,n - 1}, {i, 2, n - 1}]], -1], Flatten[Table[(cxy v[[i, j - 1]] + sxy)omega[[i, j - 1]] + (1 - 2sxy)omega[[i, j]] + (-cxy v[[i, j + 1]] + sxy)omega[[i, j + 1]] + If[i == 2, (cxy u[[i - 1, j]] + sxy)omega[[i - 1, j]],0] - If[i == n - 1, (cxy u[[i + 1, j]] - sxy) omega[[i + 1, j]], 0], {j, 2, n - 1}, {i, 2, n - 1}]]], n - 2]]; omega[[free, free]] = Partition[TridiagonalSolve[Drop[Flatten[Table[If[j == n - 1, 0, -(cxy v[[i, j]] + sxy)], {i, 2, n - 1}, {j, 2, n - 1}]], -1], Flatten[Table[1 + 2sxy, {i, 2, n - 1}, {j, 2, n - 1}]], Drop[Flatten[Table[If[j == n - 1, 0, cxy v[[i, j + 1]] - sxy], {i, 2, n - 1}, {j, 2, n - 1}]], -1], Flatten[Table[(cxy u[[i - 1, j]] + sxy)omega[[i - 1, j]] + (1 - 2sxy)omega[[i, j]] + (-cxy u[[i + 1, j]] +sxy)omega[[i + 1, j]] + If[j == 2, (cxy v[[i, j - 1]] + sxy)omega[[i, j - 1]], 0] - If[j == n - 1, (cxy v[[i, j + 1]] - sxy) omega[[i, j + 1]], 0], {i, 2, n - 1}, {j, 2, n - 1}]]], n - 2]; resid = 1; While[resid > 0.001, resid = 0; Do[resid += Abs[psi[[i, j]] - (psi[[i,j]] = lambda(psi[[i, j - 1]] + psi[[i - 1, j]] + psi[[i + 1, j]] + psi[[i, j + 1]] - dx^2omega[[i, j]])/4 + (1 - lambda)psi[[i, j]])], {i, 2, n - 1}, {j, 2, n - 1}]]; Do[u[[i, j]] = (psi[[i, j + 1]] - psi[[i, j - 1]])/(2dx); v[[i, j]] = -(psi[[i + 1, j]] - psi[[i - 1, j]])/(2dx), {i, 2, n - 1}, {j, 2, n - 1}]; ListContourPlot[Transpose[psi], PlotRange -> All, ContourShading -> False], {100}];

Link: Driven Cavity – simple Fortran program by Zheming Zheng that uses this method

Advertisements

0 Responses to “Driven Cavity”



  1. Leave a Comment

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s


Welcome !

You will find here some of my favorite hobbies and interests, especially science and art.

I hope you enjoy it!

Subscribe to the RSS feed to stay informed when I publish something new here.

I would love to hear from you! Please feel free to send me an email : bugman123-at-gmail-dot-com

Archives

Blog Stats

  • 548,971 hits

%d bloggers like this: