# Calculating a Shock using Newton method and Van der Waals

In this post we are going to calculate a vertical (compression) shock. Therefore, the equation system is solved using the Newton-Raphson method amongst other methods. Furthermore, we use the Van der Waals equation of state to calculate fluid properties.

### Introduction to Shocks

A shock is a discontinuous change of the flow condition. This phenomena can just occur in a compressible fluid. During a shock, the static pressure rises significantly, while the velocity drops. Therefore, a vertical shock transforms a supersonic flow into a subsonic flow. Hence, this is the considered case here.

### Basic equations

Notation...- – Specific gas constant in J/kg/K
- – Pressure in Pa
- – Temperature in K
- – Specific volume in m³/kg
- – Density in kg/m³
- – Critical temperature in K
- – Critical pressure in Pa
- – Velocity in m/s
- – Direction in m

First of all, let’s assume the following: **1-dimensional-**, **stationary flow** and a pipe with a **constant diameter**. Therefore, the governing equations can be described as follows:

**Continuity equation**(differential)

(1)

**Continuity equation**(algebraic)

(2)

**Energy conservation**(differential)

(3)

**Energy conservation**(algebraic)

(4)

Furthermore, you can look up the coherence between static and stagnation state here. Let us continue with the governing equations:

**Momentum conservation**(differential)

(5)

Finally, let us neglect the volume forces e.g. gravitational acceleration and the friction forces . The result is:

**Momentum conservation**(Algebraic)

(6)

### Van der Waals Equation of State

Equations of state are used to calculate fluid properties. For example, the specific volume can be calculated from pressure and temperature . The Van der Waals equation is used to take the real gas behaviour into account. Hence, it is an improvement compared to the the ideal gas law. Let us have a look at the basic form of the Van der Waals equation:

(7)

Here, we can calculate the variable from:

(8)

and from:

(9)

Afterwards, we rearrange it with respect to .

(10)

Now we can see that it is a cubic equation. Hence, we need a numerical or analytical method to solve it. Additionally, we can solve the basic equation directly for the pressure p = p(T,v):

(11)

In addition, **caloric state variables** like enthalpy , entropy and inner energy can be calculated using departure functions. The derivation is wonderfully described in this this YouTube video. Hence, the the enthalpy for a Van der Waals gas can be calculated from:

(12)

Here, is the enthalpy of an ideal gas. Hence, we can calculate it using the specific heat capacity and the temperature :

(13)

Therefore, we need to define a reference state to solve the integral. **Remember: since the reference states can differ, a single enthalpy is not worth much, always look at enthalpy differences.** Furthermore, we need temperature dependend specific heat capacity of an ideal gas . For instance, you can use a polynomial to approximate the specific heat capacity.

Do do that, we can use the Python function *integrate* from the scipy package. As an example, let us calculate the enthalpy difference of two thermodynamic states for water (steam). Therefore, we use

- The critical point of water = 647.1 K and = 220.6e5 Pa
- Specific gas constant = 461.5 J/kg/K
- Reference temperature = 300 K
- Polynomial for the specific heat capacity of water (
**mind the temperature range!**) - State 1: 500 K, 2 m³/kg
- State 2: 600 K, 2.5 m³/kg

Hence, in Python this could look like the following.

#License: MIT License (http://opensource.org/licenses/MIT)

import scipy.integrate as integrate

#Water

Tc = 647.1 #K

pc = 220.6e5 #Pa

R = 461.5 #J/kg/K

T_ref = 300 #K

def a():

return (27./64.)*R**2*Tc**2./pc

def b():

return 1./8.*R*Tc/pc/8.

def get_cp_t(T):

#Used to approximate the ideal specific heat capacity of Water (Steam)

#Applicable from 300 - 1000 K

return -0.0000001989*T**3 + 0.0005871194*T**2 - 0.1303510876*T + 1770.6334862186

def idealEnthalpy(T):

return integrate.quad(get_cp_t, T_ref, T)

def vanDerWaals_h_vt(v, T):

hIdeal = idealEnthalpy(T)[0]

return (-2.*a()/v+(b()*R*T)/(v-b())) + hIdeal

T1 = 500 #K

T2 = 600 #K

v1 = 2 #m^3/kg

v2 = 2.5 #m^3/kg

h1 = vanDerWaals_h_vt(v1, T1)

h2 = vanDerWaals_h_vt(v2, T2)

dh = h1-h2

print dh

As a result, we get the enthalpy difference .

### Solving a cubic Equation using the Newton-Raphson Method

First of all, we can use the Newton-Raphson method to solve nonlinear equations or equation systems numerically. This is achieved by finding the roots of a function. Hence, with the Newton-Raphson method the roots of a function can be successively approximated. We can describe it as follows:

(14)

### Solving a cubic Equation in Python

Let us solve the cubic equation for the specific volume using the Newton-Raphson method in Python.

Therefore, we need to know, that a cubic equation has three roots. In the two-phase area we usually have 3 real roots. Thereby, the smallest root represents the volume of the liquid, while the largest root represents the volume of the vapor. The root in between is physically not relevant.

Furthermore, we can have one real root for gas or liquid phase and two imaginary roots. Hence the real root is of relevance.

Let us apply this using Python:

#License: MIT License (http://opensource.org/licenses/MIT)

from scipy.misc import derivative

from scipy import optimize

#Water

Tc = 647.1 #K

pc = 220.6e5 #Pa

R = 461.5 #J/kg/K

def a():

return (27./64.)*R**2*Tc**2./pc

def b():

return 1./8.*R*Tc/pc/8.

def vanDerWaals_v_pt(v, p, T):

return p*v**3.-(R*T+b()*p)*v**2.+a()*v-a()*b()

def getSpecVolume(vGuess, p, T, Nmax=100, eps=1e-6):

it = 0

fv = vanDerWaals_v_pt(vGuess, p, T)

dfv = derivative(vanDerWaals_v_pt, vGuess, dx = 1e-6, args = (p, T,))

vNew = vGuess - fv/dfv

while it < Nmax:

v = vNew

fv = vanDerWaals_v_pt(v, p, T)

dfv = derivative(vanDerWaals_v_pt, v, dx = 1e-10, args = (p, T,))

vNew = v - fv/dfv

if abs(v-vNew)/v < eps:

break

it += 1

return vNew

p = 1e5 # Pa

T = 800 # K

vguess = 3 #m³/kg

v = getSpecVolume(vguess, p, T)

vTest = optimize.newton(vanDerWaals_v_pt, vguess, args=(p, T, ))

As a result, we obtain a specific volume of **3.68758 m³/kg**. Furthermore, the result from the own function *getSpecVolume* is confirmed by the *newton function* from the *scipy optimize*

package.

### Vertical (compression) shock

Furthermore, let us continue with the equations to describe shock phenomena. Therefore, let us derive the revelant equations to model the vertical compression shock. By inserting *Eq. 2* into *Eq. 4* we obtain:

(15)

Furthermore, by inserting *Eq. 2* into *Eq. 6* we obtain:

(16)

Therefore, we can solve this implicit coupled equation system for and .

### Solve the equation system using Python

The equation system *eqSystem* which contains the shock equations can be solved using the *fsolve *function from the *scipy optimize*package. In addition, this packages uses methods such as the Newton-Raphson.

#License: MIT License (http://opensource.org/licenses/MIT)

import numpy as np

from scipy import optimize

def area(D):

return math.pi*pow(D,2)*0.25

def eqSystem(x, args):

TDown = x[0]

vDown = x[1]

pUp = args[0]

hUp = args[1]

uUp = args[2]

vUp = args[3]

f = np.zeros(2)

f[0] = vanDerWaals_h_vt(vDown, TDown) - hUp+0.5*uUp**2.*((vDown/vUp)**2.-1.)

f[1] = vanDerWaals_p_vt(vDown, TDown) - pUp+uUp**2./vUp*(vDown/vUp-1.)

return f

p1 = 1e5 # Pa

T1 = 800 # K

v1guess = 3 #m³/kg

v1 = optimize.newton(vanDerWaals_v_pt, v1guess, args=(p1, T1, ))

mDot = 2 #kg/s

d = 0.1 #m

G = mDot/area(d)

h1 = vanDerWaals_h_vt(v1, T1)

u1 = G*v1

x0 = [1500, 10] #Guess values of T2, v2

args = [p1, h1, u1, v1] #Arguments (constant)

sol = optimize.fsolve(eqSystem, x0, args=args, xtol=1e-6, maxfev=200)

v2 = sol[1]

T2 = sol[0]

p2 = vanDerWaals_p_vt(v2, T2)

u2 = u1/v1*v2

### Results

In order to produce some output data, we need to specify some **input data**:

- Fluid: water
- Upstream pressure = 1e5 Pa
- Upstream temperature = 800 K
- Diameter = 0.1 m
- Variation of the velocity 700 – 1200 m/s.

Using the input data, we can run the script and plot the following figures. In the first figure we see the upstream Mach number over the downstream Mach number and the pressure ratio.

In the second figure we see the upstream Mach number over the temperature and specific volume ratio.

We can see, that we have a supersonic upstream state (Ma > 1). Furthermore, by looking at the figures, we can state the following:

- Subsonic state after the shock — < 1
- Increase in pressure after the shock —
- Increase in temperature after the shock —
- Decreasing specific volume after the shock —
- Increase in entropy —

.

### Conclusion

A shock can occur in a compressible fluid. During a shock, the static pressure rises significantly, while the velocity drops. This means, a supersonic flow is transformed into a subsonic flow. The equation system to model a shock can be obtained from a simple energy and momentum balance and the continuity equation. In order to solve this non-linear equation system we use the Newton-Raphson method which included in Python packages.

Furthermore, the Van der Waals equation can be used to receive fluid properties. Therefore, we need to solve a cubic equation using the Newton-Raphson method. Furthermore, we need departure functions to transform thermal state variables into caloric state variables. Therefore, we need to integrate fluid properties from polynomials.

Finally, in a vertical shock the pressure and Temperature increases non linear, while the Mach number and specific volume decreases.

### What did you think?

I’d like to hear what you think about this post.

Let me know by leaving a comment below.

Thank you so much for the great article, it was fluent and to the point. Cheers.

I am working on something similar using the SKR equations for the cubic EoS. I noticed that your first graph for Pr vs M1, the Pr doesnt go to 1 when M1. I had the same issue initially. I eventually figured out the built in solver for the system of equations was at fault. There were two points where the equations crossed, one corresponding to what the value should be and another I am still trying to figure out the significance of mathematically.