A Mathematica-Based

Hydrodynamics Model

George E. Hrabovsky

CSC'10 Las Vegas, 13 July 2010

What I Will Cover In This Talk

Why Hydrodynamics?

Why Hydrodynamics code?

Why Mathematica?

The Continuity Equation

The Momentum Equation

The Thermal Energy Equation

The Simple Test Model

Conclusions

The Future

References

Why Hydrodynamics?

Solids             Liquids         Gases           Plasmas
                                                                                                                                These Are Considered To Be Fluids

Many exotic forms of matter, such as superfluids, Bose-Einstein condensates, quark-gluon plasmas, and so on, can also be considered to be fluids.

Why Hydrodynamics Code?

({{Model-Type, Governing Mathematics, Applicable Solution Methods}, {Classical Mechanics of a Single Particle, F = m (d^2)x/(dt^2), Second-order differential equations where F (x, v, t)., Analytical methods appropriate for linear ODEs, 
approximation methods, and numerical methods.}, {Systems of Particles, Motion of the center of mass., Analytical methods appropriate for linear ODEs, 
approximation methods, and numerical methods.}, {Ideal Fluids "Von Neumann's Dry Water", The Continuity Equation
The Momentum Equation
The Thermal Energy Equation, These PDEs have only limited analytic solvability,
Approximation  by special functions/integral transforms/and series expansions is more common,
and numerical methods.}, {Viscous Fluids, The Continuity Equation
The Navier-Stokes Equation
The Thermal Energy Equation, These PDEs have only a few limited analytic solutions,
Approximation more often by special functions/integral transforms/and series expansions and even these are limited,
and numerical methods}})

Why Mathematica?

({{FORTRAN, MATLAB, Mathematica}, {, , }, {1 page of thirty, 1 page of six, The entire model}})

Speed of calculation is nearly that of compiled code, now Mathematica has a built-in compiler.

Built-in visualization commands and built-in/default parallel support for multiprocessing/multithreading.

The Continuity Equation

(∂_t)ρ+(∂_(Subscript[x, i]))ρ (Subscript[u, i])=0.

Time rate of      Advection of density

change of density

(∂_t)ρ+u (∂_x)ρ   = 0.

(∂_t)ρ[x,t]+(∂_x)(ρ[x,t] u[x,t])==0

The Navier-Stokes and Momentum Equations

ρ((∂_t)(Subscript[u, i])+(Subscript[u, j])(∂_(Subscript[x, j]))(Subscript[u, i]))=-(∂_(Subscript[x, i]))p+(Subscript[B, i])+(∂_(Subscript[x, j]))[μ((∂_(Subscript[x, j]))(Subscript[u, i])+(∂_(Subscript[x, i]))(Subscript[u, j])-(2/3)(Subscript[δ, ij])(∂_(Subscript[x, k]))(Subscript[u, k]))]
+(∂_(Subscript[x, i]))(ζ(∂_(Subscript[x, k]))(Subscript[u, k])).

mass × acceleration                       force

ρ((∂_t)u+u(∂_x)u)=-(∂_x)p+(Subscript[B, x])+(∂_x)(μ(4/3)(∂_x)u)+(∂_x)(ζ(∂_x)u).

(∂_t)u[x,t]+u[x,t] (∂_x)u[x,t]==-(∂_x)p[x,t]/ρ[x,t]

The Energy Equation

ρ (Subscript[c, v])((∂_t)T+(Subscript[u, i])(∂_(Subscript[x, i]))T)=-p(∂_(Subscript[x, i]))(Subscript[u, i])+κ(∂_((Subscript[x, i])(Subscript[x, i])))T-(∂_(Subscript[x, i]))(Subscript[q, i])+(Subscript[Φ, i])+q'''.

radiation           internal heat

flux vector            generation per volume

(Subscript[q, i])=0, and q'''=0.

(Subscript[Φ, i])=(Subscript[σ, ij])'(∂_(Subscript[x, i]))(Subscript[u, j]).

ρ (Subscript[c, v])((∂_t)T+(Subscript[u, i])(∂_(Subscript[x, i]))T)=-p(∂_(Subscript[x, i]))(Subscript[u, i])+κ(∂_((Subscript[x, i])(Subscript[x, i])))T+(Subscript[σ, ij])'(∂_(Subscript[x, i]))(Subscript[u, j]).

The Energy Equation 2

ρ (Subscript[c, v])((∂_t)T+u(∂_x)T)=-p(∂_x)u+κ(∂_(x,x))T+2 ((μ((∂_x)u))^2)+((ζ((∂_x)u))^2).

density/specific heat/total            pressure          thermal                      friction                 viscosity

emperature derivative                  advection        advection                   advection             advection

(∂_t)p+u (∂_x)p=((5/3p)/ρ)((∂_t)ρ+u (∂_x)ρ).

(∂_t)p[x,t]+u[x,t] (∂_x)p[x,t]==
((5/3 p[x,t])/(ρ[x,t]))((∂_t)ρ[x,t]+u[x,t] (∂_x)ρ[x,t])

A Simple Model-Boundary and Initial Conditions

Boundary Conditions:

{ρ[0,t]==ρ[30,t],
p[0,t]==p[30,t],u[0,t]==u[30,t]}

Initial Conditions:

{ρ[x,0]==2,
p[x,0]==10+Sin[2 Pi x/30],u[x,0]==0}

A Simple Model-Building the Model

soln=NDSolve[{D[ρ[x,t],t]+D[ρ[x,t] u[x,t],x]==0, 
D[u[x,t],t]+u[x,t] D[u[x,t],x]==-D[p[x,t],x]/ρ[x,t],
D[p[x,t],t]+u[x,t] D[p[x,t],x]==
((5/3 p[x,t])/(ρ[x,t]))(D[ρ[x,t],t]+u[x,t] D[ρ[x,t],x]),ρ[0,t]==ρ[30,t],
p[0,t]==p[30,t],u[0,t]==u[30,t],ρ[x,0]==2,
p[x,0]==10+Sin[2 Pi x/30],u[x,0]==0},{ρ[x,t],p[x,t],u[x,t]},
{x,0,30},{t,0,50},Method->{"ExplicitRungeKutta","DifferenceOrder"->8}]

A Simple Model-The Results

We can plot the density, pressure, and velocity.

GraphicsGrid[Plot3D[ρ[x,t]/.soln,{x,0,30},{t,0,50},PlotLabel->Density],Plot3D[u[x,t]/.soln,{x,0,30},{t,0,50},PlotLabel->Velocity],Plot3D[p[x,t]/.soln,{x,0,30},{t,0,50},PlotLabel->Pressure]]

,,

Conclusions

Good for prototyping and algorithm-development on small machines.

Good for small/medium-scale modeling on powerful desktops/work-stations.

Good for any level of modeling on multiprocessors systems/clusters.

The Future

More complicated models in 2D and 3D.

Thunderstorm/tornado modeling.

Astrophysical applications.

References

Thank You Very Much!

George E. Hrabovsky