I would like to solve a 1D Navier equation on a cylindrical imposed tubes(cartesian cordinates).
流动沿 y 方向，右室压力 p1 和左室压力 p2 初始均为 8e5 Pa。圆柱壁静止，宽度“l”的活塞移动，位移 x=asin(omegat) 给定速度 u=aomega余弦（欧米茄t) 我怀疑这会给我带来错误。
沿间隙的速度分布为 u1，在 x 空间中离散。沿 y 的压力变化假设为线性 dp/dy=(p1-p2+pf)/l，其中 pf 是由于间隙上的摩擦而造成的压力损失。
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)
但我没有得到想要的结果。速度 u1 收敛得很好。 Qgap 有方便的值。压力p1和p2不符合期望。
从物理意义上来说，p1 应该从 8e5 增加到某个值，然后减少到 8e5 以下，导致 p2 做相反的事情。 但是，当运行代码并将 p1 和 p2 绘制在一起时，显示它们都具有正弦变化，但 p2>p1。
请帮我。
如果结果不是您所期望的，但方程是正确的，则可能是积分精度存在数值问题。我建议尝试以下潜在解决方案的顺序：
npt = 501
time = np.linspace(0.25/freq,tf,npt)`
m.options.IMODE=7用于基于时间的顺序模拟。如果问题规模太大，这可以提高求解速度。它产生相同的解决方案m.options.IMODE=4但可以更快。放disp=False以避免因打印求解器输出以进行顺序模拟而导致速度减慢。delta=0.007
npx=151
dx=np.float32(delta/(npx-1))
With npx=201它达到了方程长度限制，其中符号形式为>15,000人物。除此之外，您还需要将梯形规则写为gekko方程，使用m.sum()以减小方程的大小。
m.options.NODES=4或更多（最多 6 个）以获得更高的准确性。NODES=3应该足够了，所以这是最后要尝试的事情。求解器报告成功的解，因此所有方程都满足所需的公差。
--------- 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
---------------------------------------------------
检查模型文件gk0_model.apm（文件位于文件夹中print(m.path) or m.open_folder()打开运行文件夹）确认等式：
m.Equation(Qgap==scipy.integrate.trapz(2*np.pi*xpos*np.array(u1)))`
兼容于gekko。它产生了与梯形规则相结合的象征形式。
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)))
这是最终的形式，可以解决17.8 sec并返回更准确的成功解决方案。
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)