TL;DR 我一直在实现一个 python 程序，以使用 runge-kutta 4 和拍摄方法基于特定的相似性变量来求解自然对流的数值方程。但是，当我绘制它时，我没有得到正确的解决方案。我在某个地方犯了错误吗？
你好 ！
从自然对流的一个特例开始，我们得到了这些相似方程。
第一个描述流体流动，第二个描述热流。
“Pr”代表普朗特，它基本上是流体动力学（普朗特）中使用的无量纲数：
这些方程服从以下边界值，使得板附近的温度大于边界层外的温度，并且远离边界层的流体速度为 0。
我一直在尝试使用 Runge-Kutta 4 和拍摄方法以数值方式解决这些问题，以将边界值问题转换为初始值问题。射击方法的实现方式是用牛顿法。
但是，我没有得到正确的解决方案。正如您在下面看到的，随着我们远离盘子，温度（红色）正在增加，而它应该呈指数下降。流体速度（蓝色）更加一致，但是我认为它应该上升得更快然后下降得更快。这里的曲线更平滑。
现在，事实是我们有一个 2 耦合 ODE 系统。但是，现在，我只是试图找到两个初始值之一（例如 f''(0) = a，试图找到 a），这样我们就可以解决边界值问题（拍摄方法） . 一旦找到，我想我们就有了解决整个问题的办法。
我想我也许应该管理这两个 (f''(0) = a ; theta'(0) = b) 但我不知道如何同时管理这两个。最后想提一下，如果我试图获得 theta' 的初始值（所以 theta'(0)），我没有得到正确的热量分布。
这是代码：
"""
The goal is to resolve a 3rd order non-linear ODE for the blasius problem.
It's made of 2 equations (flow / heat)
f''' = 3ff'' - 2(f')^2 + theta
3 Pr f theta' + theta'' = 0
RK4 + Shooting Method
"""
import numpy as np
import math
from scipy.integrate import odeint
from scipy.optimize import newton
from edo_solver.plot import plot
from constants import PRECISION
def blasius_edo(y, t, prandtl):
f = y[0:3]
theta = y[3:5]
return np.array([
# flow edo
f[1], # f' = df/dn
f[2], # f'' = d^2f/dn^2
- 3 * f[0] * f[2] + (2 * math.pow(f[1], 2)) - theta[0], # f''' = - 3ff'' + 2(f')^2 - theta,
# heat edo
theta[1], # theta' = dtheta/dn
- 3 * prandtl * f[0] * theta[1], # theta'' = - 3 Pr f theta'
])
def rk4(eta_range, shoot):
prandtl = 0.01
# initial values
f_init = [0, 0, shoot] # f(0), f'(0), f''(0)
theta_init = [1, shoot] # theta(0), theta'(0)
ci = f_init + theta_init # concatenate two ci
# note: tuple with single argument must have "," at the end of the tuple
return odeint(func=blasius_edo, y0=ci, t=eta_range, args=(prandtl,))
"""
if we have :
f'(t_0) = fprime_t0 ; f'(eta -> infty) = fprime_inf
we can transform it into :
f'(t_0) = fprime_t0 ; f''(t_0) = a
we define the function F(a) = f'(infty ; a) - fprime_inf
if F(a) has a root in "a",
then the solutions to the initial value problem with f''(t_0) = a
is also the solution the boundary problem with f'(eta -> infty) = fprime_inf
our goal is to find the root, we have the root...we have the solution.
it can be done with bissection method or newton method.
"""
def shooting(eta_range):
# boundary value
fprimeinf = 0 # f'(eta -> infty) = 0
# initial guess
# as far as I understand
# it has to be the good guess
# otherwise the result can be completely wrong
initial_guess = 10 # guess for f''(0)
# define our function to optimize
# our goal is to take big eta because eta should approach infty
# [-1, 1] : last row, second column => f'(eta_final) ~ f'(eta -> infty)
fun = lambda initial_guess: rk4(eta_range, initial_guess)[-1, 1] - fprimeinf
# newton method resolve the ODE system until eta_final
# then adjust the shoot and resolve again until we have a correct shoot
shoot = newton(func=fun, x0=initial_guess)
# resolve our system of ODE with the good "a"
y = rk4(eta_range, shoot)
return y
def compute_blasius_edo(title, eta_final):
ETA_0 = 0
ETA_INTERVAL = 0.1
ETA_FINAL = eta_final
# default values
title = title
x_label = "$\eta$"
y_label = "profil de vitesse $(f'(\eta))$ / profil de température $(\\theta)$"
legends = ["$f'(\eta)$", "$\\theta$"]
eta_range = np.arange(ETA_0, ETA_FINAL + ETA_INTERVAL, ETA_INTERVAL)
# shoot
y_set = shooting(eta_range)
plot(eta_range, y_set, title, legends, x_label, y_label)
compute_blasius_edo(
title="Convection naturelle - Solution de similitude",
eta_final=10
)