Matlab實驗資料處理程式大彙總

卷卷怪發表於2020-12-26
Matlab

目錄

(1) 換熱器綜合實驗

(1.1) 程式碼

clc;clear;

heatExchangeArea = 0.14; %換熱面積。
flowArea = (pi / 4) * ((16 - 1.5 * 2) * 0.001)^2; %內管流體橫截面面積。

%順流。
downStreamTi1 = [60.2 59.4 60.0 60.1 60.0]; %熱流體進口溫度（℃）。
downStreamTo1 = [54.0 52.7 52.4 51.6 50.7]; %熱流體出口溫度（℃）。
downStreamTi2 = [14.5 13.8 13.3 13.3 13.1]; %冷流體進口溫度（℃）。
downStreamTo2 = [36.4 35.2 34.6 34.0 33.1]; %冷流體出口溫度（℃）。
downStreamQ1 = [548.4 503.0 453.4 403.8 352.8]; %熱流體流量（L/h）。
downStreamQ2 = [160 160 160 160 160]; %冷流體流量（L/h）。
downStreamTmax = zeros(1,5);
downStreamTmin = zeros(1,5);
downStreamTm = zeros(1,5); %對數平均溫差（℃）。
downStreamDensity1 = [985.38 985.86 985.79 985.95 986.18]; %熱流體密度（kg/m^3）。
downStreamDensity2 = [997.00 997.24 997.37 997.44 997.57]; %冷流體密度（kg/m^3）。
downStreamHeatCapacity1 = [4.1776 4.1770 4.1771 4.1769 4.1767]; %熱流體比熱容（kJ/(kg*K)）。
downStreamHeatCapacity2 = [4.1781 4.1790 4.1794 4.1797 4.1802]; %冷流體比熱容（kJ/(kg*K)）。
downStreamHeat1 = zeros(1,5); %熱流體換熱量（kW）。
downStreamHeat2 = zeros(1,5); %冷流體換熱量（kW）。
downStreamHeatErr = zeros(1,5); %熱平衡偏差（%）。
downStreamK = zeros(1,5); %傳熱係數。
for i = 1:5
    %求對數平均溫差。
    downStreamTmax(i) = max(abs(downStreamTi1(i) - downStreamTi2(i)),abs(downStreamTo1(i) - downStreamTo2(i)));
    downStreamTmin(i) = min(abs(downStreamTi1(i) - downStreamTi2(i)),abs(downStreamTo1(i) - downStreamTo2(i)));
    downStreamTm(i) = (downStreamTmax(i) - downStreamTmin(i)) / (log(downStreamTmax(i) / downStreamTmin(i)));
    %求熱流體換熱量。
    downStreamHeat1(i) = downStreamQ1(i) * 0.001/3600 * downStreamDensity1(i) * downStreamHeatCapacity1(i) * ...
    (downStreamTi1(i) - downStreamTo1(i));
    %求冷流體換熱量。
    downStreamHeat2(i) = downStreamQ2(i) * 0.001/3600 * downStreamDensity2(i) * downStreamHeatCapacity2(i) * ...
    (downStreamTo2(i) - downStreamTi2(i));
    %求熱平衡偏差。
    downStreamHeatErr(i) = 2 * (downStreamHeat1(i) - downStreamHeat2(i)) / (downStreamHeat1(i) + ...
    downStreamHeat2(i)) * 100;
    %求傳熱係數。
    downStreamK(i) = (downStreamHeat1(i) + downStreamHeat2(i)) / 2 / heatExchangeArea / downStreamTm(i);
end
disp("順流下各工況=>");
fprintf("\t1.對數平均溫差：(℃)\n");
disp(downStreamTm);
fprintf("\t2.熱流體換熱量：（kW）\n");
disp(downStreamHeat1);
fprintf("\t3.冷流體換熱量：（kW）\n");
disp(downStreamHeat2);
fprintf("\t4.熱平衡偏差：（%%）\n");
disp(downStreamHeatErr);
fprintf("\t5.傳熱係數：\n");
disp(downStreamK);

%逆流。
upStreamTi1 = [59.6 59.8 60.3 60.2 60.4]; %熱流體進口溫度（℃）。
upStreamTo1 = [53.2 52.8 52.5 51.5 50.7]; %熱流體出口溫度（℃）。
upStreamTi2 = [13.0 13.0 13.0 13.0 13.0]; %冷流體進口溫度（℃）。
upStreamTo2 = [35.8 35.3 35.1 34.5 33.8]; %冷流體出口溫度（℃）。
upStreamQ1 = [554.3 503.9 453.5 402.6 353.2]; %熱流體流量（L/h）。
upStreamQ2 = [160 160 160 160 160]; %冷流體流量（L/h）。
upStreamTmax = zeros(1,5);
upStreamTmin = zeros(1,5);
upStreamTm = zeros(1,5); %對數平均溫差（℃）。
upStreamDensity1 = [985.70 985.74 985.70 985.95 986.09]; %熱流體密度（kg/m^3）。
upStreamDensity2 = [997.26 997.32 997.35 997.42 997.50]; %冷流體密度（kg/m^3）。
upStreamHeatCapacity1 = [4.1772 4.1772 4.1772 4.1769 4.1768]; %熱流體比熱容（kJ/(kg*K)）。
upStreamHeatCapacity2 = [4.1790 4.1793 4.1794 4.1796 4.1799]; %冷流體比熱容（kJ/(kg*K)）。
upStreamHeat1 = zeros(1,5); %熱流體換熱量（kW）。
upStreamHeat2 = zeros(1,5); %冷流體換熱量（kW）。
upStreamHeatErr = zeros(1,5); %熱平衡偏差（%）。
upStreamK = zeros(1,5); %傳熱係數。
for i = 1:5
    %求對數平均溫差。
    upStreamTmax(i) = max(abs(upStreamTi1(i) - upStreamTo2(i)),abs(upStreamTo1(i) - upStreamTi2(i)));
    upStreamTmin(i) = min(abs(upStreamTi1(i) - upStreamTo2(i)),abs(upStreamTo1(i) - upStreamTi2(i)));
    upStreamTm(i) = (upStreamTmax(i) - upStreamTmin(i)) / (log(upStreamTmax(i) / upStreamTmin(i)));
    %求熱流體換熱量。
    upStreamHeat1(i) = upStreamQ1(i) * 0.001/3600 * upStreamDensity1(i) * upStreamHeatCapacity1(i) * ...
    (upStreamTi1(i) - upStreamTo1(i));
    %求冷流體換熱量。
    upStreamHeat2(i) = upStreamQ2(i) * 0.001/3600 * upStreamDensity2(i) * upStreamHeatCapacity2(i) * ...
    (upStreamTo2(i) - upStreamTi2(i));
    %求熱平衡偏差。
    upStreamHeatErr(i) = 2 * (upStreamHeat1(i) - upStreamHeat2(i)) / (upStreamHeat1(i) + ...
    upStreamHeat2(i)) * 100;
    %求傳熱係數。
    upStreamK(i) = (upStreamHeat1(i) + upStreamHeat2(i)) / 2 / heatExchangeArea / upStreamTm(i);
end
disp("逆流下各工況=>");
fprintf("\t1.對數平均溫差：(℃)\n");
disp(upStreamTm);
fprintf("\t2.熱流體換熱量：（kW）\n");
disp(upStreamHeat1);
fprintf("\t3.冷流體換熱量：（kW）\n");
disp(upStreamHeat2);
fprintf("\t4.熱平衡偏差：（%%）\n");
disp(upStreamHeatErr);
fprintf("\t5.傳熱係數：\n");
disp(upStreamK);

%用最小二乘法擬合資料求實驗關聯式。
disp("實驗關聯式擬合=>");
downStreamVelocity = downStreamQ1 * 0.001 / 3600 / flowArea;
upStreamVelocity = upStreamQ1 * 0.001 / 3600 / flowArea;
downStreamP = polyfit(log(downStreamVelocity),log(downStreamK),1); %最小二乘擬合算出順流的n和lnc。
fprintf("\t1.順流實驗關聯式係數為：\n");
fprintf("\t\tc = %6.3f\tn = %6.3f\n",exp(downStreamP(2)),downStreamP(1));
upStreamP = polyfit(log(upStreamVelocity),log(upStreamK),1); %最小二乘擬合算出逆流的n和lnc。
fprintf("\t2.逆流實驗關聯式係數為：\n");
fprintf("\t\tc = %6.3f\tn = %6.3f\n",exp(upStreamP(2)),upStreamP(1));
%資料視覺化。
figure;
hold on;
axis tight;
grid on;
xlabel('u/(m·s^{-1})');
ylabel('k');
title('套管式換熱器傳熱係數與管內流體流速的關係');
scatter(downStreamVelocity,downStreamK,'rp');
h1 = fplot(@(x)exp(downStreamP(2)).*x.^downStreamP(1),[0.6,1.2]);
set(h1,'linestyle','-.');
scatter(upStreamVelocity,upStreamK,'b*');
h2 = fplot(@(x)exp(upStreamP(2)).*x.^upStreamP(1),[0.6,1.2]);
set(h2,'linestyle','--');
legend('順流工況資料點','順流工況擬合曲線','逆流工況資料點','逆流工況擬合曲線','location','northwest');
downStreamCurveName = strcat('k=',num2str(exp(downStreamP(2)),3),'u','^','{',num2str(downStreamP(1),3),'}');
upStreamCurveName = strcat('k=',num2str(exp(upStreamP(2)),3),'u','^','{',num2str(upStreamP(1),3),'}');
text(0.7,0.92,downStreamCurveName);
text(0.9,0.89,upStreamCurveName);

(1.2) 執行結果

順流下各工況=>
	1.對數平均溫差：(℃)
   29.4488   29.3409   29.9623   29.8573   29.8943
	2.熱流體換熱量：（kW）
    3.8879    3.8550    3.9414    3.9264    3.7540
	3.冷流體換熱量：（kW）
    4.0545    3.9637    3.9461    3.8355    3.7067
	4.熱平衡偏差：（%）
   -4.1943   -2.7821   -0.1187    2.3421    1.2688
	5.傳熱係數：
    0.9632    0.9517    0.9402    0.9284    0.8913
逆流下各工況=>
	1.對數平均溫差：(℃)
   31.2869   31.5338   31.8162   31.6701   31.8281
	2.熱流體換熱量：（kW）
    4.0574    4.0345    4.0458    4.0068    3.9197
	3.冷流體換熱量：（kW）
    4.2231    4.1310    4.0942    3.9835    3.8544
	4.熱平衡偏差：（%）
   -4.0015   -2.3653   -1.1909    0.5827    1.6787
	5.傳熱係數：
    0.9452    0.9248    0.9137    0.9011    0.8723
實驗關聯式擬合=>
	1.順流實驗關聯式係數為：
		c =  0.945	n =  0.166
	2.逆流實驗關聯式係數為：
		c =  0.921	n =  0.166

(2) 平板邊界層內的流速分佈實驗

(2.1) 程式碼

clc;clear;

%常量的定義。
Rg=287;
g=9.8;
rou_water=1000;

%對已知實驗引數的輸入。
Pa=input('請輸入大氣壓Pa：(單位Pa)');
t=input('請輸入溫度t：(單位℃)');
miu=input('請輸入流體運動粘性係數：(單位1e-7m^2/s)');
K=input('請輸入風洞係數K：');
C=input('請輸入微壓計所選的檔位係數C：');
delta_h_inf=input('請輸入風洞進口處真空度：(單位mmH2O)');
Pj0=input('請輸入流場靜壓Pj：(單位毫米酒精柱)');

%第一截面實驗資料輸入。
x10=input('請輸入第一截面距前緣的距離x：(單位mm)');
y10=input('請輸入第一截面邊界層內距離的資料(用[]括起來，元素之間用空格分隔)：(單位mm)');
p10_micromanometer=input('請輸入第一截面用微壓計測得的讀數資料(用[]括起來，元素之間用空格分隔)：(單位毫米酒精柱)');
p10_pressure_difference_transmitter=input('請輸入第一截面用壓差變送器測得的讀數資料(用[]括起來，元素之間用空格分隔)：(單位Pa)');

%第二截面實驗資料輸入。
x20=input('請輸入第二截面距前緣的距離x：(單位mm)');
y20=input('請輸入第二截面邊界層內距離的資料(用[]括起來，元素之間用空格分隔)：(單位mm)');
p20_micromanometer=input('請輸入第二截面用微壓計測得的讀數資料(用[]括起來，元素之間用空格分隔)：(單位毫米酒精柱)');
p20_pressure_difference_transmitter=input('請輸入第二截面用壓差變送器測得的讀數資料(用[]括起來，元素之間用空格分隔)：(單位Pa)');

%實驗過程中通過計算出來的、不變的資料。
rou_air=Pa/(Rg*(273.15+t));
Pj=Pa-C*Pj0*rou_water*g*0.001;
x1=x10*0.001;
y1=y10*0.001;
x2=x20*0.001;
y2=y20*0.001;
v_inf=sqrt((0.001*(2*K*g*rou_water*delta_h_inf))/rou_air);
fprintf('=>風洞進口流速為：%.3f m/s\n',v_inf);
p1_micromanometer=Pa-p10_micromanometer*C*rou_water*g*0.001;
p1_pressure_difference_transmitter=Pa-p10_pressure_difference_transmitter;
p2_micromanometer=Pa-p20_micromanometer*C*rou_water*g*0.001;
p2_pressure_difference_transmitter=Pa-p20_pressure_difference_transmitter;

%對第一截面資料進行處理
fprintf('>>>>>>>>>>**第一截面資料計算**<<<<<<<<<<\n');
Re1=(v_inf*x1)/(miu*1e-7);
fprintf('=>雷諾數 Re= %.3f\n',Re1);
fprintf('=>微壓計實際讀數計算出來的絕對壓強：(單位Pa)\n');
disp(p1_micromanometer);
u1_micromanometer=zeros(1,length(p1_micromanometer));
for i1=1:length(p1_micromanometer)
    u1_micromanometer(i1)=sqrt(2*(p1_micromanometer(i1)-Pj)/rou_air);
end
fprintf('=>用微壓計測量的資料計算得到的邊界層內的流速：（單位m/s）\n');
disp(u1_micromanometer);
u1_pressure_difference_transmitter=zeros(1,length(p1_pressure_difference_transmitter));
for j1=1:length(p1_pressure_difference_transmitter)
    u1_pressure_difference_transmitter(j1)=sqrt(2*(p1_pressure_difference_transmitter(j1)-Pj)/rou_air);
end
fprintf('=>用壓差變送器測量的資料計算得到的邊界層內的流速：（單位m/s）\n');
disp(u1_pressure_difference_transmitter);
velocity_ratio1_micromanometer=u1_micromanometer/v_inf;
fprintf('=>用微壓計測量的資料計算得到的邊界層內的速度比：\n');
disp(velocity_ratio1_micromanometer);
velocity_ratio1_pressure_difference_transmitter=u1_pressure_difference_transmitter/v_inf;
fprintf('=>用壓差變送器測量的資料計算得到的邊界層內的速度比：\n');
disp(velocity_ratio1_pressure_difference_transmitter);
subplot(2,1,1);
hold on;grid on;
plot(velocity_ratio1_micromanometer,y1,'b.','markersize',3);
plot(velocity_ratio1_pressure_difference_transmitter,y1,'g.','markersize',3);
t1_1=min(velocity_ratio1_micromanometer):0.001:max(velocity_ratio1_micromanometer);
u1_1=spline(velocity_ratio1_micromanometer,y1,t1_1);
delta1_micromanometer=spline(velocity_ratio1_micromanometer,y1,0.93);
fprintf('=>用微壓計測量得出的邊界層厚度為：%.3f mm\n',delta1_micromanometer*1000);
plot(t1_1,u1_1,'b-','linewidth',1);
t1_2=min(velocity_ratio1_pressure_difference_transmitter):0.001:max(velocity_ratio1_pressure_difference_transmitter);
u1_2=spline(velocity_ratio1_pressure_difference_transmitter,y1,t1_2);
delta1_pressure_difference_transmitter=spline(velocity_ratio1_pressure_difference_transmitter,y1,0.93);
fprintf('=>用壓差變送器測量得出的邊界層厚度為：%.3f mm\n',delta1_pressure_difference_transmitter*1000);
plot(t1_2,u1_2,'g-','linewidth',1);
y1_micromanometer=0:0.00001:delta1_pressure_difference_transmitter;
x1_micromanometer=((3*y1_micromanometer-(y1_micromanometer.^3)/...
    (delta1_pressure_difference_transmitter^2))/(2*delta1_pressure_difference_transmitter));
plot(x1_micromanometer,y1_micromanometer,'k--','linewidth',1);
y1_pressure_difference_transmitter=0:0.00001:delta1_pressure_difference_transmitter;
x1_pressure_difference_transmitter=(y1_pressure_difference_transmitter/delta1_pressure_difference_transmitter).^(1/7);
plot(x1_pressure_difference_transmitter,y1_pressure_difference_transmitter,'k-.','linewidth',1);
title('第一截面邊界層速度分佈曲線');
xlabel('速度比');
ylabel('離平板壁面距離y/m');
axis([0 1 0 10e-3]);
text(0.1,8e-3,'藍色線代表用微壓計測量得到的擬合曲線');
text(0.1,7e-3,'綠色線代表用壓差變送器測量得到的擬合曲線');
text(0.1,6e-3,'虛線代表層流理論曲線');
text(0.1,5e-3,'點畫線線代表紊流理論曲線');

%對第二截面資料進行處理。
fprintf('>>>>>>>>>>**第二截面資料計算**<<<<<<<<<<\n');
Re2=(v_inf*x2)/(miu*1e-7);
fprintf('=>雷諾數 Re= %.3f\n',Re2);
fprintf('=>微壓計實際讀數計算出來的絕對壓強：(單位Pa)\n');
disp(p2_micromanometer);
u2_micromanometer=zeros(1,length(p2_micromanometer));
for i2=1:length(p2_micromanometer)
    u2_micromanometer(i2)=sqrt(2*(p2_micromanometer(i2)-Pj)/rou_air);
end
fprintf('=>用微壓計測量的資料計算得到的邊界層內的流速：（單位m/s）\n');
disp(u2_micromanometer);
u2_pressure_difference_transmitter=zeros(1,length(p2_pressure_difference_transmitter));
for j2=1:length(p2_pressure_difference_transmitter)
    u2_pressure_difference_transmitter(j2)=sqrt(2*(p2_pressure_difference_transmitter(j2)-Pj)/rou_air);
end
fprintf('=>用壓差變送器測量的資料計算得到的邊界層內的流速：（單位m/s）\n');
disp(u2_pressure_difference_transmitter);
velocity_ratio2_micromanometer=u2_micromanometer/v_inf;
fprintf('=>用微壓計測量的資料計算得到的邊界層內的速度比：\n');
disp(velocity_ratio2_micromanometer);
velocity_ratio2_pressure_difference_transmitter=u2_pressure_difference_transmitter/v_inf;
fprintf('=>用壓差變送器測量的資料計算得到的邊界層內的速度比：\n');
disp(velocity_ratio2_pressure_difference_transmitter);
subplot(2,1,2);
hold on;grid on;
plot(velocity_ratio2_micromanometer,y2,'b.','markersize',3);
plot(velocity_ratio2_pressure_difference_transmitter,y2,'g.','markersize',3);
t2_1=min(velocity_ratio2_micromanometer):0.001:max(velocity_ratio2_micromanometer);
u2_1=spline(velocity_ratio2_micromanometer,y2,t2_1);
delta2_micromanometer=spline(velocity_ratio2_micromanometer,y2,0.93);
fprintf('=>用微壓計測量得出的邊界層厚度為：%.3f mm\n',delta2_micromanometer*1000);
plot(t2_1,u2_1,'b-','linewidth',1);
t2_2=min(velocity_ratio2_pressure_difference_transmitter):0.001:max(velocity_ratio2_pressure_difference_transmitter);
u2_2=spline(velocity_ratio2_pressure_difference_transmitter,y2,t2_2);
delta2_pressure_difference_transmitter=spline(velocity_ratio2_pressure_difference_transmitter,y2,0.93);
fprintf('=>用壓差變送器測量得出的邊界層厚度為：%.3f mm\n',delta2_pressure_difference_transmitter*1000);
plot(t2_2,u2_2,'g-','linewidth',1);
y2_micromanometer=0:0.00001:delta2_pressure_difference_transmitter;
x2_micromanometer=((3*y2_micromanometer-(y2_micromanometer.^3)/...
    (delta2_pressure_difference_transmitter^2)))/(2*delta2_pressure_difference_transmitter);
plot(x2_micromanometer,y2_micromanometer,'k--','linewidth',1);
y2_pressure_difference_transmitter=0:0.00001:delta2_pressure_difference_transmitter;
x2_pressure_difference_transmitter=(y2_pressure_difference_transmitter/delta2_pressure_difference_transmitter).^(1/7);
plot(x2_pressure_difference_transmitter,y2_pressure_difference_transmitter,'k-.','linewidth',1);
title('第二截面邊界層速度分佈曲線');
xlabel('速度比');
ylabel('離平板壁面距離y/m');
axis([0 1 0 10e-3]);
text(0.1,8e-3,'藍色線代表用微壓計測量得到的擬合曲線');
text(0.1,7e-3,'綠色線代表用壓差變送器測量得到的擬合曲線');
text(0.1,6e-3,'虛線代表層流理論曲線');
text(0.1,5e-3,'點畫線代表紊流曲線');

%流態分析。
fprintf('>>>>>>>>>>**流態分析**<<<<<<<<<<\n');
fprintf('一、對第一截面：\n');
fprintf('1.用雷諾數判定流態：\n');
if delta1_micromanometer>=3e6
    fprintf('=>（用微壓計測量所得的資料分析）流動為紊流\n');
else
    fprintf('=>（用微壓計測量所得的資料分析）流動為層流\n');
end
if delta1_pressure_difference_transmitter>=3e6
    fprintf('=>（用壓差變送器測量所得的資料分析）流動為紊流\n');
else
    fprintf('=>（用壓差變送器測量所得的資料分析）流動為層流\n');
end
fprintf('2.用邊界層厚度計算值與理論值比較判定流態：\n');
if abs(delta1_micromanometer-5*sqrt((miu*1e-7*x1)/v_inf))<0.1
    fprintf('=>(用微壓計測量所得的資料分析)流動為層流\n');
elseif abs(delta1_micromanometer-0.37*x1*((miu*1e-7)/(v_inf*x1))^0.2)<0.1
    fprintf('=>(用微壓計測量所得的資料分析)流動為紊流\n');
else
    fprintf('=>測量值與理論值相差較大，無法判定流態！\n');
end
if abs(delta1_pressure_difference_transmitter-5*sqrt((miu*1e-7*x1)/v_inf))<0.1
    fprintf('=>(用微壓計測量所得的資料分析)流動為層流\n');
elseif abs(delta1_pressure_difference_transmitter-0.37*x1*((miu*1e-7)/(v_inf*x1))^0.2)<0.1
    fprintf('=>(用微壓計測量所得的資料分析)流動為紊流\n');
else
    fprintf('=>測量值與理論值相差較大，無法判定流態！\n');
end
fprintf('二、對第二截面：\n');
fprintf('1.用雷諾數判定流態：\n');
if delta2_micromanometer>=3e6
    fprintf('=>（用微壓計測量所得的資料分析）流動為紊流\n');
else
    fprintf('=>（用微壓計測量所得的資料分析）流動為層流\n');
end
if delta2_pressure_difference_transmitter>=3e6
    fprintf('=>（用壓差變送器測量所得的資料分析）流動為紊流\n');
else
    fprintf('=>（用壓差變送器測量所得的資料分析）流動為層流\n');
end
fprintf('2.用邊界層厚度計算值與理論值比較判定流態：\n');
if abs(delta2_micromanometer-5*sqrt((miu*1e-7*x2)/v_inf))<0.1
    fprintf('=>(用微壓計測量所得的資料分析)流動為層流\n');
elseif abs(delta2_micromanometer-0.37*x2*((miu*1e-7)/(v_inf*x2))^0.2)<0.1
    fprintf('=>(用微壓計測量所得的資料分析)流動為紊流\n');
else
    fprintf('=>測量值與理論值相差較大，無法判定流態！\n');
end
if abs(delta2_pressure_difference_transmitter-5*sqrt((miu*1e-7*x2)/v_inf))<0.1
    fprintf('=>(用微壓計測量所得的資料分析)流動為層流\n');
elseif abs(delta2_pressure_difference_transmitter-0.37*x2*((miu*1e-7)/(v_inf*x2))^0.2)<0.1
    fprintf('=>(用微壓計測量所得的資料分析)流動為紊流\n');
else
    fprintf('=>測量值與理論值相差較大，無法判定流態！\n');
end

(2.2) 執行結果

輸入：
請輸入大氣壓Pa：(單位Pa)98700
請輸入溫度t：(單位℃)26
請輸入流體運動粘性係數：(單位1e-7m^2/s)8.982
請輸入風洞係數K：1.08
請輸入微壓計所選的檔位係數C：0.2
請輸入風洞進口處真空度：(單位mmH2O)22.8
請輸入流場靜壓Pj：(單位毫米酒精柱)118
請輸入第一截面距前緣的距離x：(單位mm)150
請輸入第一截面邊界層內距離的資料(用[]括起來，元素之間用空格分隔)：(單位mm)[0.45 1.45 2.45 3.45 4.45 5.45 6.45 6.95 7.45]
請輸入第一截面用微壓計測得的讀數資料(用[]括起來，元素之間用空格分隔)：(單位毫米酒精柱)[77.5 55 38.8 26 15 10.5 7.5 6.8 6.4]
請輸入第一截面用壓差變送器測得的讀數資料(用[]括起來，元素之間用空格分隔)：(單位Pa)[164.23 118.72 85.559 56.076 34.772 23.870 21.223 20.668 20.202]
請輸入第二截面距前緣的距離x：(單位mm)250
請輸入第二截面邊界層內距離的資料(用[]括起來，元素之間用空格分隔)：(單位mm)[0.45 1.45 2.45 3.45 4.45 5.45 6.45 7.45 8.45 9.45]
請輸入第二截面用微壓計測得的讀數資料(用[]括起來，元素之間用空格分隔)：(單位毫米酒精柱)[83 64 51 39.8 29 21 14 11 9.8 8]
請輸入第二截面用壓差變送器測得的讀數資料(用[]括起來，元素之間用空格分隔)：(單位Pa)[173.85 129.31 104.47 79.882 55.702 40.633 28.892 22.750 19.701 18.167]
輸出：
=>風洞進口流速為：20.490 m/s
>>>>>>>>>>**第一截面資料計算**<<<<<<<<<<
=>雷諾數 Re= 3421782.658
=>微壓計實際讀數計算出來的絕對壓強：(單位Pa)
   1.0e+04 *
  列 1 至 7
    9.8548    9.8592    9.8624    9.8649    9.8671    9.8679    9.8685
  列 8 至 9
    9.8687    9.8687
=>用微壓計測量的資料計算得到的邊界層內的流速：（單位m/s）
  列 1 至 7
   11.7516   14.6568   16.4336   17.7118   18.7408   19.1458   19.4111
  列 8 至 9
   19.4725   19.5075
=>用壓差變送器測量的資料計算得到的邊界層內的流速：（單位m/s）
  列 1 至 7
   10.8004   13.9937   15.9222   17.4588   18.4898   18.9958   19.1166
  列 8 至 9
   19.1418   19.1630
=>用微壓計測量的資料計算得到的邊界層內的速度比：
  列 1 至 7
    0.5735    0.7153    0.8020    0.8644    0.9146    0.9344    0.9474
  列 8 至 9
    0.9504    0.9521
=>用壓差變送器測量的資料計算得到的邊界層內的速度比：
  列 1 至 7
    0.5271    0.6830    0.7771    0.8521    0.9024    0.9271    0.9330
  列 8 至 9
    0.9342    0.9353
=>用微壓計測量得出的邊界層厚度為：5.211 mm
=>用壓差變送器測量得出的邊界層厚度為：5.736 mm
>>>>>>>>>>**第二截面資料計算**<<<<<<<<<<
=>雷諾數 Re= 5702971.096
=>微壓計實際讀數計算出來的絕對壓強：(單位Pa)
   1.0e+04 *
  列 1 至 7
    9.8537    9.8575    9.8600    9.8622    9.8643    9.8659    9.8673
  列 8 至 10
    9.8678    9.8681    9.8684
=>用微壓計測量的資料計算得到的邊界層內的流速：（單位m/s）
  列 1 至 7
   10.9246   13.5696   15.1150   16.3295   17.4207   18.1868   18.8316
  列 8 至 10
   19.1013   19.2081   19.3672
=>用壓差變送器測量的資料計算得到的邊界層內的流速：（單位m/s）
  列 1 至 7
    9.9957   13.3192   14.8532   16.2294   17.4774   18.2120   18.7644
  列 8 至 10
   19.0470   19.1857   19.2551
=>用微壓計測量的資料計算得到的邊界層內的速度比：
  列 1 至 7
    0.5332    0.6623    0.7377    0.7970    0.8502    0.8876    0.9191
  列 8 至 10
    0.9322    0.9375    0.9452
=>用壓差變送器測量的資料計算得到的邊界層內的速度比：
  列 1 至 7
    0.4878    0.6500    0.7249    0.7921    0.8530    0.8888    0.9158
  列 8 至 10
    0.9296    0.9364    0.9398
=>用微壓計測量得出的邊界層厚度為：7.133 mm
=>用壓差變送器測量得出的邊界層厚度為：7.491 mm
>>>>>>>>>>**流態分析**<<<<<<<<<<
一、對第一截面：
1.用雷諾數判定流態：
=>（用微壓計測量所得的資料分析）流動為層流
=>（用壓差變送器測量所得的資料分析）流動為層流
2.用邊界層厚度計算值與理論值比較判定流態：
=>(用微壓計測量所得的資料分析)流動為層流
=>(用微壓計測量所得的資料分析)流動為層流
二、對第二截面：
1.用雷諾數判定流態：
=>（用微壓計測量所得的資料分析）流動為層流
=>（用壓差變送器測量所得的資料分析）流動為層流
2.用邊界層厚度計算值與理論值比較判定流態：
=>(用微壓計測量所得的資料分析)流動為層流
=>(用微壓計測量所得的資料分析)流動為層流

在這裡插入圖片描述

(3) 空氣橫掠單管強制對流換熱實驗

(3.1) 程式碼

clc;clear;
d = 3.73e-3;
l = 99.6e-3;
tw = [40.57 42.41 44.27 46.48 49.72];
tf = [21.12 21.59 21.77 21.77 21.99];
deltaP = [213.5 180.3 144.0 108.5 74.4];
U = [0.3946 0.3979 0.3985 0.3992 0.4004];
I = [14.601 14.689 14.688 14.687 14.685];
fprintf('1. 定性溫度（℃）：');
tm = (tw + tf)/2
fprintf('1.1 以定性溫度查的密度（kg/m^3）：');
rou = (1.165-1.128)/(30-40)*(tm - 30)+1.165
fprintf('1.2 以定性溫度查的導熱係數（W/(m·K)）：');
lamda = (2.76e-2-2.67e-2)/(40-30)*(tm - 30)+2.67e-2
fprintf('1.3 以定性溫度查的運動粘度（m^2/s）：');
miu = ((16.96-16.00)/(40-30)*(tm - 30)+16.00).*1e-6
fprintf('1.4 以定性溫度查的普朗特數：');
Pr = (0.699-0.701)/(40-30)*(tm - 30)+0.701
fprintf('2. 流體速度（m/s）：');
u = sqrt(2*deltaP./rou)
fprintf('3. 雷諾數：');
Re = u.*d./miu
fprintf('4. 實驗測得努塞爾數：');
h = (U.*I)./(pi*d*l)./(tw-tf);
Nu = h.*d./lamda
Nu_cal = zeros(1,5);
for i = 1:5
    if (Re(i)>=40)&&(Re(i)<=4000)
        Nu_cal(i) = 0.683*(Re(i)^0.466)*(Pr(i)^(1/3));
    elseif (Re(i)>4000)&&(Re(i)<=40000)
        Nu_cal(i) = 0.193*(Re(i)^0.618)*(Pr(i)^(1/3));
    end
end
fprintf('5. 努塞爾數理論計算值：');
Nu_cal
fprintf('6. 相對偏差（%%）：');
e = 100*abs(Nu - Nu_cal)./Nu_cal
fprintf('7. 對流換熱量計算（單位J）：');
qConvection = U.*I
fprintf('8. 輻射換熱量估算（取表面發射率為0.5，單位J）：');
qRadiation = 0.5*5.67e-8*(tw.^4-tf.^4)*(pi*d*l)

(3.2) 執行結果

1. 定性溫度（℃）：
tm =
   30.8450   32.0000   33.0200   34.1250   35.8550
1.1 以定性溫度查的密度（kg/m^3）：
rou =
    1.1619    1.1576    1.1538    1.1497    1.1433
1.2 以定性溫度查的導熱係數（W/(m·K)）：
lamda =
    0.0268    0.0269    0.0270    0.0271    0.0272
1.3 以定性溫度查的運動粘度（m^2/s）：
miu =
   1.0e-04 *
    0.1608    0.1619    0.1629    0.1640    0.1656
1.4 以定性溫度查的普朗特數：
Pr =
    0.7008    0.7006    0.7004    0.7002    0.6998
2. 流體速度（m/s）：
u =
   19.1705   17.6495   15.7989   13.7382   11.4081
3. 雷諾數：
Re =
   1.0e+03 *
    4.4466    4.0658    3.6176    3.1254    2.5693
4. 實驗測得努塞爾數：
Nu =
   35.3561   33.3770   30.8240   28.0113   24.8891
5. 努塞爾數理論計算值：
Nu_cal =
   30.8018   29.1405   27.6115   25.7897   23.5354
6. 相對偏差（%）：
e =
   14.7859   14.5381   11.6346    8.6141    5.7518
7. 對流換熱量計算（單位J）：
qConvection =
    5.7616    5.8448    5.8532    5.8631    5.8799
8. 輻射換熱量估算（取表面發射率為0.5，單位J）：
qRadiation =
   1.0e-03 *
0.0831    0.0999    0.1197    0.1470    0.1945

(4) 水平管外自然對流換熱實驗

(4.1) 程式碼

%用最小二乘擬合計算出實驗關聯式
clc;clear;
syms LnGrPr;
GrPr = [458934.6 128178.6 81973.53 47893.61 1903764 1523845 241316.5 586097.9 811834.9 1291265];
Nu = [12.77 9.53 9.42 8.02 17.51 17.01 11.39 15.16 13.1013 15.54];
ln_GrPr = log(GrPr);
ln_Nu = log(Nu);
p = polyfit(ln_GrPr,ln_Nu,1); %1次多項式擬合（即線性擬合）
LnNu = vpa(poly2sym(p,LnGrPr),5)
n = vpa(p(1),5)
C = vpa(exp(p(2)),5) %求n和C

%雙對數座標圖繪製曲線：
clc;clear;
GrPr = [458934.6 128178.6 81973.53 47893.61 1903764 1523845 241316.5 586097.9 811834.9 1291265];
Nu = [12.77 9.53 9.42 8.02 17.51 17.01 11.39 15.16 13.1013 15.54];
figure;
loglog(GrPr,Nu,'*');
hold on;
grid on;
p = polyfit(log(GrPr),log(Nu),1);
x = min(GrPr):0.1:max(GrPr);
y = exp(polyval(p,log(x)));
loglog(x,y,'-');
%legend('原始資料點','擬合曲線','location','northwest');
xlabel('GrPr');
ylabel('Nu');
title('水平管外自然對流換熱實驗曲線');
scatter(GrPr(9),Nu(9),'o');
legend('原始資料點','擬合曲線','3#組資料點','location','northwest');
gtext(['3#','(',num2str(GrPr(9)),',',num2str(Nu(9)),')']);

(4.2) 執行結果

LnNu =
0.20561*LnGrPr - 0.12204 
n = 
0.20561
C =
0.88511

在這裡插入圖片描述

原始碼作者：
  Aiden Lee
特別宣告：
  文章僅供學習參考，轉載請註明出處，嚴禁盜用！

相關文章