本文介绍

• 本文主要用于记录实验数据处理过程 采用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_{pi}^2})=\ln (\frac{A_{k}R}{E_K})-\frac{E_k}{R{T}_{pi}^2}[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:UsersAdministratorDesktopdsc数据运算数据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]:

最后

以上就是优美翅膀为你收集整理的热动力数据MATLAB代码分享本文介绍补充说明热动力学介绍总结的全部内容，希望文章能够帮你解决热动力数据MATLAB代码分享本文介绍补充说明热动力学介绍总结所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。