Navier Stokes Reformulation










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

MOMENTUM FULL VERSION

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


MOMENTUM FULL VERSION

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
got rid of 1/R

MOMENTUM FULL VERSION

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)

use above version of momentum perhaps

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)
inverted /p/(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
inverted /p/(R*T)

********************************************************************

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

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

*******************************************************************

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) =

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

+ k*(d2 T/dx2 + d2T/dy2 + d2T/dz2)

+ (R*T*nu/p)*{2*(du/dx)^2 + 2 *(dv/dy)^2 + 2 *(dw/dz)^2

+ (dv/dx + du/dy)^2

+ (dw/dy + dv/dz)^2


+ lambda*(du/dx + dv/dy + dw/dz)^2

*******************************************************************

Combining the DT/Dt terms and simplifying:

Energy

0 = ((p/T)- (Cp*p)/(R*T))*dT/dt

+ ((2*p))/T - (Cp*p)/(R*T)) * (u*dT/dx + v*dT/dy + w*dT/dz)

-2*p* (du/dx + dv/dy + dw/dz)

- (u*dp/dx + v*dp/dy + w*dp/dz)

+ k*(d2 T/dx2 + d2T/dy2 + d2T/dz2)

+ lambda*(du/dx + dv/dy + dw/dz)^2


+ (R*T*nu/p)*{2*(du/dx)^2 + 2 *(dv/dy)^2 + 2 *(dw/dz)^2

+ (dv/dx)^2 + 2*(dv/dx)*(du/dy) + (du/dy)^2

+ (dw/dy)^2 + 2*(dw/dy)*(dv/dz) + (dv/dz)^2

*******************************************************************

Multiplying by p and pulling out p where indicated

Energy

0 =p^2 * [ ((1/T)- (Cp)/(R*T))*dT/dt

+ (2/T - Cp/(R*T) * (u*dT/dx + v*dT/dy + w*dT/dz)

-2* (du/dx + dv/dy + dw/dz)]

+p*[- (u*dp/dx + v*dp/dy + w*dp/dz)

+ k*(d2 T/dx2 + d2T/dy2 + d2T/dz2)

+ lambda*(du/dx + dv/dy + dw/dz)^2]


+ R*T*nu*{2*(du/dx)^2 + 2 *(dv/dy)^2 + 2 *(dw/dz)^2

+ (dv/dx)^2 + 2*(dv/dx)*(du/dy) + (du/dy)^2

+ (dw/dy)^2 + 2*(dw/dy)*(dv/dz) + (dv/dz)^2

*******************************************************************

Now we insert dp/dx, dp/dy, dp/dz from the momentum equations

Energy

0 =p^2 * [ ((1/T)- (Cp)/(R*T))*dT/dt

+ (2/T - Cp/(R*T) * (u*dT/dx + v*dT/dy + w*dT/dz)

-2* (du/dx + dv/dy + dw/dz)]

+p*[- u*[ -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)]

- v* [-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)]

- w*[ -p*(dw/dt + u*dw/dx + v*dw/dy + w*dw/dz - gz)/(R*T) + (R*T/p)*nu*(d2wdx2 + d2w/dy2 + d2w/dz2)]

+ k*(d2 T/dx2 + d2T/dy2 + d2T/dz2)

+ lambda*(du/dx + dv/dy + dw/dz)^2]

+ R*T*nu*{2*(du/dx)^2 + 2 *(dv/dy)^2 + 2 *(dw/dz)^2

+ (dv/dx)^2 + 2*(dv/dx)*(du/dy) + (du/dy)^2

+ (dw/dy)^2 + 2*(dw/dy)*(dv/dz) + (dv/dz)^2

*******************************************************************

Collecting terms of p and simplifying

Energy

0 =p^2 * [ ((1/T)- (Cp)/(R*T))*dT/dt

+ (2/T - Cp/(R*T) * (u*dT/dx + v*dT/dy + w*dT/dz)

- 2*(du/dx + dv/dy + dw/dz)

+ u*(du/dt + u*du/dx + v*du/dy + w*du/dz - gx)/(R*T)

+v*(dv/dt + u*dv/dx + v*dv/dy + w*dv/dz - gy)/(R*T)

+w* (dw/dt + u*dw/dx + v*dw/dy + w*dw/dz - gz)/(R*T) ]

+p*[ k*(d2 T/dx2 + d2T/dy2 + d2T/dz2)

+ lambda*(du/dx + dv/dy + dw/dz)^2]

+ R*T*nu*{

-u*(d2u/dx2 + d2u/dy2 + d2u/dz2)

- v*nu*(d2v/dx2 + d2v/dy2 + d2v/dz2)

- w*(d2wdx2 + d2w/dy2 + d2w/dz2)

+ 2*(du/dx)^2 + 2 *(dv/dy)^2 + 2 *(dw/dz)^2

+ (dv/dx)^2 + 2*(dv/dx)*(du/dy) + (du/dy)^2

+ (dw/dy)^2 + 2*(dw/dy)*(dv/dz) + (dv/dz)^2

*******************************************************************

Pull out the 1/R*T term

Energy

0 =p^2 *(1/(R*T))* [ (R- Cp)*dT/dt

+ (2*R - Cp) * (u*dT/dx + v*dT/dy + w*dT/dz)

- 2*(du/dx + dv/dy + dw/dz)*R*T

+ u*(du/dt + u*du/dx + v*du/dy + w*du/dz - gx)

+v*(dv/dt + u*dv/dx + v*dv/dy + w*dv/dz - gy)

+w* (dw/dt + u*dw/dx + v*dw/dy + w*dw/dz - gz) ]

+p*[ k*(d2 T/dx2 + d2T/dy2 + d2T/dz2)

+ lambda*(du/dx + dv/dy + dw/dz)^2]

+ R*T*nu*{

-u*(d2u/dx2 + d2u/dy2 + d2u/dz2)

- v*(d2v/dx2 + d2v/dy2 + d2v/dz2)

- w*(d2wdx2 + d2w/dy2 + d2w/dz2)

+ 2*(du/dx)^2 + 2 *(dv/dy)^2 + 2 *(dw/dz)^2

+ (dv/dx)^2 + 2*(dv/dx)*(du/dy) + (du/dy)^2

+ (dw/dy)^2 + 2*(dw/dy)*(dv/dz) + (dv/dz)^2

*******************************************************************

Collect Terms in the p^2 coefficient

Energy

0 =p^2*(1/(R*T))* [

u*du/dt + v*dv/dt + w*dw/dt + (R- Cp)*dT/dt

+ (2*R - Cp) * (u*dT/dx + v*dT/dy + w*dT/dz)

- 2*(du/dx + dv/dy + dw/dz)*R*T

+ u^2*du/dx + v^2*dv/dy + w^2*dw/dz

+ u* v*(du/dy+dv/dx) + u*w*(du/dz+dw/dx) + v*w*(dv/dz + dw/dy)

- u*gx - v*gy - w*gz ]

+p*[ k*(d2 T/dx2 + d2T/dy2 + d2T/dz2)

+ lambda*(du/dx + dv/dy + dw/dz)^2]

+ R*T*nu*[

-u*(d2u/dx2 + d2u/dy2 + d2u/dz2)

- v*(d2v/dx2 + d2v/dy2 + d2v/dz2)

- w*(d2wdx2 + d2w/dy2 + d2w/dz2)

+ 2*(du/dx)^2 + 2 *(dv/dy)^2 + 2 *(dw/dz)^2

+ (dv/dx)^2 + 2*(dv/dx)*(du/dy) + (du/dy)^2

+ (dw/dy)^2 + 2*(dw/dy)*(dv/dz) + (dv/dz)^2

**********************************************************************

The crowning achievment is to solve the quadratic equation for p=[-b+-(b^2-4ac)^.5 ]/ 2a


p ={ -[ k*(d2 T/dx2 + d2T/dy2 + d2T/dz2) + lambda*(du/dx + dv/dy + dw/dz)^2]

(+ -)[

[ k*(d2 T/dx2 + d2T/dy2 + d2T/dz2) + lambda*(du/dx + dv/dy + dw/dz)^2]^2

-4*[

u*du/dt + v*dv/dt + w*dw/dt + (R- Cp)*dT/dt

+ (2*R - Cp) * (u*dT/dx + v*dT/dy + w*dT/dz)

- 2*(du/dx + dv/dy + dw/dz)*R*T

+ u^2*du/dx + v^2*dv/dy + w^2*dw/dz

+ u* v*(du/dy+dv/dx) + u*w*(du/dz+dw/dx) + v*w*(dv/dz + dw/dy)

- u*gx - v*gy - w*gz ] / (R*T)]

*[ R*T*nu*[

-u*(d2u/dx2 + d2u/dy2 + d2u/dz2)

- v*(d2v/dx2 + d2v/dy2 + d2v/dz2)

- w*(d2wdx2 + d2w/dy2 + d2w/dz2)

+ 2*(du/dx)^2 + 2 *(dv/dy)^2 + 2 *(dw/dz)^2

+ (dv/dx)^2 + 2*(dv/dx)*(du/dy) + (du/dy)^2

+ (dw/dy)^2 + 2*(dw/dy)*(dv/dz) + (dv/dz)^2


/ 2*[

u*du/dt + v*dv/dt + w*dw/dt + (R- Cp)*dT/dt

+ (2*R - Cp) * (u*dT/dx + v*dT/dy + w*dT/dz)

- 2*(du/dx + dv/dy + dw/dz)*R*T

+ u^2*du/dx + v^2*dv/dy + w^2*dw/dz

+ u* v*(du/dy+dv/dx) + u*w*(du/dz+dw/dx) + v*w*(dv/dz + dw/dy)

- u*gx - v*gy - w*gz ] / (R*T)] **********************************************************************

Simplifying the quadratic equation for p=[-b+-(b^2-4ac)^.5 ]/ 2a


p ={ -[ k*(d2 T/dx2 + d2T/dy2 + d2T/dz2) + lambda*(du/dx + dv/dy + dw/dz)^2]

(+ -)[

[ k*(d2 T/dx2 + d2T/dy2 + d2T/dz2) + lambda*(du/dx + dv/dy + dw/dz)^2]^2

-4*[

u*du/dt + v*dv/dt + w*dw/dt + (R- Cp)*dT/dt

+ (2*R - Cp) * (u*dT/dx + v*dT/dy + w*dT/dz)


+( - 2*R*T + u^2) * du/dx + (-2*R*T + v^2) * dv/dy + (-2*R*T + w^2)* dw/dz


+ u* v*(du/dy+dv/dx) + u*w*(du/dz+dw/dx) + v*w*(dv/dz + dw/dy)

- u*gx - v*gy - w*gz ]]

*[nu*[

-u*(d2u/dx2 + d2u/dy2 + d2u/dz2)

- v*(d2v/dx2 + d2v/dy2 + d2v/dz2)

- w*(d2wdx2 + d2w/dy2 + d2w/dz2)

+ 2*(du/dx)^2 + 2 *(dv/dy)^2 + 2 *(dw/dz)^2

+ (dv/dx)^2 + 2*(dv/dx)*(du/dy) + (du/dy)^2

+ (dw/dy)^2 + 2*(dw/dy)*(dv/dz) + (dv/dz)^2


] ^ .5 }

/ 2*[

u*du/dt + v*dv/dt + w*dw/dt + (R- Cp)*dT/dt

+ (2*R - Cp) * (u*dT/dx + v*dT/dy + w*dT/dz)

- 2*(du/dx + dv/dy + dw/dz)*R*T

+ u^2*du/dx + v^2*dv/dy + w^2*dw/dz

+ u* v*(du/dy+dv/dx) + u*w*(du/dz+dw/dx) + v*w*(dv/dz + dw/dy)

- u*gx - v*gy - w*gz ] / (R*T)]