Problem on Solving a 1D Navier Stokes with Compressible Mass Conservation (Hydraulic Damper)

Question

I would like to solve a 1D Navier equation on a cylindrical imposed tubes(cartesian cordinates).

The flow is along y direction, with right chamber having pressure p1 and left chamber pressure p2 both initially at 8e5 Pa. The Cylindrical wall is stationary, and the piston of width 'l' is moving with a displacement x=asin(omegat) given a velocity of u=aomegacos(omegat) which I doubt to give me an error.

The velocity profile along the gap is u1, discretised in x space. Pressure variation along y is assumed linear dp/dy=(p1-p2+pf)/l with pf is pressure loss due to friction on the gap.

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate

amp=10e-3
freq=10
tf = 2/freq
omega=2*np.pi*freq
npt = 101
time = np.linspace(0.25/freq,tf,npt)

xglob=amp*np.sin(omega*time)
uglob=amp*omega*np.cos(omega*time)

delta=0.007
npx=101
dx=np.float32(delta/(npx-1))

r1=0.01
r2=0.004
xpos = np.float32(np.linspace(r1,r1+delta,npx))
Ap=np.pi*(r1**2-r2**2)

rho=850
beta=5e8
L=0.06
l=0.01

pressf=12*amp*0.707*uglob*np.cos(2*np.pi*freq*time)/(delta**2)\
        *630*(1+(0.05*(0.707*uglob/delta)**2))**(-0.23)

m = GEKKO()

m.time = time
t=m.Param(m.time)

x=m.MV(ub=amp+1)
x.value=np.ones(npt)*xglob
u=m.MV(ub=10)
u.value=np.ones(npt)*uglob

pf=m.MV(lb=20)
pf.value=np.ones(npt)*pressf

u1=[m.Var(0) for i in range(npx)]
p1=m.Var(800000)
p2=m.Var(800000)

# mu.value[:]=(630 * (1+(0.05*(dudx[:]))**2)**(-0.23))
# mu=m.Intermediate(630 * (1+(0.05*(dudx))**2)**(-0.23))

Qgap=m.Var(value=0)

m.Equation(u1[0].dt()==-1*(p1-p2+pf)/l  \
  + ((630 * (1+(0.05*((u-u1[0])/dx))**2)**(-0.23))/rho)\
  *((u-2*u1[0]+u1[1])/(dx*dx)))
m.Equations([(u1[i].dt()==-1*(p1-p2+pf)/l  \
  + ((630 * (1+(0.05*((u1[i+1]-u1[i-1])/(2*dx)))**2)**(-0.23))/rho)\
  *((u1[i+1]-2*u1[i]+u1[i-1])/(dx*dx))) for i in range(1,npx-1)])
m.Equation(u1[-1].dt()==-1*(p1-p2+pf)/l  \
  + (((630 * (1+(0.05*((u1[-1]-0)/dx))**2)**(-0.23)))/rho)\
  *((u1[-2]-2*u1[-1]+0)/(dx*dx)))
m.Equation(Qgap==scipy.integrate.trapz(2*np.pi*xpos*np.array(u1)))

m.Equation(p1.dt()==(beta/(Ap*(L-x)))*(-Qgap+Ap*u))
m.Equation(p2.dt()==(beta/(Ap*(L+x)))*(Qgap-Ap*u))

# simulation
m.options.IMODE = 4
m.options.solver=1
m.options.nodes=3

m.solve(disp=True)

But I don't get the desired results. The velocity u1 converges quite well. Qgap has convenient values. The pressures p1 and p2 are not as desired.

In a physical sense, p1 should increase from 8e5 to a certain value, then decrease below 8e5 resulting to p2 to do the opposite. But when running the code and plotting p1 and p2 together, it shows they both have got sinusoidal variation but p2>p1.

Please help me.

Answer 1

If the results are not what you expect but the equations are correct, then it could be numerical issues with the integration accuracy. I recommend trying this order of potential solutions:

1. Increase the number of time points.

npt = 501
time = np.linspace(0.25/freq,tf,npt)`

2. Use m.options.IMODE=7 for time-based sequential simulation. This can improve the solution speed if the problem size is too large. It produces the same solution as m.options.IMODE=4 but can be faster. Set disp=False to avoid a slowdown due to printing the solver output for sequential simulation.
3. Increase the number of position points.

delta=0.007
npx=151
dx=np.float32(delta/(npx-1))

With npx=201 it reaches an equation length limit where the symbolic form is >15,000 characters. Beyond this, you would need to write the trapezoidal rule as a gekko equation, using m.sum() to reduce the equation size.

4. Use m.options.NODES=4 or more (up to 6) for more accuracy. NODES=3 should be sufficient so this is a last thing to try.

The solver reports a successful solution so all of the equations are satisfied to the required tolerance.

 --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :            0
   Variables    :          108
   Intermediates:            0
   Connections  :            0
   Equations    :          104
   Residuals    :          104
 
 Number of state variables:          41700
 Number of total equations: -        41700
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :              0
 
 ----------------------------------------------
 Dynamic Simulation with APOPT Solver
 ----------------------------------------------
 
 Iter    Objective  Convergence
    0  3.36523E-08  1.26261E+04
    1  2.02851E-08  4.39248E+03

   40  4.41072E-24  3.38415E-06
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :    29.1031000000003      sec
 Objective      :   0.000000000000000E+000
 Successful solution
 ---------------------------------------------------

Checking the model file gk0_model.apm (file is in folder print(m.path) or m.open_folder() to open the run folder) confirms that the equation:

m.Equation(Qgap==scipy.integrate.trapz(2*np.pi*xpos*np.array(u1)))`

is compatible with gekko. It produces a symbolic form of integration with the trapezoidal rule.

    v104=((((((((((((((((((((((((((((((((((((((((((((((((((((((((((
(((((((((((((((((((((((((((((((((((((((((((((1.0)*((((0.06327167898416519)*
(v2))+((0.06283185631036758)*(v1))))))/(2.0))+((((1.0)*
((((0.0637115016579628)*(v3))+((0.06327167898416519)*(v2))))))/(2.0)))+
((((1.0)*((((0.0641513243317604)*(v4))+((0.0637115016579628)*
(v3))))))/(2.0)))+((((1.0)*((((0.06459114700555801)*(v5))+
...
+((((1.0)*((((0.10681416094303131)*(v101))+((0.1063743382692337)*
(v100))))))/(2.0)))

Here is a final form that solves in 17.8 sec and returns a successful solution with more accuracy.

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate

amp=10e-3
freq=10
tf = 2/freq
omega=2*np.pi*freq
npt = 501
time = np.linspace(0.25/freq,tf,npt)

xglob=amp*np.sin(omega*time)
uglob=amp*omega*np.cos(omega*time)

delta=0.007
npx=151
dx=np.float32(delta/(npx-1))

r1=0.01
r2=0.004
xpos = np.float32(np.linspace(r1,r1+delta,npx))
Ap=np.pi*(r1**2-r2**2)

rho=850
beta=5e8
L=0.06
l=0.01

pressf=12*amp*0.707*uglob*np.cos(2*np.pi*freq*time)/(delta**2)\
        *630*(1+(0.05*(0.707*uglob/delta)**2))**(-0.23)

m = GEKKO(remote=False)

m.time = time
t=m.Param(m.time)

x=m.MV(ub=amp+1)
x.value=np.ones(npt)*xglob
u=m.MV(ub=10)
u.value=np.ones(npt)*uglob

pf=m.MV(lb=20)
pf.value=np.ones(npt)*pressf

u1=[m.Var(0) for i in range(npx)]
p1=m.Var(800000)
p2=m.Var(800000)

# mu.value[:]=(630 * (1+(0.05*(dudx[:]))**2)**(-0.23))
# mu=m.Intermediate(630 * (1+(0.05*(dudx))**2)**(-0.23))

Qgap=m.Var(value=0)

m.Equation(u1[0].dt()==-1*(p1-p2+pf)/l  \
  + ((630 * (1+(0.05*((u-u1[0])/dx))**2)**(-0.23))/rho)\
  *((u-2*u1[0]+u1[1])/(dx*dx)))
m.Equations([(u1[i].dt()==-1*(p1-p2+pf)/l  \
  + ((630 * (1+(0.05*((u1[i+1]-u1[i-1])/(2*dx)))**2)**(-0.23))/rho)\
  *((u1[i+1]-2*u1[i]+u1[i-1])/(dx*dx))) for i in range(1,npx-1)])
m.Equation(u1[-1].dt()==-1*(p1-p2+pf)/l  \
  + (((630 * (1+(0.05*((u1[-1]-0)/dx))**2)**(-0.23)))/rho)\
  *((u1[-2]-2*u1[-1]+0)/(dx*dx)))

m.Equation(p1.dt()==(beta/(Ap*(L-x)))*(-Qgap+Ap*u))
m.Equation(p2.dt()==(beta/(Ap*(L+x)))*(Qgap-Ap*u))

m.Equation(Qgap==scipy.integrate.trapz(2*np.pi*xpos*np.array(u1)))

# simulation
m.options.IMODE = 7
m.options.solver=1
m.options.nodes=3

#m.open_folder()
m.solve(disp=False)

print(m.options.SOLVETIME)
print(m.options.APPSTATUS)