## Archive for the 'potential flow' Category

17
Nov
07

### Periodic Steady Vortices in a Stagnation Point Flow

 Periodic Steady Vortices in a Stagnation Point Flow. I got the idea for this potential flow animation from this paper. The entrained fluid was integrated using the 4th order Runge-Kutta method.
11
Nov
07

### Schwarz-Christoffel Flow

Here’s a potential flow of a vortex near some stair steps. This was accomplished by mapping a flat wall to stair steps using a Schwarz-Christoffel transformation:

```(* runtime: 0.7 second *) Clear[z]; gamma = 1; z0 = 1 + I; z1 = I;z2 = 0; z3 = 1; zeta1 = -1; zeta2 = 0; zeta3 = 1; theta1 = 3Pi/2; theta2 = Pi/2; theta3 = 3Pi/2; z = z[zeta] /. DSolve[{z'[zeta] == K(zeta - zeta1)^(theta1/Pi - 1)(zeta - zeta2)^(theta2/Pi - 1)(zeta -zeta3)^(theta3/Pi - 1), z[zeta1] == z1}, z[zeta], zeta][[1]]; z = z /. K -> (K /. Solve[(z /. zeta -> zeta2) == z2, K][[1]]); zeta0 = zeta /. FindRoot[z == z0, {zeta, zeta2 + I}]; Z[w_] = z /. zeta -> (zeta /. Solve[w == (I gamma/(2Pi))(Log[zeta - zeta0] - Log[zeta - Conjugate[zeta0]]), zeta][[1]]); ToPoint[z_] := {Re[z], Im[z]}; grid = Table[ToPoint[Z[phi + I psi]], {phi, -1, 0, 0.05}, {psi, -1*^-6, -1, -0.02}]; Show[Graphics[{Map[Line, grid], Map[Line, Transpose[grid]]}, AspectRatio -> Automatic, PlotRange -> {{-1, 2}, {-1, 2}}]]```

08
Nov
07

### Joukowski Airfoil

These animations were created using a conformal mapping technique called the Joukowski Transformation. A Joukowski airfoil can be thought of as a modified Rankine oval. It assumes inviscid incompressible potential flow (irrotational). Potential flow can account for lift on the airfoil but it cannot account for drag because it does not account for the viscous boundary layer (D’Alembert’s paradox). In these animations, red represents regions of low pressure. The left animation shows what the surrounding fluid looks like when the Kutta condition is applied. Notice that the fluid separates smoothly at the trailing edge of the airfoil and a low pressure region is produced on the upper surface of the wing, resulting in lift. The lift is proportional to the circulation around the airfoil. The right animation shows what the surrounding fluid looks like when there is no circulation around the airfoil (stall). Notice the sharp singularity at the trailing edge of the airfoil.

Joukowski Airfoil – C++, 11/8/07

Here is an animation that shows how the streamlines change when you increase the circulation around the airfoil. (Please note: The background fluid motion in this animation is just for effect and is not accurate!) Here is some Mathematica code to plot the streamlines and pressure using Bernoulli’s equation:
```(* runtime: 13 seconds *) U = rho = 1; chord = 4; thk = 0.5; alpha = Pi/9; y0 = 0.2; x0 = -thk/5.2; L = chord/4; a = Sqrt[y0^2 + L^2]; gamma = 4Pi a U Sin[alpha + ArcCos[L/a]]; w[z_, sign_] := Module[{zeta = (z + sign Sqrt[z^2 - 4 L^2])/2}, zeta = (zeta - x0 - I y0)Exp[-I alpha]/Sqrt[(1 - x0)^2 + y0^2]; U(zeta + a^2/zeta) + I gamma Log[zeta]/(2Pi)]; sign[z_] := Sign[Re[z]]If[Abs[Re[z]] < chord/2 && 0 < Im[z] < 2y0(1 - (2Re[z]/chord)^2), -1, 1];w[z_] := w[z, sign[z]]; V[z_] = D[w[z, sign], z] /. sign -> sign[z]; << Graphics`Master`; DisplayTogether[DensityPlot[-0.5rho Abs[V[(x + I y)Exp[I alpha]]]^2, {x, -3, 3}, {y, -3, 3}, PlotPoints -> 275, Mesh -> False, Frame -> False, ColorFunction -> (If[# == 1, Hue[0, 0, 0], Hue[(5# - 1)/6]] &)],ContourPlot[Im[w[(x + I y)Exp[I alpha]]], {x, -3, 3}, {y, -3, 3}, Contours -> Table[x^3 + 0.0208, {x, -2, 2, 0.1}], PlotPoints -> 100, ContourShading -> False], AspectRatio -> Automatic];```

24
Oct
07

### Magnus Effect

A Curveball is a spinning ball that accelerates sideways due to the Magnus Effect. Here is some Mathematica code to plot the streamlines around a lifting cylinder, assuming inviscid incompressible potential flow (irrotational) :
```(* runtime: 0.3 second *) U = 1; gamma = 4; a = 1; w[z_] := U(z + a^2/z) + I gamma Log[z]/(2Pi); ContourPlot[Im[w[x + I y]], {x, -2, 2}, {y, -2, 2}, Contours -> Table[psi, {psi, -3, 3, 0.25}], PlotPoints -> 100, ContourShading -> False, AspectRatio -> Automatic]```

Here is some Mathematica code to plot the entrained fluid using the 4th order Runge-Kutta method:
```(* runtime: 17 seconds, increase n for better resolution *) Clear[v]; n = 40; U = 1; gamma = 4; a = 1; tmax = 5.0; dt = tmax/100; v[z_] := If[Abs[z] < 1, 0, Conjugate[U (1 - (a/z)^2) + I gamma/(2Pi z)]]; image = Table[z = x + I y; Do[k1 = dt v[z]; k2 = dt v[z + k1/2]; k3 = dt v[z + k2/2]; k4 =dt v[z + k3]; z -= (k1 + 2 k2 + 2 k3 + k4)/6, {t, tmax, 0, -dt}]; z, {y, -2, 2, 4.0/n}, {x, -2, 2, 4.0/n}]; ListDensityPlot[Map[Abs[Mod[#, 1]] &, image, {2}], Mesh -> False, Frame -> False, AspectRatio -> 1]```

Link: My Educational Project – some basic potential flows using Matlab and Mathematica

31
Dec
05

### Leap-Frogging Bubble Rings

Here is a simple technique for modeling vortices assuming inviscid incompressible potential flow (irrotational). This technique was inspired from Kerry Mitchell’s paper based on Ultra Fractal. You can also see this here on Dynamical Systems’ Picture Gallery.

Here is some Mathematica code to plot the entrained fluid:
``` (* runtime: 25 seconds, increase n for better resolution *) tmax = 0.85; rcore = 0.1; Klist = {1, -1, 1, -1}; zlist0 = {-1 - 0.5I, -1 + 0.5I, -0.5 - 0.5I, -0.5 + 0.5I}; m = Length[zlist0]; v[K_, z_, z0_] := Module[{r2 = Abs[z - z0]^2}, I K(z0 - z)/r2(1 - Exp[-r2/rcore^2])]; zlist = NDSolve[Flatten[Table[{Subscript[z,i]'[t] == Sum[If[i == j, 0, v[Klist[[j]], Subscript[z, i][t], Subscript[z, j][t]]], {j, 1, m}], Subscript[z, i][0] ==zlist0[[i]]}, {i, 1, m}]], Table[Subscript[z, i][t], {i, 1, m}], {t, 0, tmax}][[1, All, 2]]; n = 23; image = Table[NDSolve[{z'[t] == Sum[v[Klist[[i]], z[t], zlist[[i]]], {i, 1, m}], z[tmax] == x + I y}, z[t], {t, 0, tmax}, MaxSteps -> 5000][[1, 1, 2]] /. t -> 0, {y, -1.125, 1.125, 2.25/n}, {x, -0.35, 1.9, 2.25/n}]; ListDensityPlot[Map[Sign[Im[#]]Arg[#] &, image, {2}], Mesh -> False, Frame -> False, AspectRatio -> Automatic]```

16
Dec
05

### NACA Airfoil

Here is a 9415 NACA airfoil. It’s almost shaped the same as the Joukowski airfoil (but it’s a little different):
```(* runtime: 0.01 second *) n = 40; c = 0.09; p = 0.4; t = 0.15; Clear[x]; camber := c If[x < p, (2p x - x^2)/p^2, ((1 - 2p) + 2p x - x^2)/(1 - p)^2]; theta = ArcTan[D[camber, x]]; p = Table[x = 0.5(1 - Cos[Pi s]); x1 = 1.00893x; thk = 5t(0.2969Sqrt[x1] - 0.126x1 - 0.3516x1^2 + 0.2843x1^3 - 0.1015x1^4); {x, camber} + Sign[s]thk {-Sin[theta],Cos[theta]}, {s, -1, 1, 2/(n - 1)}]; ListPlot[p, PlotJoined -> True, AspectRatio -> Automatic];```

We can approximate the pressure profile using the vortex panel method assuming inviscid incompressible potential flowirrotational). The following Mathematica code was adapted from my Fortran program for my 4/25/99 research project on optimal airfoil design:``` (* runtime: 0.3 second *) alpha = Pi/9; pc = Table[(p[[i]] + p[[i + 1]])/2, {i, 1, n - 1}]; s = Table[v = p[[i + 1]] - p[[i]];Sqrt[v.v], {i, 1, n - 1}]; theta = Table[v = p[[i + 1]] - p[[i]]; ArcTan[v[[1]], v[[2]]], {i, 1, n - 1}]; sin = Sin[theta]; cos = Cos[theta]; Cn1 = Cn2 = Ct1 = Ct2 =Table[0, {n - 1}, {n - 1}]; Do[If[i == j, Cn1[[i, j]] = -1; Cn2[[i, j]] = 1; Ct1[[i, j]] = Ct2[[i, j]] = Pi/2,v = pc[[i]] - p[[j]]; a = -v.{cos[[j]], sin[[j]]}; b = v.v;t = theta[[i]] - theta[[j]];c = Sin[t]; d = Cos[t]; e = v.{sin[[j]], -cos[[j]]}; f =Log[1 + s[[j]](s[[j]] + 2a)/b]; g = ArcTan[b + a s[[j]], e s[[j]]];t = theta[[i]] - 2theta[[j]]; p1 = v.{Sin[t], Cos[t]}; q = v.{Cos[t], -Sin[t]}; Cn2[[i, j]] = d + (0.5q f - (a c + d e)g)/s[[j]]; Cn1[[i, j]] = 0.5d f + c g - Cn2[[i, j]];Ct2[[i, j]] = c + 0.5p1 f/s[[j]] + (a d - c e)g/s[[j]]; Ct1[[i, j]] = 0.5c f - d g - Ct2[[i, j]]], {i, 1, n - 1}, {j, 1, n - 1}]; gamma = LinearSolve[Table[If[i == n, If[j == 1 || j == n, 1, 0], If[j == n, 0, Cn1[[i, j]]] + If[j == 1, 0, Cn2[[i, j - 1]]]], {i, 1, n}, {j, 1, n}], Table[If[i ==n, 0, Sin[theta[[i]] - alpha]], {i, 1, n}]]; ListPlot[Table[q = Cos[theta[[i]] - alpha] + Sum[(If[j == n, 0, Ct1[[i, j]]] + If[j == 1, 0, Ct2[[i, j - 1]]])gamma[[j]], {j, 1, n}]; {pc[[i, 1]], 1 - q^2}, {i, 1, n - 1}], PlotJoined -> True];```

Here is some code for a pretty pressure plot:
```(* runtime: 33 seconds *) w[z_] := z Exp[-I alpha] + I Sum[s[[j]]gamma[[j]]Log[z - pc[[j, 1]] - I pc[[j, 2]]], {j, 1, n - 1}]; V[z_] = D[w[z], z]; DensityPlot[-Abs[V[(x + I y)Exp[I alpha]]]^2, {x, -0.25, 1.25}, {y, -0.75, 0.75}, PlotPoints -> 275, Mesh -> False, Frame -> False, ColorFunction -> (Hue[(5# - 1)/6] &), AspectRatio -> Automatic];```

## 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