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

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