python - 用射击法求解自然对流方程（热量和流动）

问题描述

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
)

标签： pythonnumerical-methodsodefluid-dynamics

python - 用射击法求解自然对流方程（热量和流动）

问题描述

解决方案

推荐阅读