热工数值方法-改进尤拉方法

首先，什么是尤拉方法？尤拉方法的思想是利用近似值y_n代替y(x_n)的值，然后利用曲线倒数乘以一个设定的步长得出y_n+1的值，当我们设定的步长比较小的时候，所得的拟合曲线便会有足够的精度，然而为进一步提高计算的精度，结合梯形方法即可以形成改进尤拉格式，这样在进行热工数值计算时能够使精度足够且计算量较小。

在这种改进的尤拉格式中，先用尤拉方法求得一个初步的近似值，记y_n+1，这一步即预报，然后再通过梯形方法对这一预报值进行修正，得到较为精确的y_n+1值，这是一种热工数值计算中常用的方法。

这一方法可以表示为以下的两式：

• 预报: y_n+1 = y_n + hf (x_n , y_n )
• 校正: y_n+1 = y_n + h*[f (x_n , y_n )+f (x_n+1 , y_n+1 )]/2

在计算机求解的时候，我们需要给出一个y关于x的函数，同时设定步长h，这个函数不管是显式或者隐式的函数都可以。

Matlab中利用改进尤拉方法计算传热问题举例：

题：薄壁物体温度响应的数值求解。如图所示，T_∞=1500℃，T_w=500 ℃。从某时刻起，T₀=20℃，直径1mm的热电偶结点置于高温气体中，气流对热电偶的表面换热系数 h=200W(m²/K)，结点表面发射率 ε=0.2，高温气流对金属的氧化作用使热电偶结点产生一定的热量，其强度为 q=1000W/m3，求热电偶结点的温度响应。结点材料：ρ=8000kg/m³，c=0.38kJ/(kg*K)，λ=50W/(m*K)

解：我们在这里使用改进尤拉方法进行热工数值计算。

Matlab源代码如下：

function [x,y]=eulerg(x,y,h,N)
x=zeros(1,N+1);
y=zeros(1,N+1);
x(1)=0;
y(1)=293;
for n=1:N
x(n+1)=x(n)+h;
ybar=y(n)+h*dyfun(x(n),y(n));
y(n+1)=y(n)+(h/2)*(dyfun(x(n),y(n))+dyfun(x(n),ybar));
end
plot(x,y,'red');
function m=dyfun(x,y)
m=1/304+75/19*(1773-y)+(6.804e-8)/304*(1773^4-y^4);

在 matlab 中新建文件写入代码保存后，输入例如：eulerg(0,293,0.001,1000) 即可以得出热电偶温度y随时间x的变化规律曲线。