热动力数据MATLAB代码分享

最新推荐文章于 2024-06-06 20:37:20 发布

热力学负一定理

最新推荐文章于 2024-06-06 20:37:20 发布

阅读量2.9k

点赞数 5

分类专栏：实验数据分析文章标签： matlab 安全恰饭

本文链接：https://blog.csdn.net/qq_33980829/article/details/107189642

版权

实验数据分析专栏收录该内容

6 篇文章 2 订阅

订阅专栏

DSC 数据拟合和转化率计算

本文介绍
补充说明
热动力学介绍
总结

本文介绍

本文主要用于记录实验数据处理过程采用matlab进行数据处理分析
本文主要针基础过程进行代码分享，代码以解决问题为主没有过多注释和完善
代码极其粗野请紧张的往下看
代码主要涉及基础的拟合方法，作图格式设置方法
MATLAB的基础交互过程

补充说明

代码 jisuan_1 函数是里面注释部分，运行时额外粘贴一下。（2021年3月3日）

热动力学介绍

你好！建议[点击进行查询][1]

拟合采用kissinger方程

特殊形式
$\frac{\partial \alpha }{\partial t}=A \exp (-E/RT) (1-\alpha)^{n}$
通用形式
$\frac{\partial \alpha }{\partial t}=A \exp (-E/RT)f(\alpha )$
通过推导最后建立出反应放热峰值温度与升温速率的关系

$\ln(\frac{\beta _{i}}{T^{2}_{pi}})=\ln(\frac{A_{k}R}{E_{K}})- \frac{E_{k}}{RT^{2}_{pi}} {[i=1, 2,3 ..6]}$
根据公式可以知道最后反应归结为升温速率与峰值温度的关系。只需要将数据代入进行直线拟合就能算出表观活化能 $E_{K}$ 和反应频率 $A$ ,

拟合过程代码

%% 对dsc数据进行拟合  k方法
% 输入放热峰值
% 输入温度梯度
% 结果保存
% 2020/7/4 
% 版本 0.1 尚帝
%% 清理变量空间
clc
clear  
close all
%% 数据输入与预处理
x_1=[10 15 20 25 30];  % 温度梯度变化 
y_x=[687.64 703.4 716.21 726.87 730.55]; %样本数据的峰值温度 T_P
y_y=[740.10  757.52  780.78  793.49  797.40]; %样本数据的峰值温度  T_P
y_z=[725.51 748.39 759.47 770.79 774.49];  %样本数据的峰值温度  T_P
%% 拟合过程
[E_K_z, A_K_z,fitresult_z, gof_z]=jisuan_1(x_1,y_z,'Z的拟合结果');
[E_K_x, A_K_x,fitresult_x, gof_x]=jisuan_1(x_1,y_x,'X的拟合结果');
[E_K_y, A_K_y,fitresult_y, gof_y]=jisuan_1(x_1,y_y,'Y的拟合结果');
A_K_ALL=[A_K_x,A_K_y,A_K_z];
E_K_ALL=[E_K_x,E_K_y,E_K_z];
%% 结果汇总
rsquare_all=[gof_x.rsquare,gof_y.rsquare,gof_z.rsquare];
all_num=[E_K_ALL' A_K_ALL' rsquare_all'];
all_num=array2table(all_num, 'VariableNames',{'E_K','A','R^2'});
writetable(all_num,'拟合结果.xlsx')
%% 图像处理部分
function [E_K, A_K, fitresult, gof]=jisuan_1(x_1,y_1,name)
x=[x_1;y_1]';
gradient(x);
y1=log(x(1:5,1)./x(1:5,2).^2);
x1=1./x(1:5,2);
[xData, yData] = prepareCurveData( x1, y1 );
% Set up fittype and options.
ft = fittype( 'poly1' );
[fitresult, gof] = fit( xData, yData, ft );
figure( 'Name',name );
%% 图像形式处理
% 创建 axes
axes1 = axes('Parent',gcf);
hold(axes1,'on');
set(axes1,'FontName','Times New Roman','FontSize',13);
h = plot( fitresult, xData, yData );
legend( h, 'y1 vs. x1', '拟合方程', 'Location', 'NorthEast' );
% Label axes
xlabel({'x1'},'FontWeight','bold','FontSize',13,'FontName','Times New Roman');
ylabel({'y1'},'FontWeight','bold','FontSize',13,'FontName','Times New Roman');

title({name},'FontWeight','bold','FontSize',15,'FontName','宋体');

grid off
saveas(gcf,strcat(name,".bmp")); 
E_K=-fitresult.p1*8.314/10^3;
A_K=fitresult.p2+log(-fitresult.p1);
disp(strcat("反应的E_K=",num2str(E_K)))
disp(strcat("反应的A_K=",num2str(A_K)))
disp("-----------------------------------------")
end

自定义拟合函数

%% 图像处理部分
function [E_K, A_K, fitresult, gof]=jisuan_1(x_1,y_1,name)
x=[x_1;y_1]';
gradient(x);
y1=log(x(1:5,1)./x(1:5,2).^2);
x1=1./x(1:5,2);
[xData, yData] = prepareCurveData( x1, y1 );
% Set up fittype and options.
ft = fittype( 'poly1' );
[fitresult, gof] = fit( xData, yData, ft );
figure( 'Name',name );
%% 图像形式处理
% 创建 axes
axes1 = axes('Parent',gcf);
hold(axes1,'on');
set(axes1,'FontName','Times New Roman','FontSize',13);
h = plot( fitresult, xData, yData );
legend( h, 'y1 vs. x1', '拟合方程', 'Location', 'NorthEast' );
% Label axes
xlabel({'x1'},'FontWeight','bold','FontSize',13,'FontName','Times New Roman');
ylabel({'y1'},'FontWeight','bold','FontSize',13,'FontName','Times New Roman');

title({name},'FontWeight','bold','FontSize',15,'FontName','宋体');

grid off
saveas(gcf,strcat(name,".bmp")); 
E_K=-fitresult.p1*8.314/10^3;
A_K=fitresult.p2+log(-fitresult.p1);
disp(strcat("反应的E_K=",num2str(E_K)))
disp(strcat("反应的A_K=",num2str(A_K)))
disp("-----------------------------------------")
end

拟合结果展示

通过计算，输出反应结果和拟合的r^2大小

拟合结果展示

反应转率介绍

通过查阅文献，反应速率是通过定义为 [DSC反应过程下TA曲线变化对应热流曲线下的面积][2]可以定义为下面的公式
$\alpha _{i}=\int_{T1}^{Ti}\frac{\mathrm{d}\alpha }{\mathrm{d}T} \: \: [i=1,2...n] [0<\alpha _{i}=<1]$

公式中的 $T_{i}$ 表示启示反应温度点， $T_{n}$ 表示反应结束温度，可以通过热重曲线确定反应的开始和结束。反应转化率表示为随温度升高反应产物的积累量。

反应产物转化计算程序

导入DSC数据

%% 时间 2020/7/3 
% 输入DSC数据，求取反应转化过程
% 0.1版
% 对转化率进行拟合 
% 输入dsc数据 输出转化率与温度的变化 输出反应时间变化
% 尚帝
%% 清理空间
clc
clear 
close all
%% 文件读取设置
name_U='Z10'; % 计算的名字
%% 读取文件内容
opts = spreadsheetImportOptions("NumVariables", 9);
% 指定工作表和范围
opts.Sheet = name_U;
opts.DataRange = "A1:I11755";
% 指定列名称和类型
opts.VariableNames = ["VarName1", "VarName2", "VarName3", "VarName4", "VarName5", "VarName6", "VarName7", "VarName8", "VarName9"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double", "double"];
% 导入数据
x1 = readtable("C:\Users\Administrator\Desktop\dsc数据\运算数据\Z转化率数据.xlsx", opts, "UseExcel", false);
%% 转换为输出类型
x1 = table2array(x1);
%% 清除临时变量
clear opts
%% 数据预处理
p_h=x1(:,4); %热流密度
p_t=x1(:,2); %温度曲线
figure 
plot(p_t,p_h);
% 输入起始与结束点
[x_s1,y_s1]=ginput(1);
hold on
plot(x_s1,y_s1,'o')
[x_s2,y_s2]=ginput(1);
hold on
plot(x_s2,y_s2,'o')
plot([x_s1,x_s2],[y_s1,y_s2])
hold off
% 判断点选取的合理性
if y_s1>max(p_h)&&y_s1<min(p_h)
    disp('选择的点不存在于线上')
    x_s1=655;
    x_s2=737.6;
    y_s1= p_h(max(p_t<=x_s1&p_t>=x_s1-0.1));
    y_s2= p_h(max(p_t<=x_s2&p_t>=x_s2-0.1)) ;
end
x=[x_s1;x_s2];
y=[y_s1;y_s2];
% 截断后的数据
p_t1=p_t(p_t>=x_s1 & p_t<= x_s2); 
p_h1=p_h(p_t>=x_s1 & p_t<= x_s2);
% 求取两点的直线
[xData, yData] = prepareCurveData( x, y );
ft = fittype( 'poly1' );
[fitresult, gof] = fit( xData, yData, ft );
% 截断数据减去直线
figure
plot(fitresult,x,y)
hold on
plot(p_t1,p_h1);
[fitresult1, gof1] = fit_p_1(p_t1, p_h1);
hold off
figure
p_h2=p_h1-fitresult(p_t1);
plot(p_t1,p_h2);
%% 计算转化率
p_h2_d=diff(p_h2); %对热流算差分
p_t1_d=diff(p_t1); %对温度算差分
a_m=zeros(size(p_t1_d,1),1);
for i=1:size(p_t1_d,1)
    if i==1
        a_m(i)=(p_h2(i)+p_h2(i+1))*p_t1_d(i)*0.5;
    else
        a_m(i)=(p_h2(i)+p_h2(i+1))*p_t1_d(i)*0.5+a_m(i-1);
    end
end
a_m_gui=mapminmax(a_m',0,1); % 归一化
figure
p_t1_end=p_t1(1:end-2);
p_t1_end_1=p_t1(1:end-1);
plot(p_t1(1:end-1),a_m_gui);
xlabel('温度')
ylabel('转化率')
a_m_u=diff(a_m_gui);
 a_m_u_1=smooth(a_m_u);
 figure
 plot(a_m_u_1)
 [fitresult2, gof2] = fit_p_1(p_t1_end, a_m_u_1);
web("www.baidu.com")
%% 保存数据 
% dsc曲线
name_z10=a_m_u_1;
save(strcat(name_U,"_dsc_1.mat") ,'name_z10') 
% 转化率
name_z10_p=a_m_gui;
save(strcat(name_U,"_tr.mat") ,'name_z10_p') 
% 温度范围
name_z10_t=p_t1_end_1;
save(strcat(name_U,"_tr_t.mat") ,'name_z10_t')

代码运行过程展示

交互过程选取反应起始和结束范围

交互过程展示
在这里插入图片描述

不同升温速率梯度下的结果汇总

在这里插入图片描述

根据反应区间范围计算反应转化率代码

%% 清理空间
clc
clear 
close all
a=1;
%% 数据读取
name_p='y粒径的反应';
for i=10:5:30
    name=strcat("y",num2str(i),'_tr.mat');
    name1=strcat("y",num2str(i),'_tr_t.mat');
    load(name)
    load(name1)
    y=eval(strcat('name_y',num2str(i),"_p"));
    x=eval(strcat('name_y',num2str(i),"_t"));
    f=plot(x,y,'DisplayName',strcat('反应速率',num2str(i)));
    legend('show');
    hold on
end
xlabel({'t'},'FontWeight','bold','FontSize',13,'FontName','Times New Roman');
ylabel({'a'},'FontWeight','bold','FontSize',13,'FontName','Times New Roman');
title({name_p},'FontWeight','bold','FontSize',15,'FontName','宋体');
saveas(gcf,strcat(name_p,".bmp"))
hold off

name_p='x粒径的反应';
for i=10:5:30
    name=strcat("x",num2str(i),'_tr.mat');
    name1=strcat("x",num2str(i),'_tr_t.mat');
    load(name)
    load(name1)
    y=eval(strcat('name_x',num2str(i),"_p"));
    x=eval(strcat('name_x',num2str(i),"_t"));
    f=plot(x,y,'DisplayName',strcat('反应速率',num2str(i)));
    legend('show');
    hold on
end
xlabel({'t'},'FontWeight','bold','FontSize',13,'FontName','Times New Roman');
ylabel({'a'},'FontWeight','bold','FontSize',13,'FontName','Times New Roman');
title({name_p},'FontWeight','bold','FontSize',15,'FontName','宋体');
saveas(gcf,strcat(name,".bmp"))
hold off
name_p='x粒径的反应';
for i=10:5:30
    name=strcat("z",num2str(i),'_tr.mat');
    name1=strcat("z",num2str(i),'_tr_t.mat');
    load(name)
    load(name1)
    y=eval(strcat('name_z',num2str(i),"_p"));
    x=eval(strcat('name_z',num2str(i),"_t"));
    f=plot(x,y,'DisplayName',strcat('反应速率',num2str(i)));
    legend('show');
    hold on
end
xlabel({'t'},'FontWeight','bold','FontSize',13,'FontName','Times New Roman');
ylabel({'a'},'FontWeight','bold','FontSize',13,'FontName','Times New Roman');
title({name_p},'FontWeight','bold','FontSize',15,'FontName','宋体');
saveas(gcf,strcat(name_p,".bmp"))
hold off

总结

*咕咕咕咕本文主要为做动力学曲线分析
[1]: www.baidu.com
[2]:
[3]: /
[4]:

热力学负一定理

关注

5
点赞
踩
20

收藏

觉得还不错? 一键收藏
16
评论
热动力数据MATLAB代码分享

DSC 数据拟合和转化率计算本文介绍热动力学介绍拟合采用kissinger方程拟合过程代码拟合结果展示通过计算，输出反应结果和拟合的r^2大小反应转率介绍反应产物转化计算程序导入DSC数据过程展示交互过程选取反应起始和结束范围不同升温速率梯度下的结果汇总根据反应区间范围计算反应转化率代码总结本文介绍本文主要用于记录实验数据处理过程采用matlab进行数据处理分析本文主要针基础过程进行代码分享，代码以解决问题为主没有过多注释和完善代码极其粗野请紧张的往下看代码主要涉及基础的拟合方法，作图格式设
复制链接

扫一扫