Prandtl模型流的Python solve_ivp

2024-03-29 10:55:26 发布

您现在位置:Python中文网/ 问答频道 /正文

我正试图解决偏微分方程enter image description here

$\frac{\partial \textbf{u}}{\partial t}=b\sin{\alpha}+ \nu \frac{\partial^2 \textbf{u}}{\partial z^2}$
$\frac{\partial \textbf{b}}{\partial t}=-\textbf{u} N^2\sin{\alpha}+ \frac{\nu}{Pr} \frac{\partial^2 \textbf{b}}{\partial z^2}$

我将PED转换为一个ODE系统,但它无论如何都不会收敛。我使用了每一种解决ivp的方法,但仍然没有收敛,我的代码有什么问题吗?因为它应该有解决方案

我尝试了许多方法,但大多数方法都不收敛,并且会显示溢出

`

def func(t,y):

n = len(y)
# set u,b,k,w as variables
u=y
b=y

# 2 equations initialization
f1=np.zeros(nz)
f2=np.zeros(nz)

#1st equation with BCs
f1[1]= b[1]*m.sin(alpha)+(Nu*(u[2]-2*u[1]+0)/dz**2)
f2[1]= -u[1]*N**2*m.sin(alpha)+(Nu*(b[2]-2*b[1]+bs)/dz**2)
for i in range(2,nz-2):

    f1[i]=b[i]*m.sin(alpha)+((Nu)*(u[i+1]-2*u[i]+u[i-1])/dz**2)

    f2[i]=-u[i]*N**2*m.sin(alpha)+((Nu)*(b[i+1]-2*b[i]+b[i-1])/dz**2)

f1[nz-2]= b[nz-2]*m.sin(alpha)+(Nu*(0-2*u[nz-2]+u[nz-3])/dz**2)
f2[nz-2]= -u[nz-2]*N**2*m.sin(alpha)+(Nu*(0-2*b[nz-2]+b[nz-3])/dz**2)

f1=f1[1:]  
f1=f1[:-1]

f2=f2[1:]  
f2=f2[:-1]
return np.append(f1,f2)

global Beta, alpha, Nu, dz, N, BetaS, Beta, AlphaS, SigmaS, bs #needed in equations

H     = 800               
T     = 1200               
alpha = m.radians(90)      # sloping angle [rad]
bs    = 8.26e-04           # surface buoyancy [m/s^2]
N     = 0.006              # buoyancy frequency [Hz]
Nu    = 1e-05              # kinematic viscosity [m^2/s]
Pr    = 1                  # Prandtl number  


nz = 1200              # number of elements
nt = 6000                # number of time steps


dz = H/nz                  # grid stencil [m]        
dt = T/nt                 # time step [s]




z=np.linspace(0,H,nz)

print(dt)

# I.Cs

S0 = m.sqrt(N*N*m.sin(alpha)**2/Pr)
SC = m.sqrt(S0/2/Kb)#*100.25))
u_0=np.zeros(nz)   
b_0=np.zeros(nz)  

"""Initial condition for U,b   Prantdl model analytical solution with constant viscosity"""
for n in range(nz):            
    u_0[n] = (bs/N/m.sqrt(Pr))*m.exp(-SC*z[n])*m.sin(SC*z[n]);
    b_0[n] = bs*m.exp(-SC*z[n])*m.cos(SC*z[n]);


u_0=u_0[1:]  
u_0=u_0[:-1]

b_0=b_0[1:]  
b_0=b_0[:-1]


s=z[1:]  
s=s[:-1]

plt.plot(u_0, s, label='u')
plt.legend()
plt.show()    

plt.plot(b_0, s, label='b')
plt.legend()
plt.show()  
IC=np.append(u_0,b_0)
#Record start time
start_time=time.time()    

YS=solve_ivp(func, [0,T], IC,method='LSODA', t_eval=np.arange(0,T,dt),  rtol=1e-3)

#Record end time
end_time=time.time() 

默认方法将立即显示溢出,我唯一可以使用的方法是'LSODA',但它将运行很长时间而不会结束