零维系统热自燃问题求解

wibibaboo

已于 2022-03-02 12:37:05 修改

阅读量1k

点赞数 1

分类专栏：杂货铺文章标签： matlab fortran 编程语言

于 2021-01-15 13:33:10 首次发布

本文链接：https://blog.csdn.net/dasdkhaksd/article/details/112657982

版权

杂货铺专栏收录该内容

10 篇文章 0 订阅

订阅专栏

零维系统热自燃问题求解

(1) 物理问题
(2) 求解目标
(3) Euler法求解微分方程
(4) FORTRAN求解微分方程
(5) MATLAB数据可视化

(1) 物理问题

$1 . 热平衡方程$

$V{\rho}c_{v}\frac{dT}{d\tau}=Q_{1}-Q_{2}$

产热 $Q_{1}$ 为： $Q_{1}=k_{0}e^{-\frac{E}{RT}}{\cdot}C^{n}{\cdot}V{\cdot}Q$
散热 $Q_{2}$ 为： $Q_{2}=h{\cdot}A{\cdot}(T-T_{0})$
联立解得： $\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}}$

这个方程手工求解很困难，但是用计算机编程的方法求解则很容易。

取初值如下：

$k_{0}$ ——频率因子，初值取100.0
$E$ ——活化能，初值取1.E5
$R$ ——气体常数，初值取8.314
$C$ ——可燃混合物中反应物浓度，初值取1.0
$n$ ——反应级数，初值取1.0
$V$ ——容器体积，初值取1.0
$Q$ ——可燃混合物的燃烧热，初值取2.E7
$T$ ——容器内可燃混合物温度，初值取与 $T_{0}$ 相等的数值。
$h$ ——散热系数，初值取5.0
A——容器壁散热面积，初值取6.0
$T_{0}$ ——容器壁温度，初值取800.0
$\rho$ ——混合气密度，初值取1.0
$c_{v}$ ——比热，初值取1.E3

初始条件为： $\tau=0,T=T_{0}$

(2) 求解目标

分别令系统初温 $T_{0}$ =300K、400K、500K、600K、700K、800K、900K、1000K、1100K（其他参数不变），获得不同初始温度下的系统升温曲线。
分别令散热系数 $h$ ＝1.0～10.0（ $T_{0}$ 保持800K），获得不同初始温度下的系统升温曲线。

(3) Euler法求解微分方程

设有一个微分方程：
$\frac{dy}{dx}=f(x,y(x))$

定解条件为：
$x=x_{0},y=y_{0}$

取向前差分：
$\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})$

由此可以得到一个递推算法：
$\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$

(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版)》.徐通模.惠世恩.机械工业出版社.

$源码作者：Aiden\ Lee$

$博客创作：Aiden\ Lee$

特别声明：文章仅供学习参考，转载请注明出处，严禁盗用！

wibibaboo

关注

1
点赞
踩
3

收藏

觉得还不错? 一键收藏
1
评论
零维系统热自燃问题求解

目录(1) 物理问题(2) 求解目标(3) Euler法求解微分方程(4) Fortran求解微分方程(5) Matlab数据可视化(1) 物理问题热平衡方程： VρcvdTdτ=Q1−Q2V{\rho}c_{v}\frac{dT}{d\tau}=Q_{1}-Q_{2}VρcvdτdT=Q1−Q2产热Q1Q_{1}Q1为： Q1=k0e−ERT⋅Cn⋅V⋅QQ_{1}=k_{0}e^{-\frac{E}{RT}}{\cdot}C^{n}{\cdot}V{\cdot}QQ1=k0e−R
复制链接

扫一扫