A Baby Hydrodynamics Model

George E. Hrabovsky

Christopher J. Winfield

James Firmiss

Introduction

This is the first attempt at a model of hydrodynamics in Mathematica for the MAST hydrocode. This is a very simple problem and is more of a gesture than a full-blown model. I will make no attempt to discuss stability issues, or code artifacts, or accuracy and error at this stage. Consider it a means of getting used to the apparatus in a lab before performing a full-blown experiment.

Some Fluid Dynamics

These are the standard conventions and equations for fluid dynamics [1]. Here is a table of variables.

({{Symbol, Meaning}, {u(x,t), Velocity field depending upon the position x and the time t.}, {p(x,t), Pressure field depending upon the position x and the time t.}, {ρ(x,t), Density field depending upon the position x and the time t.}, {∇, The nabla operator, for Cartesian coordinates this is ((∂/(∂x)),(∂/(∂y)),(∂/(∂z))).}, {(D/Dt)[], Total time derivative of the field [], (∂/(∂t))[]+u·∇[].}, {(Subscript[Γ, 1]), The ratio of the specific heats at constant pressure and at constant volume.}})

An idealized fluid element, essentially a point that moves with the bulk velocity of the fluid u is called a material point. This concept deserves more description [2]. A material point, or fluid element, is a region over which we can define some variable, and is a sort of mathematical convenience. The following criteria apply:

The first equation of motion we will write is the continuity equation.

((∂ρ)/(∂t))+∇·(ρu)=0.

See appendix 1 for details. The second equation of motion is the equation of motion,

((∂u)/(∂t))+u·∇u=-((∇p)/ρ).

details of this equation can be found in appendix 2. The third equation is the thermal energy equation,

((∂p)/(∂t))+u·∇p=(((Subscript[Γ, 1])p)/ρ)(((∂ρ)/(∂t))+u·∇ρ).

Details of this can be found in Appendix 3.

Coding the Equations

In one spatial dimension, the equations of motion become, in Mathematica,

Γ1=((cp/cv))
coneq=D[ρ[x,t],t]+D[ρ[x,t] u[x,t],x]==0
eom=D[u[x,t],t]+u[x,t] D[u[x,t],x]==-D[p[x,t],x]/ρ[x,t]
theq=D[p[x,t],t]+u[x,t] D[p[x,t],x]==
((Γ1 p[x,t])/(ρ[x,t]))(D[ρ[x,t],t]+u[x,t] D[ρ[x,t],x])

The gamma expression needs to be calculated for each case.

The Boundary Conditions for the Problem

We will have periodic boundary conditions, with sinusoidal initial conditions of pressure and density, initial velocity will be zero everywhere.

Coding the Boundary Conditions

The boundary conditions are:

ρ[0,t]==ρ[10,t]
p[0,t]==p[10,t]
u[0,t]==u[10,t]

the initial conditions are:

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

The Model

We will assume (Subscript[Γ, 1])=5/3. Putting the equations together with the boundary and initial conditions we have the following:

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

The Results

Plot3D[ρ[x,t]/.soln,{x,0,30},{t,0,50},PlotPoints->100]

Plot3D[p[x,t]/.soln,{x,0,30},{t,0,50},PlotPoints->100]

Plot3D[u[x,t]/.soln,{x,0,30},{t,0,50},PlotPoints->100]

The Model 2

We will assume (Subscript[Γ, 1])=5/3. Putting the equations together with the boundary and initial conditions we have the following:

delta[x_,a_]:=1/(a Sqrt[π]) Exp[-x^2/a^2]

soln2=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+delta[(x-15),0.5],p[x,0]==2+delta[(x-15),0.5],u[x,0]==2+delta[(x-15),0.5]},{ρ[x,t],p[x,t],u[x,t]},{x,0,30},{t,0,5},Method->{"ExplicitRungeKutta","DifferenceOrder"->8},PrecisionGoal->2,MaxStepFraction->1]

Options[NDSolve]

soln21=NDSolve[(*continuity equation*){D[ρ[x,t],t]+D[ρ[x,t] u[x,t],x]==0,(*equation of motion*)D[u[x,t],t]+u[x,t] D[u[x,t],x]==-D[p[x,t],x]/ρ[x,t],(*thermal equation*)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]),(*boundary conditions*)ρ[0,t]==0,ρ[30,t]==0,p[0,t]==0,p[30,t]==0,u[0,t]==0,u[30,t]==0,(*initial conditions*)ρ[x,0]==5+delta[(x-15),0.5],p[x,0]==5+delta[(x-15),0.5],u[x,0]==5+delta[(x-15),0.5]},(*plot parameters*){ρ[x,t],p[x,t],u[x,t]},{x,0,30},{t,0,50},Method->{"ExplicitRungeKutta","DifferenceOrder"->8},PrecisionGoal->2,MaxStepFraction->1]

The Results 2

Plot3D[ρ[x,t]/.soln2,{x,0,30},{t,0,5}]

Plot3D[ρ[x,t]/.soln21,{x,0,30},{t,0,15}]

Plot3D[p[x,t]/.soln2,{x,0,30},{t,0,5}]

Plot3D[p[x,t]/.soln21,{x,0,30},{t,0,15}]

Plot3D[u[x,t]/.soln2,{x,0,30},{t,0,5}]

Plot3D[u[x,t]/.soln21,{x,0,30},{t,0,15}]

This is taken, with modifications, from [3]. Given an initial volume (Subscript[V, 0]) the mass of fluid within that volume will be,

∫ρdV

with ρ the density. The mass will flow in unit time through a surface element surrounding the volume, dS. This mass will be,

ρ u · dS,

where the magnitude of the vector dS is equal to the surface area. By convention we take dS to be along the outward normal. Thus, ρ u · dS is positive for mass flowing out of the volume, and negative flowing into it. The total mass flowing out of the volume is then

∮ρ u · dS.

Now decrease, per unit time, the mass in the volume. (4) then becomes,

-(∂/(∂t))∫ρdV.

Since this must also equal the mass flowing out of the volume, we now equate this to (5),

(∂/(∂t))∫ρdV=-∮ρ u · dS.

This surface integral is transformed via Green's formula to a volume integral,

∮ρ u · dS=∫∇·(ρu)dV.

We can rewrite (6),

(∂/(∂t))∫ρdV+∫∇·(ρu)dV=0.

or,

∫[((∂ρ)/(∂t))+∇·(ρu)]dV=0.

Since this holds for any volume then the integrand must vanish, so we have the condition,

((∂ρ)/(∂t))+∇·(ρu)=0.

This is the equation of continuity.

We begin by stating Newton's second law of motion,

m ((du)/(dt))=F

-∮pdS

the integral of the pressure taken over the surface bounding the volume of the fluid element. Transforming this into a volume integral we have,

-∮pdS=-∫∇p dV.

ρ ((Du)/Dt)=-∇p.

((Du)/Dt)=((∂u)/(∂t))+u·∇u.

ρ ((∂u)/(∂t))+ρu·∇u=-∇p.

((∂u)/(∂t))+u·∇u=-((∇p)/ρ).

See reference [4].