(1) 物理问题
1. 热 平 衡 方 程 1. 热平衡方程 1.热平衡方程
V ρ c v d T d τ = Q 1 − Q 2 V{\rho}c_{v}\frac{dT}{d\tau}=Q_{1}-Q_{2} VρcvdτdT=Q1−Q2
- 产热 Q 1 Q_{1} Q1为: Q 1 = k 0 e − E R T ⋅ C n ⋅ V ⋅ Q Q_{1}=k_{0}e^{-\frac{E}{RT}}{\cdot}C^{n}{\cdot}V{\cdot}Q Q1=k0e−RTE⋅Cn⋅V⋅Q
- 散热 Q 2 Q_{2} Q2为: Q 2 = h ⋅ A ⋅ ( T − T 0 ) Q_{2}=h{\cdot}A{\cdot}(T-T_{0}) Q2=h⋅A⋅(T−T0)
- 联立解得: d T d τ = k 0 e − E R T ⋅ C n ⋅ V ⋅ Q − h ⋅ A ⋅ ( T − T 0 ) V ρ c v \frac{dT}{d{\tau}}=\frac{k_{0}e^{-\frac{E}{RT}}{\cdot}C^{n}{\cdot}V{\cdot}Q-h{\cdot}A{\cdot}(T-T_{0})}{V{\rho}c_{v}} dτdT=Vρcvk0e−RTE⋅Cn⋅V⋅Q−h⋅A⋅(T−T0)
这个方程手工求解很困难,但是用计算机编程的方法求解则很容易。
取初值如下:
- k 0 k_{0} k0——频率因子,初值取100.0
- E E E——活化能,初值取1.E5
- R R R——气体常数,初值取8.314
- C C C——可燃混合物中反应物浓度,初值取1.0
- n n n——反应级数,初值取1.0
- V V V——容器体积,初值取1.0
- Q Q Q——可燃混合物的燃烧热,初值取2.E7
- T T T——容器内可燃混合物温度,初值取与 T 0 T_{0} T0相等的数值。
- h h h——散热系数,初值取5.0
- A——容器壁散热面积,初值取6.0
- T 0 T_{0} T0——容器壁温度,初值取800.0
- ρ \rho ρ——混合气密度,初值取1.0
- c v c_{v} cv——比热,初值取1.E3
初始条件为: τ = 0 , T = T 0 \tau=0,T=T_{0} τ=0,T=T0
(2) 求解目标
- 分别令系统初温 T 0 T_{0} T0=300K、400K、500K、600K、700K、800K、900K、1000K、1100K(其他参数不变),获得不同初始温度下的系统升温曲线。
- 分别令散热系数 h h h=1.0~10.0( T 0 T_{0} T0保持800K),获得不同初始温度下的系统升温曲线。
(3) Euler法求解微分方程
设有一个微分方程:
d y d x = f ( x , y ( x ) ) \frac{dy}{dx}=f(x,y(x)) dxdy=f(x,y(x))定解条件为:
x = x 0 , y = y 0 x=x_{0},y=y_{0} x=x0,y=y0取向前差分:
d y d x = y n + 1 − y n x n + 1 − x n → y n + 1 = y n + ( x n + 1 − x n ) ⋅ f ( x n , y n ) \frac{dy}{dx}=\frac{y_{n+1}-y_{n}}{x_{n+1}-x_{n}} \rightarrow y_{n+1}=y_{n}+(x_{n+1}-x_{n}){\cdot}f(x_{n},y_{n}) dxdy=xn+1−xnyn+1−yn→yn+1=yn+(xn+1−xn)⋅f(xn,yn)由此可以得到一个递推算法:
{ x = x 0 y = y 0 ⇒ f ( x 0 , y 0 ) ⇒ y 1 , { x = x 1 y = y 1 ⇒ f ( x 1 , y 1 ) ⇒ y 2 ⋯ ⇒ y n \begin{cases}x=x_0\\y=y_0\\\end{cases} \Rightarrow f(x_0,y_0) \Rightarrow y_1,\begin{cases}x=x_1\\y=y_1\\\end{cases} \Rightarrow f(x_1,y_1) \Rightarrow y_2 \dots \Rightarrow y_n {x=x0y=y0⇒f(x0,y0)⇒y1,{x=x1y=y1⇒f(x1,y1)⇒y2⋯⇒yn
(4) FORTRAN求解微分方程
program main
implicit none
real,parameter :: k0 = 100.0, E = 1e5, R = 8.314, C = 1.0,&
&n = 1.0, V = 1.0, Q = 2e7,A = 6.0,&
&rou = 1.0, cv = 1e3 !常量
real :: T0 = 300.0
real :: h = 5.0
real :: tao = 0, T = 300.0 !设定初值
real :: taoTemp, TTemp !临时交换变量
integer :: i
real :: deltaTao = 0.05 !差分步长
real,parameter :: taoMax = 100 !研究前100秒内的变化
real,parameter :: TMax = 1500 !规定零维燃烧的最高温度不超过1500K
open(unit = 1,file = 'e:/tao.csv') !将时间坐标读入文件
do i = 0,int(taoMax/deltaTao)
write(1,"(f6.2)") tao
tao = tao + deltaTao
end do
close(1)
tao = 0
open(unit = 2,file = 'e:/differentInitialTemperatureT.csv') !将温度坐标读入文件
!研究初始温度T0 = 300、400、500、600、700、800、900、1000、1100K时候的变化
do while(T0 <= 1100)
do while(T <= 1500) !T0温度下对应的变化
TTemp = T + (deltaTao*k0*exp(-E/R/T)*C**n*V*Q-A*h*(tao-T0)) / (V*rou*cv)
taoTemp = tao + deltaTao
tao = taoTemp
T = TTemp
write(2,"(f7.2)") T
end do
write(2,*) -273.15 !用来界定不同初温的数据
tao = 0
T = T0
T0 = T0 + 100.0
end do
close(2)
tao = 0
T0 = 800.0
T = 800.0
h = 1.0
open(unit = 3,file = 'e:/differentHeatDissipationCoefficientT.csv') !将温度坐标读入文件
!研究散热系数h = 1.0、2.0...10.0时候的变化
do while(h <= 10.0)
do while(T <= 1500) !每一个h下对应的温度变化
TTemp = T + (deltaTao*k0*exp(-E/R/T)*C**n*V*Q-A*h*(tao-T0)) / (V*rou*cv)
taoTemp = tao + deltaTao
tao = taoTemp
T = TTemp
write(3,"(f7.2)") T
end do
write(3,*) -273.15 !用来界定不同除温的数据
tao = 0
h = h + 1.0
T = 800.0
end do
close(3)
end program main
运行后生成三个.csv文件:
(5) MATLAB数据可视化
clc;clear;
%读取.csv文件中的数据
fileTao = 'e:\tao.csv';
fileDIT = 'e:\differentInitialTemperatureT.csv';
fileDDC = 'e:\differentHeatDissipationCoefficientT.csv';
tao = csvread(fileTao);
DIT_T = csvread(fileDIT);
DDC_T = csvread(fileDDC);
DITT = cell(1,9);
DITTao = cell(1,9);
DDCT = cell(1,10);
DDCTao = cell(1,10);
start = 1;
count = 1;
for i = 1:length(DIT_T)
if(DIT_T(i) == -273.15)
DITT{1,count} = DIT_T(start:i-1);
DITTao{1,count} = tao(1:length(DITT{1,count}));
count = count + 1;
start = i + 1;
end
end
start = 1;
count = 1;
for i = 1:length(DDC_T)
if(DDC_T(i) == -273.15)
DDCT{1,count} = DDC_T(start:i-1);
DDCTao{1,count} = tao(1:length(DDCT{1,count}));
count = count + 1;
start = i + 1;
end
end
%数据可视化
figure;
hold on;
grid on;
axis tight;
xlabel('time/s');
ylabel('temperature/K');
title('不同初始温度下零维燃烧温度变化曲线');
for i = 1:9
t = 0:0.005:max(DITTao{1,i});
u = interp1(DITTao{1,i},DITT{1,i},t,'spline');
plot(t,u);
end
legend('T_{0} = 300K','T_{0} = 400K','T_{0} = 500K','T_{0} = 600K','T_{0} = 700K','T_{0} = 800K',...
'T_{0} = 900K','T_{0} = 1000K','T_{0} = 1100K','location','southeast');
figure;
hold on;
grid on;
axis tight;
xlabel('time/s');
ylabel('temperature/K');
title('不同散热系数下零维燃烧温度变化曲线');
for i = 1:10
t = 0:0.005:max(DDCTao{1,i});
u = interp1(DDCTao{1,i},DDCT{1,i},t,'spline');
plot(t,u);
end
legend('h = 1.0 W/(m^{2}·K)','h = 2.0 W/(m^{2}·K)','h = 3.0 W/(m^{2}·K)','h = 4.0 W/(m^{2}·K)',...
'h = 5.0 W/(m^{2}·K)','h = 6.0 W/(m^{2}·K)','h = 7.0 W/(m^{2}·K)','h = 8.0 W/(m^{2}·K)',...
'h = 9.0 W/(m^{2}·K)','h = 10.0 W/(m^{2}·K)','location','southeast');
参 考 资 料 来 源 参考资料来源 参考资料来源
- 《燃烧学(第2版)》.徐通模.惠世恩.机械工业出版社.
源 码 作 者 : A i d e n L e e 源码作者:Aiden\ Lee 源码作者:Aiden Lee
博 客 创 作 : A i d e n L e e 博客创作:Aiden\ Lee 博客创作:Aiden Lee
特别声明:文章仅供学习参考,转载请注明出处,严禁盗用!