Archive for the 'fluids' Category

10
Dec
08

Smoothed Particle Hydrodynamics (SPH)

Smoothed Particle Hydrodynamics (SPH)

Links

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

Links

22
Oct
07

Leap-Frogging Vortex Rings

Mathematica 4.2 version: 12/31/05;

Here is some Mathematica code to numerically solve this using the 4th order Runge-Kutta method:

(* runtime: 21 seconds, increase n for better resolution *)
n = 23; tmax = 0.85; dt = tmax/100; rcore = 0.1; Klist = {1, -1, 1, -1}; zlist = {-1 - 0.5I, -1 + 0.5I, -0.5 - 0.5I, -0.5 + 0.5I}; m = Length[zlist];
v[K_, z_, z0_] := Module[{r2 = Abs[z - z0]^2}, I K(z0 - z)/r2(1 - Exp[-r2/rcore^2])];
Runge[z_] := Module[{}, k1 = dt vtot[z]; k2 = dt vtot[z + k1/2]; k3 = dt vtot[z + k2/2]; k4 = dt vtot[z + k3]; (k1 + 2 k2 + 2 k3 + k4)/6];
vtot[zlist_] := Table[Sum[If[i == j, 0, v[Klist[[j]], zlist[[i]], zlist[[j]]]], {j, 1, m}], {i, 1, m}]; zlists = Transpose[Table[zlist += Runge[zlist], {t, 0, tmax, dt}]];
vtot[z_] := Plus @@ v[Klist, z, zlists[[All, t]]]; image = Table[z = x + I y; Do[z -= Runge[z], {t, 100, 1, -1}]; z, {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]

Links

15
Mar
07

1D Shock Tube


A shock tube is a tube containing high and low pressure gas separated by a thin diaphragm. A shock wave is produced when the diaphragm is quickly removed. The color in the upper plot shows the pressure. The lower plot shows the density. The following Mathematica code solves Euler’s equations using the finite volume method with the Jameson-Schmidt-Turkel (JST) scheme and Runge-Kutta time stepping.


(* runtime: 5 seconds *)
gamma = 1.4;
R[W_] := Module[{}, rho = W[[All, 1]]; u = W[[All, 2]]/rho; p = (gamma - 1)(W[[All, 3]] - rho u^2/2); F = u W + Transpose[{Table[0, {n}], p, u p}]; h = Table[(F[[Min[n, i + 1]]] + F[[i]])/2, {i,1, n}]; Q = Table[h[[Max[i, 2]]] - h[[Max[i, 2] - 1]], {i, 1, n}]; nu = Table[i = Max[2, Min[n - 1, i]]; Abs[(p[[i + 1]] - 2p[[i]] + p[[i - 1]])/(p[[i + 1]] + 2p[[i]] + p[[i - 1]])], {i, 1, n}]; S = Table[Max[nu[[Min[n, i + 1]]], nu[[i]]], {i, 1, n}]; alpha1 = 1/2; beta1 = 1/4;alpha2 = alpha1; beta2 = beta1; epsilon2 = Map[Min[alpha1, alpha2#] &, S];epsilon4 = Map[Max[0, beta1 - beta2#] &, epsilon2];dW = Table[W[[Min[n - 1, i] + 1]] - W[[Min[n - 1, i]]], {i, 1, n}];dW3 = Table[i = Max[2, Min[n - 2, i]]; -W[[i - 1]] + 3W[[i]] - 3W[[i + 1]] + W[[i + 2]], {i, 1, n}];d = (epsilon2 dW - epsilon4 dW3)(Abs[u] + a); Dflux = Table[d[[Max[2, i]]] - d[[Max[2, i] - 1]], {i, 1, n}]; (Q - Dflux)/dx];
n = 50; dx = 1.0/n; a = 1.0; dt = dx/(1.0 + a); W = Table[{If[i > n/2, 0.125, 1], 0, If[i > n/2, 0.1, 1]/(gamma - 1)}, {i, 1, n}];
Do[W -= dt R[W - dt R[W - dt R[W]/4]/3]/2;ListPlot[W[[All, 1]], PlotJoined -> True,PlotRange -> {0, 1}, AxesLabel -> {"i", "rho"}], {t, 0, 100dt, dt}];

01
Feb
07

Supersonic Flow

Here are some weird tests of CFL3D. The color represents the Mach number.

Supersonic Flow – CFL3D (Computational Fluids Laboratory 3D), Tecplot 360

Links

16
Nov
06

Kelvin-Helmholtz Instability Waves

Kelvin-Helmholtz Instability Waves are formed when one fluid layer is on top of another fluid layer moving with a different velocity. These instabilities take the form of small waves that can eventually grow into vortex rollers. This is a purely potential flowirrotational) phenomenon. The following Mathematica code was adapted from Zheming Zheng’s Fortran program. It uses vortex blobs to simulate smooth rollers and periodic point vortices to simulate the far field. A simple predictor-corrector scheme is used to integrate the differential equations:

(* runtime: 9 minutes *) n = 40; dt = 0.05; L = 2Pi; gamma = L/n; rcore = 0.5; plist = Table[{x, 0.01 Sin[2Pi x/L]}, {x, 0, L(1 - 1/n), L/n}]; vcalc[plist_] := Table[Sum[{x, y} = plist[[j]] - plist[[i]]; gamma/(2Pi) Sum[If[i == j && k == 0, 0, r2 = (x + k L)^2 +y^2; {-y, x + k L}/r2(1 - Exp[-r2/rcore^2] - 1)], {k, -3, 3}] + If[i ==j, 0, gamma/(2L) {-Sinh[2Pi y/L], Sin[2Pi x/L]}/(Cosh[2Pi y/L] - Cos[2Pi x/L])], {j, 1, n}], {i, 1, n}]; Do[vlist = vcalc[plist]; vlist0 =vlist; vlist = vcalc[plist + dt vlist]; plist += dt(vlist0 + vlist)/2; ListPlot[plist, PlotJoined -> True, PlotRange -> {-2, 2}, AspectRatio -> Automatic], {i, 1, 400}];

Here is some Mathematica code to plot streamlines assuming periodic point vortices: (* runtime: 2 seconds *) Clear[psi]; psi[x_, y_] := Plus @@ Map[-gamma/(4 Pi) Re[Log[Cos[2Pi(x - #[[1]])/L] – Cosh[2Pi (y - #[[2]])/L]]] &, plist]; ContourPlot[psi[x, y], {x, 0, L}, {y, -L/3, L/3}, PlotRange -> All, ContourShading -> False, PlotPoints -> 20, AspectRatio -> Automatic]

Link: Kelvin-Helmholtz waves in the clouds

Here is another variation showing vortex pairing and pathlines.

Link: movie – famous experiment by Bernal

03
Nov
06

Folding Wing Aircraft

This is a preliminary test of a proposed project to simulate the wing fluttering of an aircraft with folding wings. This was solved using the finite volume method for 3D, steady, compressible, viscous, laminar flow on a curvilinear grid. This image shows the pressure distribution over the wings.

Folding Wing Aircraft calculated in CFL3D (Computational Fluids Laboratory 3D), rendered in POV-Ray 3.6.1, 11/3/06

Links

21
Oct
06

Vortices

Here is another code to solve Euler’s equation (inviscid flow) and plot the pathlines. This method was adapted from Stephen Montgomery-Smith’sEuler2D program. The vortex effects are found using the Biot-Savart law and the differential equations are solved using the Adams Bashforth method. Click here to download some Matlab code for this. Shown below is some Mathematica code:
(* runtime: 19 seconds *)
ToPoint[z_] := {Re[z], Im[z]};dt = 0.01; rcore = 0.1; dz0 = 0; Klist = {1, -1, 1, -1};
zlist = Join[{-1 - 0.5I, -1 + 0.5I, -0.5 - 0.5I, -0.5 + 0.5I}, Flatten[Table[x + I y, {x, -1, 1, 0.1}, {y, -1, 1, 0.1}], 1]];
paths = Map[# Table[1, {10}] &, zlist];
Do[dz = Map[Sum[zj = zlist[[j]]; If[# != zj, r2 = Abs[# - zj]^2; I Klist[[j]](zj - #)/r2 (1 -Exp[-r2/rcore^2]), 0], {j, 1, Length[Klist]}] &, zlist]; zlist += dt (1.5dz - 0.5dz0); dz0 = dz; paths = Table[Prepend[Delete[paths[[i]], -1], zlist[[i]]], {i, 1, Length[zlist]}];Show[Graphics[Map[Line[Map[ToPoint, #]] &, paths], PlotRange -> {{-1, 1}, {-1, 1}}, AspectRatio -> Automatic]], {85}];

Link: math paper by Stephen Montgomery-Smith, seems fairly advanced to me

28
Apr
06

Flapping Wing

This flapping wing was calculated using the unsteady vortex panel methodAlan Lai’s Fortran code. It assumes inviscid incompressible potential flow (irrotational). I also have a working Mathematica version of this code, but it is a little lengthy to show here.

Flapping Wing – Fortran 90, rendered in POV-Ray 3.6.1, 4/28/06

Links

22
Apr
06

Water Caustics

This water-like surface was generated in Mathematica using frequency filtered random noise, and then it was raytraced in POV-Ray and water caustics were added using Henrik Jensen’sphoton mapping technique. The smooth_triangle command was used for phong normal interpolation of the water’s surface.

Link: Water Caustics Generator – C++ program by Kjell Andersson, uses Heron’s Formula to calculate the caustics




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

  • 490,812 hits

Follow

Get every new post delivered to your Inbox.