TL;DR I've been implementing a python program to solve numerically equations for natural convection based on a particular similarity variable using runge-kutta 4 and the shooting method. However I don't get the right solutions when I plot it. Did I make a mistake somewhere ?
嗨
从自然对流的特殊情况出发,我们得到了这些相似方程。
第一个描述流体流动,第二个描述热流。
“Pr”是指Prandtl,它基本上是流体动力学中使用的无量纲数(Prandtl):
这些方程受以下边界值的约束,因此板附近的温度大于边界层外的温度,并且流体速度远离边界层为0
我一直在尝试用Runge Kutta 4和打靶法数值求解这些问题,将边值问题转化为初值问题。打靶法的实现方式是牛顿法
然而,我没有得到正确的解决方案。 正如你在下面所看到的,当我们离开极板时,温度(红色)正在升高,而它应该呈指数下降。 它对于流体速度(蓝色)更一致,但是我认为它应该上升得更快,然后下降得更快。这里的曲线更平滑
事实上,我们有一个由2个耦合的ODE组成的系统。然而,现在,我只想找到两个初始值中的一个(例如,f''(0)=a,试图找到a),这样我们就有了一个边值问题的解决方案(shooting method)。一旦找到,我想我们就有了解决整个问题的办法
我想我应该管理这两个(f'(0)=a;θ(0)=b),但我不知道如何并行管理这两个。 最后想一提的是,如果我试图得到θ′(所以θ′(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
)
您需要使用android.print包中的类
您需要使打印文档适配器类派生自PrintDocumentAdapter 这个类负责你打印的文档内容
要执行打印,请执行以下操作:
注:安卓。打印要求API级别19及更高
相关问题 更多 >
编程相关推荐