We want to reformulate the traditional compressible Navier Stokes equations and the continuity, energy and ideal gas laws in (6) equations and unknowns as (4) equations and unknowns by eliminating variables, etc. The goal is not necessarily to form a smaller equation set, but to cross eliminate factors with the hope of reducing the number of approximations and simplifications that take place in differencing techniques.
Since this is incredibly tedious and error prone, we hope intrepid souls will check the algebra and calculus of this step by step treatment.
We start with the coordinate variables t=time, x, y, z, velocity = (u, v, w) p=pressure rho=density T=temperature nu=viscosity gx, gy, gz=gravity and the usual standard variable nomenclature. Note: the (d) is used here to indicate partial derivatives for lack of a delta symbol, sorry.
|
The standard equations are: Ideal gas p=rho * R * T Continuity d rho/dt + d(rho*u)/dx + d(rho*v)/dy + d(rho*w)/dz = 0 X Momentum rho*(du/dt + u*du/dx + v*du/dy + w*du/dz) = rho*gx - dp/dx + d(2*mu(du/dx) + lambda*divV)/dx + d[mu(du/dy + dv/dx)]/dy + d([mu(dw/dx + du/dz)/dz]
Y Momentum rho*(dv/dt + u*dv/dx + v*dv/dy + w*dv/dz) = rho*gy - dp/dy + d(mu*(dvdx+dudy))/dx + d(2*mu(dvdy + lambda * div V))/dy + d(mu*(dvdz+dw/dy)) Z Momentum rho*(dw/dt + u*dw/dx + v*dw/dy + w*dw/dz) = rho*gz - dp/dz + mu*(d2wdx2 + d2w/dy2 + d2w/dz2)/rho Energy rho*Dh/Dt = Dp/Dt + div(k del T) + phi) where phi= nu*[2*(du/dx)^2 + 2 *(dv/dy)^2 + 2 *(dw/dz)^2 + (dv/dx + du/dy)^2 + (dw/dy + dv/dz)^2 + (du/dz + dw/dx)^2]/rho + lambda*(du/dx + dv/dy + dw/dz)^2 |
|
******************************************************************** The obvious first move is to use the ideal gas law to delete rho from all equations. We could replace p or pressure, but it is experimentally difficult to measure rho. Thus, replacing rho with p/R*T we get: Continuity d (p/(R*T)/dt + d(p/(R*T)*u)/dx + d(p/(R*T)*v)/dy + d(p/(R*T)*w)/dz = 0
rho*gx - dp/dx + mu*[2*d2u/dx2 + (lambda/mu)*(du2/dx2 + dv2/dydx + dw2dzdx) + d2u/dy2 + d2v/dxdy + d2u/dz2 + d2w/dxdz] = rho*(du/dt + u*du/dx + v*du/dy + w*du/dz)
X Momentum p*gx/(R*T) - dp/dx + nu*(d2u/dx2 + d2u/dy2 + d2u/dz2)/(p/(R*T)) = p/(R*T)*(du/dt + u*du/dx + v*du/dy + w*du/dz) Y Momentum p*gy/(R*T) - dp/dy + nu*(d2v/dx2 + d2v/dy2 + d2v/dz2)/(p/(R*T)) = p/(R*T)*(dv/dt + u*dv/dx + v*dv/dy + w*dv/dz) Z Momentum p*gz/(R*T) - dp/dz + nu*(d2wdx2 + d2w/dy2 + d2w/dz2)/(p/(R*T)) = (p/(R*T))*(dw/dt + u*dw/dx + v*dw/dy + w*dw/dz) Energy p*(Dh/Dt)/(R*T) = Dp/Dt + div(k del T) + phi) where phi= nu*[2*(du/dx)^2 + 2 *(dv/dy)^2 + 2 *(dw/dz)^2 + (dv/dx + du/dy)^2 + (dw/dy + dv/dz)^2 + (du/dz + dw/dx)^2]/(p/(R*T)) + lambda*(du/dx + dv/dy + dw/dz)^2 |
|
******************************************************************** Flipping a bunch of division terms around, pulling out constants from derivatives and simplifying some:
Continuity d (p/T)/dt + d((p/T)*u)/dx + d((p/(T)*v)/dy + d((p/(T)*w)/dz = 0
rho*gx - dp/dx + mu*[(2+lambda/mu)*d2u/dx2 + d2u/dy2 + d2u/dz2 + (1+lambda/mu)*(dv2/dydx + dw2dzdx) ] = rho*(du/dt + u*du/dx + v*du/dy + w*du/dz)
X Momentum p*gx/(R*T) - dp/dx + R*T*nu*(d2u/dx2 + d2u/dy2 + d2u/dz2)/p = p*(du/dt + u*du/dx + v*du/dy + w*du/dz)/(R*T)
Y Momentum p*gy/(R*T) - dp/dy + R*T*nu*(d2v/dx2 + d2v/dy2 + d2v/dz2)/p = p*(dv/dt + u*dv/dx + v*dv/dy + w*dv/dz)/(R*T) Z Momentum p*gz/(R*T) - dp/dz + R*T*nu*(d2wdx2 + d2w/dy2 + d2w/dz2)/p = p*(dw/dt + u*dw/dx + v*dw/dy + w*dw/dz)/(R*T) Energy p*(dh/dt + u*dh/dx + v*dh/dy + w*dh/dz)/(R*T) = dp/dt + u*dp/dx + v*dp/dy + w*dp/dz + div(k del T) + phi) expanded substantive derivitives where phi= R*T*nu{2*(du/dx)^2 + 2 *(dv/dy)^2 + 2 *(dw/dz)^2 + (dv/dx + du/dy)^2 + (dw/dy + dv/dz)^2 + (du/dz + dw/dx)^2]/p + lambda*(du/dx + dv/dy + dw/dz)^2
|
|
******************************************************************** Expand continuity derivitives, multiply momentum by R*T, substitute Cp*T for h in energy Continuity (dp/dt)/T + p*d(1/T)/dt + d((p*u)/T)/dx + d((p*v)/T)/dy + d((p*w)/T)/dz = 0
X Momentum p*gx - (R*T)*dp/dx + R^2*T^2*nu*(d2u/dx2 + d2u/dy2 + d2u/dz2)/p = p*(du/dt + u*du/dx + v*du/dy + w*du/dz) Y Momentum p*gy - (R*T)*dp/dy + R^2*T^2*nu*(d2v/dx2 + d2v/dy2 + d2v/dz2)/p = p*(dv/dt + u*dv/dx + v*dv/dy + w*dv/dz) Z Momentum p*gz - (R*T)*dp/dz + R^2*T^2*nu*(d2wdx2 + d2w/dy2 + d2w/dz2)/p = p*(dw/dt + u*dw/dx + v*dw/dy + w*dw/dz) Energy Cp*p*(dT/dt + u*dT/dx + v*dT/dy + w*dT/dz)/(R*T) = dp/dt + u*dp/dx + v*dp/dy + w*dp/dz + div(k del T) + phi) where phi= R*T*nu{2*(du/dx)^2 + 2 *(dv/dy)^2 + 2 *(dw/dz)^2 + (dv/dx + du/dy)^2 + (dw/dy + dv/dz)^2 + (du/dz + dw/dx)^2]/p + lambda*(du/dx + dv/dy + dw/dz)^2 |
|
******************************************************************** Expand continuity derivitives again, multiply momentum by p, expand div(k del T) in energy and substitute phi.
Continuity (dp/dt)/T + p*d(1/T)/dt + p*d(u/T)/dx + u*d(p/T)/dx + (1/T)*d(p*u)/dx + p*d(v/T)/dy + v*d(p/T)/dy + (1/T)*d(p*v)/dy + p*d(w/T)/dz + w*d(p/T)/dz + (1/T)*d(p*w)/dz = 0
X Momentum p^2*gx - (p*R*T)*dp/dx + R^2*T^2*nu*(d2u/dx2 + d2u/dy2 + d2u/dz2) = p^2*(du/dt + u*du/dx + v*du/dy + w*du/dz) Y Momentum p^2*gy - (p*R*T)*dp/dy + R^2*T^2*nu*(d2v/dx2 + d2v/dy2 + d2v/dz2) = p^2*(dv/dt + u*dv/dx + v*dv/dy + w*dv/dz) Z Momentum p^2*gz - (p*R*T)*dp/dz + R^2*T^2*nu*(d2wdx2 + d2w/dy2 + d2w/dz2) = p^2*(dw/dt + u*dw/dx + v*dw/dy + w*dw/dz)
Energy Cp*p*(dT/dt + u*dT/dx + v*dT/dy + w*dT/dz)/(R*T) = dp/dt + u*dp/dx + v*dp/dy + w*dp/dz + k*(d2 T/dx2 + d2T/dy2 + d2T/dz2) + R*T*nu{2*(du/dx)^2 + 2 *(dv/dy)^2 + 2 *(dw/dz)^2 + (dv/dx + du/dy)^2 + (dw/dy + dv/dz)^2 + (du/dz + dw/dx)^2]/p + lambda*(du/dx + dv/dy + dw/dz)^2 |
|
******************************************************************** Expand continuity derivitives again, rearrange momentum Continuity (1/T)*(dp/dt) + p*d(1/T)/dt + p*u*d(1/T)/dx + p*(1/T)*du/dx + u*p*d(1/T)/dx + u*(1/T)*dp/dx+ (1/T)*p*du/dx + (1/T)*u*dp/dx + p*v*d(1/T)/dy + p*(1/T)*dv/dy + v*p*d(1/T)/dy + v*(1/T)*dp/dy+ (1/T)*p*dv/dy + (1/T)*v*dp/dy + p*w*d(1/T)/dz + p*(1/T)*dw/dz + w*p*d(1/T)/dz + w*(1/T)*dp/dz+ (1/T)*p*dw/dz + (1/T)*w*dp/dz = 0
X Momentum -p^2*(du/dt + u*du/dx + v*du/dy + w*du/dz - gx) - p*(R*T)*dp/dx + R^2*T^2*nu*(d2u/dx2 + d2u/dy2 + d2u/dz2) = 0 Y Momentum -p^2*(dv/dt + u*dv/dx + v*dv/dy + w*dv/dz - gy) - p*(R*T)*dp/dy + R^2*T^2*nu*(d2v/dx2 + d2v/dy2 + d2v/dz2) = 0 Z Momentum -p^2*(dw/dt + u*dw/dx + v*dw/dy + w*dw/dz - gz) - p*(R*T)*dp/dz + R^2*T^2*nu*(d2wdx2 + d2w/dy2 + d2w/dz2) = 0
Expand continuity derivitives again, move pRT*(dp/dn) terms in momentum to other side. Continuity (1/T)*(dp/dt) - (1/T^2)*p*dT/dt
- (1/T^2)*p*u*dT/dx + p*(1/T)*du/dx - (1/T^2)*u*p*dT/dx + u*(1/T)*dp/dx+ (1/T)*p*du/dx + (1/T)*u*dp/dx
- (1/T^2)p*v*dT/dy + p*(1/T)*dv/dy - (1/T^2)*v*p*dT/dy + v*(1/T)*dp/dy+ (1/T)*p*dv/dy + (1/T)*v*dp/dy
- (1/T^2)*p*w*dT/dz + p*(1/T)*dw/dz - (1/T^2)*w*p*dT/dz + w*(1/T)*dp/dz+ (1/T)*p*dw/dz + (1/T)*w*dp/dz
X Momentum -p^2*(du/dt + u*du/dx + v*du/dy + w*du/dz - gx) + R^2*T^2*nu*(d2u/dx2 + d2u/dy2 + d2u/dz2) = p*(R*T)*dp/dx Y Momentum -p^2*(dv/dt + u*dv/dx + v*dv/dy + w*dv/dz - gy) + R^2*T^2*nu*(d2v/dx2 + d2v/dy2 + d2v/dz2) = p*(R*T)*dp/dy Z Momentum -p^2*(dw/dt + u*dw/dx + v*dw/dy + w*dw/dz - gz) + R^2*T^2*nu*(d2wdx2 + d2w/dy2 + d2w/dz2) = p*(R*T)*dp/dz |
|
******************************************************************** Multiply continuity by T, divide momentum by p*R*T Continuity dp/dt - (1/T)*p*dT/dt
- (1/T)*p*u*dT/dx + p*du/dx - (1/T)*u*p*dT/dx + u*dp/dx+ p*du/dx + u*dp/dx
- (1/T)*p*v*dT/dy + p*dv/dy - (1/T)*v*p*dT/dy + v*dp/dy+ p*dv/dy + v*dp/dy
- (1/T)*p*w*dT/dz + p*dw/dz - (1/T)*w*p*dT/dz + w*dp/dz+ p*dw/dz + w*dp/dz = 0
X Momentum -p*(du/dt + u*du/dx + v*du/dy + w*du/dz - gx)/(R*T) + (R*T/p)*nu*(d2u/dx2 + d2u/dy2 + d2u/dz2) = dp/dx Y Momentum -p*(dv/dt + u*dv/dx + v*dv/dy + w*dv/dz - gy)/(R*T) + (R*T/p)*nu*(d2v/dx2 + d2v/dy2 + d2v/dz2) = dp/dy Z Momentum -p*(dw/dt + u*dw/dx + v*dw/dy + w*dw/dz - gz)/(R*T) + (R*T/p)*nu*(d2wdx2 + d2w/dy2 + d2w/dz2) = dp/dz |
|
******************************************************************** Rearrange continuity by partial derivatives, Continuity dp/dt - (1/T)*p*dT/dt
- 2*(1/T)*p*u*dT/dx + 2*u*dp/dx+ 2*p*du/dx
- 2*(1/T)*p*v*dT/dy + 2*v*dp/dy+ 2*p*dv/dy
- 2*(1/T)*p*w*dT/dz + 2*w*dp/dz+ 2*p*dw/dz = 0 X Momentum -p*(du/dt + u*du/dx + v*du/dy + w*du/dz - gx)/(R*T) + (R*T/p)*nu*(d2u/dx2 + d2u/dy2 + d2u/dz2) = dp/dx Y Momentum -p*(dv/dt + u*dv/dx + v*dv/dy + w*dv/dz - gy)/(R*T) + (R*T/p)*nu*(d2v/dx2 + d2v/dy2 + d2v/dz2) = dp/dy Z Momentum -p*(dw/dt + u*dw/dx + v*dw/dy + w*dw/dz - gz)/(R*T) + (R*T/p)*nu*(d2wdx2 + d2w/dy2 + d2w/dz2) = dp/dz |
|
******************************************************************** Rearrange continuity by partial derivatives, Continuity - 2 * (1/T) * p * ( (1/2)*dT/dt + u*dT/dx + v*dT/dy + w*dT/dz) + 2 * ((1/2) * dp/dt + u*dp/dx + v*dp/dy + w*dp/dz) + 2 * p * (du/dx + dv/dy + dw/dz) = 0 |
|
******************************************************************** Divide continuity by two and by p, Continuity - (1/T) * ( (1/2) * dT/dt + u * dT/dx + v * dT/dy + w * dT/dz) + (1/p) * ( (1/2) * dp/dt + u * dp/dx + v * dp/dy + w * dp/dz) + (du/dx + dv/dy + dw/dz) = 0 |
|
******************************************************************** Start to solve continuity for dp/dt, Continuity - (1/T) * ( (1/2) * dT/dt + u * dT/dx + v * dT/dy + w * dT/dz) + (1/p) * ( u * dp/dx + v * dp/dy + w * dp/dz) + (du/dx + dv/dy + dw/dz) = -1/(2*p)*dp/dt |
|
******************************************************************** Multiply continuity by -2*p, Continuity -2*p* [ - (1/T) * ( (1/2) * dT/dt + u * dT/dx + v * dT/dy + w * dT/dz) + (1/p) * ( u * dp/dx + v * dp/dy + w * dp/dz) + (du/dx + dv/dy + dw/dz)] = dp/dt |
|
******************************************************************** Separate continuity by partials of p, Continuity -2*p* [ - (1/T) * ( (1/2) * dT/dt + u * dT/dx + v * dT/dy + w * dT/dz) + (du/dx + dv/dy + dw/dz)] -2*( u * dp/dx + v * dp/dy + w * dp/dz) = dp/dt Energy |
|
******************************************************************* Inserting dp/dt from continuity into the energy equation Energy Cp*p*(dT/dt + u*dT/dx + v*dT/dy + w*dT/dz)/(R*T) = |
|
******************************************************************* Combining the DT/Dt terms and simplifying: Energy 0 = ((p/T)- (Cp*p)/(R*T))*dT/dt
|
|
******************************************************************* Multiplying by p and pulling out p where indicated Energy 0 =p^2 * [ ((1/T)- (Cp)/(R*T))*dT/dt
|
|
******************************************************************* Now we insert dp/dx, dp/dy, dp/dz from the momentum equations Energy 0 =p^2 * [ ((1/T)- (Cp)/(R*T))*dT/dt |
|
******************************************************************* Collecting terms of p and simplifying Energy 0 =p^2 * [ ((1/T)- (Cp)/(R*T))*dT/dt |
|
******************************************************************* Pull out the 1/R*T term Energy 0 =p^2 *(1/(R*T))* [ (R- Cp)*dT/dt |
|
******************************************************************* Collect Terms in the p^2 coefficient Energy 0 =p^2*(1/(R*T))* [ |
|
********************************************************************** The crowning achievment is to solve the quadratic equation for p=[-b+-(b^2-4ac)^.5 ]/ 2a
|
|
Simplifying the quadratic equation for p=[-b+-(b^2-4ac)^.5 ]/ 2a
|