TL; DR Я реализовал программу python для численного решения уравнений естественной конвекции на основе конкретной переменной подобия, используя рунге-кутта 4 и метод стрельбы. Однако я не получаю правильных решений, когда строю его. Я где-то ошибся?
Привет!
Исходя из частного случая естественной конвекции, мы получаем эти уравнения подобия.
Первое описывает поток жидкости, второе - тепловой поток.
«Pr» для Прандтля в основном безразмерное число, используемое в гидродинамике ( Прандтль ):
Эти уравнения подчиняются следующим граничным значениям, так что температура вблизи пластины выше, чем температура за пределами пограничного слоя и такая, что скорость жидкости равна 0 вдали от пограничного слоя.
I ' Мы пытались решить их численно с помощью Рунге-Кутта 4 и метода стрельбы, чтобы преобразовать краевую задачу в задачу с начальным значением. Метод стрельбы реализован с помощью метода Ньютона.
Однако я не получаю правильных решений. Как вы можете видеть ниже, температура (выделенная красным) увеличивается по мере удаления от пластины, тогда как она должна уменьшаться по экспоненте. Он более согласован для скорости жидкости (выделен синим цветом), однако я думаю, что скорость должна go увеличиваться быстрее, чем go быстрее снижаться. Здесь кривая более плавная.
Теперь дело в том, что у нас есть система из двух связанных ОДУ. Однако прямо сейчас я пытаюсь найти одно из двух начальных значений (например, 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
)