天天看點

2018年高教社杯A題 高溫作業專用服裝設計

  • 題目

    在高溫環境下工作時,人們需要穿着專用服裝以避免灼傷。專用服裝通常由三層織物材料構成,記為I、II、III層,其中I層與外界環境接觸,III層與皮膚之間還存在空隙,将此空隙記為IV層。

    為設計專用服裝,将體内溫度控制在37ºC的假人放置在實驗室的高溫環境中,測量假人皮膚外側的溫度。為了降低研發成本、縮短研發周期,請你們利用數學模型來确定假人皮膚外側的溫度變化情況,并解決以下問題:

(1) 專用服裝材料的某些參數值由附件1給出,對環境溫度為75ºC、II層厚度為6 mm、IV層厚度為5 mm、工作時間為90分鐘的情形開展實驗,測量得到假人皮膚外側的溫度(見附件2)。建立數學模型,計算溫度分布,并生成溫度分布的Excel檔案(檔案名為problem1.xlsx)。

(2) 當環境溫度為65ºC、IV層的厚度為5.5 mm時,确定II層的最優厚度,確定工作60分鐘時,假人皮膚外側溫度不超過47ºC,且超過44ºC的時間不超過5分鐘。

(3) 當環境溫度為80時,确定II層和IV層的最優厚度,確定工作30分鐘時,假人皮膚外側溫度不超過47ºC,且超過44ºC的時間不超過5分鐘。

附件1. 專用服裝材料的參數值

附件2. 假人皮膚外側的測量溫度

題目下載下傳位址:http://special.univs.cn/service/jianmo/sxjmtmhb/2018/0709/1179372.shtml

  • 解答如下:
  • 假設:

    假設 1: 假設在整個過程中外部溫度保持不變;

    假設 2: 假設在整個的過程中隻考慮熱傳導這一種熱傳遞方式, 熱對流以及熱輻射可以

    忽略;

    假設 3: 假設織物材料無褶皺且其表面溫度處處相同;

    假設 4: 假設在實驗中熱傳導的方向與假人人體垂直;

    假設 5: 假設假人可以承受實驗中所達到的所有可能溫度;

    假設 6: 假設每一層的參數固定不變, 即不會随着溫度的改化發生變化;

    假設 7: 假設專用服裝的初始溫度與假人體内溫度一樣, 為 37℃;

    問題(1)

    為了求出溫度的分布, 我們需要通過建立溫度——外部到假人皮膚之間的距離——時間的數學模型(示意圖如下圖)

    2018年高教社杯A題 高溫作業專用服裝設計

    用微元法來建立熱防護服系統的熱量傳遞模型

    微元[x,x+Δx]在時間段[t,t+Δt]内吸收的熱量為:

    2018年高教社杯A題 高溫作業專用服裝設計
    Fourier熱傳導定律
    2018年高教社杯A題 高溫作業專用服裝設計
    Q1=Q2
    2018年高教社杯A題 高溫作業專用服裝設計

    上述方程在各層之間的區間内成立

    邊界條件:

    左邊界

    2018年高教社杯A題 高溫作業專用服裝設計
    右邊界
    2018年高教社杯A題 高溫作業專用服裝設計
    交界面
    2018年高教社杯A題 高溫作業專用服裝設計
    熱防護服系統熱量傳遞的偏微分方程模型:
    2018年高教社杯A題 高溫作業專用服裝設計

    再通過有限差分格式求解:

    ks為第四層與皮膚之間的熱交換系數 ks=8.36

    ke為第I層與高溫環境之間的熱交換系數 ke=117.40

    2018年高教社杯A題 高溫作業專用服裝設計
    2018年高教社杯A題 高溫作業專用服裝設計

    (2)這裡用的二分法,在滿足條件的情況下,求最小的厚度

    L2=17.6609mm

    (3)這問要求第二層,第四層的最優厚度,雙目标最優函數,看了一下各路大神的思路,無非就是遺傳算法,暴力窮舉,個人認為這裡将主要考慮的是産品的成本,皮膚外側溫度随着第二層、第四層厚度的增加而降低的,這個是可以驗證的,第四層越厚,第二層必然就會薄一些, 成本自然就低一些,而第四層厚度的增加,所需要的服裝材料幾乎沒有變化。是以我這裡将第四層的厚度取最大,将第四層厚度确定之後,這個問題就轉化成第二問的問題了,這樣就很簡單了。

    L2=19.2813mm

    L4=6.4mm

  • 具體代碼如下:

問題(1)

%問題一的主函數 T1.m

clc
clear;
T=1;      %時間步長
m=[6 60 36 50];%将每一層細分為0.1mm每塊
%m=[30 300 180 250] %将每一層細分為0.02mm每塊
ms=[0 0 0 0];          % ms表示總的分層數
ms(1,1)=m(1,1);
for j=2:4
    ms(1,j)=ms(1,j-1)+m(1,j);
end
l=[0.6 6.6 10.2 15.2]*10^(-3);        %l1 l2 l3 l4z的坐标
p=[300 862 74.2 1.18];       %各層密度
c=[1377 2100 1726 1005];     %比熱
k=[0.082 0.37 0.045 0.028];  %熱傳導率
u0=37.0;          %初始溫度
ue=75;
%求空間步長
x=zeros(1,4);   
x(1,1)=l(1,1)/m(1,1);
for j=2:4
    x(1,j)=(l(1,j)-l(1,j-1))/m(1,j);
end
 % M(1,j) 從第1層到第j層的等分數
M(1,1)=m(1,1);
for j=2:4
    M(1,j)=M(1,j-1)+m(1,j);     
end
%求式(10)裡的λ,用b表示λ
for j=1:4
    b(1,j)=(k(1,j)*T)/(c(1,j)*p(1,j)*(x(1,j)^2));
end
ke=117.4035;    
e=(ke*x(1,1)/k(1,1));   %求式(12)的ue 用e表示
ks=8.36;
us=(ks*x(1,4)/k(1,4));   %求式(13)的us
for j=1:4
    v(1,j)=k(1,j)/x(1,j);   %求式(14)的v
end
u=zeros(5400/T+1,ms(1,4)+1);   %建立一個以距離為x軸,時間為y軸的坐标軸
u(1,:)=37;           %時間t=0時,各點溫度為37度
f=zeros(ms(1,4)+1,1);        %AX=f 
x01=e+1;               %l0的邊界條件
x02=-1;  
x1=[-b(1,1) 1+2*b(1,1) -b(1,1)];    %第一層  
x2=[-b(1,2) 1+2*b(1,2) -b(1,2)];    %第二層       
x3=[-b(1,3) 1+2*b(1,3) -b(1,3)];   %第三層   
x4=[-b(1,4) 1+2*b(1,4) -b(1,4)];   %第四層 
x12=[-v(1,1) v(1,1)+v(1,2) -v(1,2)];  %12層邊界
x23=[-v(1,2) v(1,2)+v(1,3) -v(1,3)];%23層邊界
x34=[-v(1,3) v(1,3)+v(1,4) -v(1,4)];%34層邊界
x51=-1;              %l4的邊界條件
x52=1+us;
A=zeros(ms(1,4)+1,ms(1,4)+1);
A(1,1)=x01;
A(1,2)=x02;
f(1,1)=e*ue;
for i=2:ms(1,1)
    A(i,i-1:i+1)=x1;
    f(i,1)=37;
end
A(ms(1,1)+1,ms(1,1):ms(1,1)+2)=x12;
f(ms(1,1)+1,1)=0;
for i=ms(1,1)+2:ms(1,2)
     A(i,i-1:i+1)=x2;
     f(i,1)=37;
end
A(ms(1,2)+1,ms(1,2):ms(1,2)+2)=x23;
f(ms(1,2)+1,1)=0;
for i=ms(1,2)+2:ms(1,3)
     A(i,i-1:i+1)=x3;
     f(i,1)=37;
end
A(ms(1,3)+1,ms(1,3):ms(1,3)+2)=x34;
f(ms(1,3)+1,1)=0;
for i=ms(1,3)+2:ms(1,4)
     A(i,i-1:i+1)=x4;
     f(i,1)=37;
end
A(ms(1,4)+1,ms(1,4))=x51;
A(ms(1,4)+1,ms(1,4)+1)=x52;
f(ms(1,4)+1,1)=us*u0;
 
%追趕法
a=zeros(1,ms(1,4)+1);%下對角元素
a(1,1)=0;
for i=2:ms(1,4)+1
    a(1,i)=A(i,i-1);
end
c=zeros(1,ms(1,4)+1);%上對角元素
c(1,ms(1,4)+1)=0;
for i=1:ms(1,4)
   c(1,i)=A(i,i+1);   
end
b=zeros(1,ms(1,4)+1);
for i=1:ms(1,4)+1
    b(1,i)=A(i,i);
end
y=chase(a,b,c,f);
ans=zeros(5400/T+1,ms(1,4)+1); %結果矩陣
ans(1,:)=u0;
ans(2,:)=y;
%循環求結果
for i=2:5400/T
f=y';
f(1,1)=e*ue;
f(ms(1,1)+1,1)=0;
f(ms(1,2)+1,1)=0;
f(ms(1,3)+1,1)=0;
f(ms(1,4)+1,1)=us*u0;
y=chase(a,b,c,f);
ans(i+1,:)=y;
end
xlswrite('problem1.xlsx',ans,'Sheet2','A1');
figure(1)
%溫度随時間、空間的分布圖
t=0:T:5400;
y=0:ms(1,4);
[b1,a1]=meshgrid(y,t);
mesh(b1,a1,ans);
title('溫度随時間、空間的分布圖');
%xlabel('空間/0.02mm');
xlabel('空間/0.1mm');
ylabel('時間/s');
zlabel('溫度/℃');
figure(2)
%分界面的溫度随時間的變化圖
t=0:T:5400;
%左邊界
l1=ans(:,1);
%第一層和第二層邊界
l2=ans(:,ms(1,1)+1);
%第二層和第三層邊界
l3=ans(:,ms(1,2)+1);
%第三層和第四層邊界
l4=ans(:,ms(1,3)+1);
%右邊界I、II、III、IV
l5=ans(:,ms(1,4)+1);
plot(t,l1,'m-',t,l2,'g-',t,l3,'r-',t,l4,'c-',t,l5,'k-');
xlabel('時間/s');
ylabel('溫度/℃');
title('分界面的溫度随時間的變化');
legend('左邊界','I、II層邊界','II、III層邊界','III、IV層邊界','右邊界');
           

%追趕法 chase.m

function x=chase(a,b,c,f) 
%求解線性方程組 Ax=f, 其中 A是三對角陣
%a是矩陣 A的下對角線元素 a(1)=0
%b是矩陣 A的對角線元素
%c是矩陣 A的上對角線元素 c(n)=0
%f是方程組的右端向量
n=length(f); 
x=zeros(1,n);
y=zeros(1,n); 
d=zeros(1,n);
z= zeros(1,n); 
%預處理
d(1)=b(1); 
for i=1:n-1 
z(i)=c(i)/d(i); 
d(i+1)=b(i+1)-a(i+1)*z(i); 
end
%追的過程
y(1)=f(1)/d(1); 
for i=2:n 
 y(i)=(f(i)-a(i)*y(i-1))/d(i); 
end
%趕的過程
x(n)=y(n); 
for i=n-1:-1:1 
x(i)=y(i)-z(i)*x(i+1); 
end
           

% 求ke,ks quedingxishu.m

clc,clear;
a=load('TT.txt'); %TT.txt 存儲的是附件2的所有數字資訊
ue=75; %環境溫度
us=37; %皮膚内測溫度
ks=8.36; %與皮膚的熱交換系數
keq=1;%最合适的 與環境的熱交換系數
ksq=1;%最合适的 與皮膚的熱交換系數
min=100000;  %初始值
k=[0.082 0.37 0.045 0.028];  %熱傳導率
l=[0.6 6.6 10.2 15.2]*10^-3;
sum=l(1,1)/k(1,1);
for j=2:4  
    sum=sum+(l(1,j)-l(1,j-1))/k(1,j);
end
um=48.08;    %穩态溫度
time=a(:,2);
time=time';
for ks=7.8:0.02:8.8
    ke=1/((ue-um)/(ks*(um-us))-sum); %與外界的熱交換系數
    w=ul4(ke,ks);
    min1=0;
    for i=1:5041
        min1=min1+1/2*(time(1,i)-w(1,i))^2;
    end
    if min1<min
        min=min1;
        keq=ke;
        ksq=ks;
    end
end
min
keq %最合适的 與環境的熱交換系數
ksq %最合适的 與皮膚的熱交換系數
           

%問題1 傳回穩态的溫度用于求ke,ks的值 ul4.m

%問題1 傳回穩态的溫度用于求ke,ks ul4.m
function [q] = ul4( ke,ks )
T=1;      %時間步長
m=[6 60 36 50];%将每一層細分為0.1mm每塊
%m=[30 300 180 250] %将每一層細分為0.02mm每塊
ms=[0 0 0 0];          % ms表示總的分層數
ms(1,1)=m(1,1);
for j=2:4
    ms(1,j)=ms(1,j-1)+m(1,j);
end
l=[0.6 6.6 10.2 15.2]*10^(-3);  %l1 l2 l3 l4z的坐标
p=[300 862 74.2 1.18];       %各層密度
c=[1377 2100 1726 1005];     %比熱
k=[0.082 0.37 0.045 0.028];  %熱傳導率
u0=37.0;          %初始溫度
ue=75;
%求空間步長
x=zeros(1,4);   
x(1,1)=l(1,1)/m(1,1);
for j=2:4
    x(1,j)=(l(1,j)-l(1,j-1))/m(1,j);
end
 % M(1,j) 從第1層到第j層的等分數
M(1,1)=m(1,1);
for j=2:4
    M(1,j)=M(1,j-1)+m(1,j);     
end
%求式(10)裡的λ,用b表示λ
for j=1:4
    b(1,j)=(k(1,j)*T)/(c(1,j)*p(1,j)*(x(1,j)^2));
end
e=(ke*x(1,1)/k(1,1));   %求式(12)的ue 用e表示
us=(ks*x(1,4)/k(1,4));   %求式(13)的us
for j=1:4
    v(1,j)=k(1,j)/x(1,j);   %求式(14)的v
end
u=zeros(5400/T+1,ms(1,4)+1);   %建立一個以距離為x軸,時間為y軸的坐标軸
u(1,:)=37;           %時間t=0時,各點溫度為37度

f=zeros(ms(1,4)+1,1);        %AX=f 
x01=e+1;               %l0的邊界條件
x02=-1;  

x1=[-b(1,1) 1+2*b(1,1) -b(1,1)];    %第一層  
x2=[-b(1,2) 1+2*b(1,2) -b(1,2)];    %第二層       
x3=[-b(1,3) 1+2*b(1,3) -b(1,3)];   %第三層   
x4=[-b(1,4) 1+2*b(1,4) -b(1,4)];   %第四層 
x12=[-v(1,1) v(1,1)+v(1,2) -v(1,2)];  %12層邊界
x23=[-v(1,2) v(1,2)+v(1,3) -v(1,3)];%23層邊界
x34=[-v(1,3) v(1,3)+v(1,4) -v(1,4)];%34層邊界

x51=-1;              %l4的邊界條件
x52=1+us;
A=zeros(ms(1,4)+1,ms(1,4)+1);
A(1,1)=x01;
A(1,2)=x02;
f(1,1)=e*ue;
for i=2:ms(1,1)
    A(i,i-1:i+1)=x1;
    f(i,1)=37;
end
A(ms(1,1)+1,ms(1,1):ms(1,1)+2)=x12;
f(ms(1,1)+1,1)=0;
for i=ms(1,1)+2:ms(1,2)
     A(i,i-1:i+1)=x2;
     f(i,1)=37;
end
A(ms(1,2)+1,ms(1,2):ms(1,2)+2)=x23;
f(ms(1,2)+1,1)=0;
for i=ms(1,2)+2:ms(1,3)
     A(i,i-1:i+1)=x3;
     f(i,1)=37;
end
A(ms(1,3)+1,ms(1,3):ms(1,3)+2)=x34;
f(ms(1,3)+1,1)=0;
for i=ms(1,3)+2:ms(1,4)
     A(i,i-1:i+1)=x4;
     f(i,1)=37;
end
A(ms(1,4)+1,ms(1,4))=x51;
A(ms(1,4)+1,ms(1,4)+1)=x52;
f(ms(1,4)+1,1)=us*u0;

%追趕法
a=zeros(1,ms(1,4)+1);%下對角元素
a(1,1)=0;
for i=2:ms(1,4)+1
    a(1,i)=A(i,i-1);
end
c=zeros(1,ms(1,4)+1);%上對角元素
c(1,ms(1,4)+1)=0;
for i=1:ms(1,4)
   c(1,i)=A(i,i+1);   
end
b=zeros(1,ms(1,4)+1);
for i=1:ms(1,4)+1
    b(1,i)=A(i,i);
end
y=chase(a,b,c,f);
ans=zeros(5400/T+1,ms(1,4)+1); %結果矩陣
ans(1,:)=u0;
ans(2,:)=y;
%循環求結果
for i=2:5400/T
f=y';
f(1,1)=e*ue;
f(ms(1,1)+1,1)=0;
f(ms(1,2)+1,1)=0;
f(ms(1,3)+1,1)=0;
f(ms(1,4)+1,1)=us*u0;
y=chase(a,b,c,f);
ans(i+1,:)=y;
end
q=ans(:,ms(1,4)+1);
q=q';
end
           

問題(2)

%問題2的主函數 T2.m

clc,clear;
len=0.1; %二分法結束的條件
min=0.6;
max=25;
ans=0;%結果
q=fun2(min);%求解最小厚度情況下是否滿足限制條件
if q(1,3601)<=47&&q(1,3301)<=44
   ans=min;
   fprintf('最小厚度情況下滿足限制條件');
end
q=fun2(max);%求解最大厚度情況下是否滿足限制條件
if q(1,3601)>47||q(1,3301)>44
   fprintf('最大厚度情況下不滿足限制條件,此題無解'); 
end 
while(max-min>len)
  mid=(min+max)/2;
  q=fun2(mid);
  if q(1,3601)<=47&&q(1,3301)<=44
  max=mid;
  else
  min=mid;
  end
end
ans=max
           

%對于不同的L2,傳回皮膚外側溫度變化 fun2.m

function [q]=fun2(l2)
%l2 傳入II層的厚度
T=1;%時間步長
s=ceil(l2/0.1);
m=zeros(1,4);
%将每一層細分為0.1mm每塊
m(1,1)=6;
m(1,2)=s;
m(1,3)=36;
m(1,4)=55; 
ms=[0 0 0 0];          % ms表示總的分層數
ms(1,1)=m(1,1);
for j=2:4
    ms(1,j)=ms(1,j-1)+m(1,j);
end
l=[0.6 l2+0.6 l2+4.2 l2+9.7]*10^(-3);        %l1 l2 l3 l4的坐标
p=[300 862 74.2 1.18];       %各層密度
c=[1377 2100 1726 1005];     %比熱
k=[0.082 0.37 0.045 0.028];  %熱傳導率
u0=37.0;          %初始溫度
ue=65;            %環境溫度
%求空間步長
x=zeros(1,4);   
x(1,1)=l(1,1)/m(1,1);
for j=2:4
    x(1,j)=(l(1,j)-l(1,j-1))/m(1,j);
end
 % M(1,j) 從第1層到第j層的等分數
M(1,1)=m(1,1);
for j=2:4
    M(1,j)=M(1,j-1)+m(1,j);     
end
%求式(10)裡的λ,用b表示λ
for j=1:4
    b(1,j)=(k(1,j)*T)/(c(1,j)*p(1,j)*(x(1,j)^2));
end
ke=117.4035;    
e=(ke*x(1,1)/k(1,1));   %求式(12)的ue 用e表示
ks=8.36;
us=(ks*x(1,4)/k(1,4));   %求式(13)的us
for j=1:4
    v(1,j)=k(1,j)/x(1,j);   %求式(14)的v
end
u=zeros(5400/T+1,ms(1,4)+1);   %建立一個以距離為x軸,時間為y軸的坐标軸
u(1,:)=37;           %時間t=0時,各點溫度為37度
 
f=zeros(ms(1,4)+1,1);        %AX=f 
x01=e+1;               %l0的邊界條件
x02=-1;  
x1=[-b(1,1) 1+2*b(1,1) -b(1,1)];    %第一層  
x2=[-b(1,2) 1+2*b(1,2) -b(1,2)];    %第二層       
x3=[-b(1,3) 1+2*b(1,3) -b(1,3)];   %第三層   
x4=[-b(1,4) 1+2*b(1,4) -b(1,4)];   %第四層 
x12=[-v(1,1) v(1,1)+v(1,2) -v(1,2)];  %12層邊界
x23=[-v(1,2) v(1,2)+v(1,3) -v(1,3)];%23層邊界
x34=[-v(1,3) v(1,3)+v(1,4) -v(1,4)];%34層邊界
x51=-1;              %l4的邊界條件
x52=1+us;
A=zeros(ms(1,4)+1,ms(1,4)+1);
A(1,1)=x01;
A(1,2)=x02;
f(1,1)=e*ue;
for i=2:ms(1,1)
    A(i,i-1:i+1)=x1;
    f(i,1)=37;
end
A(ms(1,1)+1,ms(1,1):ms(1,1)+2)=x12;
f(ms(1,1)+1,1)=0;
for i=ms(1,1)+2:ms(1,2)
     A(i,i-1:i+1)=x2;
     f(i,1)=37;
end
A(ms(1,2)+1,ms(1,2):ms(1,2)+2)=x23;
f(ms(1,2)+1,1)=0;
for i=ms(1,2)+2:ms(1,3)
     A(i,i-1:i+1)=x3;
     f(i,1)=37;
end
A(ms(1,3)+1,ms(1,3):ms(1,3)+2)=x34;
f(ms(1,3)+1,1)=0;
for i=ms(1,3)+2:ms(1,4)
     A(i,i-1:i+1)=x4;
     f(i,1)=37;
end
A(ms(1,4)+1,ms(1,4))=x51;
A(ms(1,4)+1,ms(1,4)+1)=x52;
f(ms(1,4)+1,1)=us*u0;
%追趕法
a=zeros(1,ms(1,4)+1);%下對角元素
a(1,1)=0;
for i=2:ms(1,4)+1
    a(1,i)=A(i,i-1);
end
c=zeros(1,ms(1,4)+1);%上對角元素
c(1,ms(1,4)+1)=0;
for i=1:ms(1,4)
   c(1,i)=A(i,i+1);   
end
b=zeros(1,ms(1,4)+1);
for i=1:ms(1,4)+1
    b(1,i)=A(i,i);
end
y=chase(a,b,c,f);
ans=zeros(5400/T+1,ms(1,4)+1); %結果矩陣
ans(1,:)=u0;
ans(2,:)=y;
%循環求結果
for i=2:5400/T
f=y';
f(1,1)=e*ue;
f(ms(1,1)+1,1)=0;
f(ms(1,2)+1,1)=0;
f(ms(1,3)+1,1)=0;
f(ms(1,4)+1,1)=us*u0;
y=chase(a,b,c,f);
ans(i+1,:)=y;
end
q=ans(:,ms(1,4)+1);
q=q';
end
           

%求假人皮膚外側穩定溫度随II層厚度的變化 draw.m

clc,clear;
i=0;
for l2=0.6:0.1:25
    q=fun2(l2);
    i=i+1;
    a(1,i)=q(1,length(q));   
end
l2=0.6:0.1:25;
plot(l2,a,'c');
title('假人皮膚外側穩定溫度随II層厚度的變化');
xlabel('II層厚度/mm');
ylabel('假人皮膚外側穩定溫度/攝氏度');
           

問題(3)

%問題3的主函數 T3.m

clc,clear;
len=0.1; %二分法結束的條件
min=0.6;
max=25;
ans=0;%結果
q=fun3(min);%求解最小厚度情況下是否滿足限制條件
if q(1,1801)<=47&&q(1,1501)<=44
   ans=min;
   fprintf('最小厚度情況下滿足限制條件');
end
q=fun3(max);%求解最大厚度情況下是否滿足限制條件
if q(1,1801)>47||q(1,1501)>44
   fprintf('最大厚度情況下不滿足限制條件,此題無解'); 
end
while(max-min>len)
  mid=(min+max)/2;
  q=fun3(mid);
  if q(1,1801)<=47&&q(1,1501)<=44
  max=mid;
  else
  min=mid;
  end
end
ans=max
           

%傳入不同的L2,傳回皮膚外側的溫度變化矩陣 fun3.m

function [q]=fun3(l2)
%l2 傳入II層的厚度
T=1;%時間步長
s=ceil(l2/0.1);
m=zeros(1,4);
%将每一層細分為0.1mm每塊
m(1,1)=6;
m(1,2)=s;
m(1,3)=36;
m(1,4)=64; %第四層取6.4mm
ms=[0 0 0 0];          % ms表示總的分層數
ms(1,1)=m(1,1);
for j=2:4
    ms(1,j)=ms(1,j-1)+m(1,j);
end
l=[0.6 l2+0.6 l2+4.2 l2+10.6]*10^(-3);        %l1 l2 l3 l4的坐标
p=[300 862 74.2 1.18];       %各層密度
c=[1377 2100 1726 1005];     %比熱
k=[0.082 0.37 0.045 0.028];  %熱傳導率
u0=37.0;          %初始溫度
ue=80;            %環境溫度
%求空間步長
x=zeros(1,4);
x(1,1)=l(1,1)/m(1,1);
for j=2:4
    x(1,j)=(l(1,j)-l(1,j-1))/m(1,j);
end
 % M(1,j) 從第1層到第j層的等分數
M(1,1)=m(1,1);
for j=2:4
    M(1,j)=M(1,j-1)+m(1,j);     
end
%求式(10)裡的λ,用b表示λ
for j=1:4
    b(1,j)=(k(1,j)*T)/(c(1,j)*p(1,j)*(x(1,j)^2));
end
ke=117.4035;    
e=(ke*x(1,1)/k(1,1));   %求式(12)的ue 用e表示
ks=8.36;
us=(ks*x(1,4)/k(1,4));   %求式(13)的us
for j=1:4
    v(1,j)=k(1,j)/x(1,j);   %求式(14)的v
end
u=zeros(5400/T+1,ms(1,4)+1);   %建立一個以距離為x軸,時間為y軸的坐标軸
u(1,:)=37;           %時間t=0時,各點溫度為37度
 
f=zeros(ms(1,4)+1,1);        %AX=f 
x01=e+1;               %l0的邊界條件
x02=-1;  
x1=[-b(1,1) 1+2*b(1,1) -b(1,1)];    %第一層  
x2=[-b(1,2) 1+2*b(1,2) -b(1,2)];    %第二層       
x3=[-b(1,3) 1+2*b(1,3) -b(1,3)];   %第三層   
x4=[-b(1,4) 1+2*b(1,4) -b(1,4)];   %第四層 
x12=[-v(1,1) v(1,1)+v(1,2) -v(1,2)];  %12層邊界
x23=[-v(1,2) v(1,2)+v(1,3) -v(1,3)];%23層邊界
x34=[-v(1,3) v(1,3)+v(1,4) -v(1,4)];%34層邊界
x51=-1;              %l4的邊界條件
x52=1+us;
A=zeros(ms(1,4)+1,ms(1,4)+1);
A(1,1)=x01;
A(1,2)=x02;
f(1,1)=e*ue;
for i=2:ms(1,1)
    A(i,i-1:i+1)=x1;
    f(i,1)=37;
end
A(ms(1,1)+1,ms(1,1):ms(1,1)+2)=x12;
f(ms(1,1)+1,1)=0;
for i=ms(1,1)+2:ms(1,2)
     A(i,i-1:i+1)=x2;
     f(i,1)=37;
end
A(ms(1,2)+1,ms(1,2):ms(1,2)+2)=x23;
f(ms(1,2)+1,1)=0;
for i=ms(1,2)+2:ms(1,3)
     A(i,i-1:i+1)=x3;
     f(i,1)=37;
end
A(ms(1,3)+1,ms(1,3):ms(1,3)+2)=x34;
f(ms(1,3)+1,1)=0;
for i=ms(1,3)+2:ms(1,4)
     A(i,i-1:i+1)=x4;
     f(i,1)=37;
end
A(ms(1,4)+1,ms(1,4))=x51;
A(ms(1,4)+1,ms(1,4)+1)=x52;
f(ms(1,4)+1,1)=us*u0;
%追趕法
a=zeros(1,ms(1,4)+1);%下對角元素
a(1,1)=0;
for i=2:ms(1,4)+1
    a(1,i)=A(i,i-1);
end
c=zeros(1,ms(1,4)+1);%上對角元素
c(1,ms(1,4)+1)=0;
for i=1:ms(1,4)
   c(1,i)=A(i,i+1);   
end
b=zeros(1,ms(1,4)+1);
for i=1:ms(1,4)+1
    b(1,i)=A(i,i);
end
y=chase(a,b,c,f);
ans=zeros(5400/T+1,ms(1,4)+1); %結果矩陣
ans(1,:)=u0;
ans(2,:)=y;
%循環求結果
for i=2:5400/T
f=y';
f(1,1)=e*ue;
f(ms(1,1)+1,1)=0;
f(ms(1,2)+1,1)=0;
f(ms(1,3)+1,1)=0;
f(ms(1,4)+1,1)=us*u0;
y=chase(a,b,c,f);
ans(i+1,:)=y;
end
q=ans(:,ms(1,4)+1);
q=q';
end
           
代碼都是自己花時間碼出來的,問題應該是沒有的,不足就是重複性比較高。非常感謝大家的關注與提問,因本題是2018國賽的題目,網上的優秀論文也比較多,花個幾天時間去了解,算法也比較容易實作。因為作者在準備考研,時間比較緊張,不能及時回複大家,還請多多見諒。
作者有幸參加了2019年高教社杯數學模組化,運氣比較好,獲得了湖南省一等獎、全國二等獎的成績,就自己最大的感受而言,抛開團隊及其他因素,前期的積澱是比較重要的,數學模組化算法及應用這本書大部分内容要去弄懂,代碼盡量去實作,如果學校有教育訓練,盡量去參加,不要劃水,我所了解到的劃水的最後至多得了個省二省三,該付出時間的你一個都跑不掉。今年國賽我也報名了,希望我和各位都能取得一個滿意的成績。

繼續閱讀