解决1D Navier Stokes具有可压缩质量保护（液压阻尼器）的问题
我想在圆柱施加的管（笛卡尔山脉）上求解一个1D Navier方程。
沿y方向的流动为p1，右腔室最初都在8e5 pa处。圆柱壁是固定的，宽度'L'的活塞正在用A移动位移x = asin（omega t）给定u = a omega cos（omega t）的速度，我怀疑会给我一个错误。
沿间隙的速度轮廓为U1，在X空间中离散。沿Y沿Y的压力变化假定PF的线性DP/DY =（P1-P2+PF）/L是由于间隙上的摩擦而导致的压力损失。
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＆gt; p1。
请帮我。
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.
如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱，如果其他用户或者作者回复了您的评论，将不能在第一时间通知您！
发布评论
评论（1）
如果结果不是您的期望，而是方程式是正确的，那么集成精度可能是数值问题。我建议尝试此潜在解决方案的顺序：
m.options.imode = 7用于基于时间的顺序模拟。如果问题大小太大，则可以提高解决方案速度。它产生与m.options.imode = 4的解决方案相同，但可以更快。 SETdisp = false避免由于打印求解器输出以进行顺序模拟而避免放缓。使用
npx = 201它达到方程式长度限制，其中符号表格为＆gt; 15,000字符。除此之外，您还需要使用m.sum（）将梯形规则编写为gekko等式，以减少方程式。m.options.nodes = 4或更多（最多6），以提高准确性。nodes = 3应该足够，因此这是最后一件事。求解器报告了成功的解决方案，因此所有方程式都可以满足所需的公差。
检查模型文件
gk0_model.apm（文件在文件夹中）确认方程：与
gekko兼容。它产生与梯形规则的符号形式。这是在
17.8 sec中求解的最终形式，并以更准确的速度返回了成功的解决方案。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:
m.options.IMODE=7for time-based sequential simulation. This can improve the solution speed if the problem size is too large. It produces the same solution asm.options.IMODE=4but can be faster. Setdisp=Falseto avoid a slowdown due to printing the solver output for sequential simulation.With
npx=201it reaches an equation length limit where the symbolic form is>15,000characters. Beyond this, you would need to write the trapezoidal rule as agekkoequation, usingm.sum()to reduce the equation size.m.options.NODES=4or more (up to 6) for more accuracy.NODES=3should 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.
Checking the model file
gk0_model.apm(file is in folderprint(m.path)orm.open_folder()to open the run folder) confirms that the equation:is compatible with
gekko. It produces a symbolic form of integration with the trapezoidal rule.Here is a final form that solves in
17.8 secand returns a successful solution with more accuracy.