2023年华数杯数学建模A题隔热材料的结构优化控制研究解题全过程文档及程序

2023年华数杯全国大学生数学建模

A题 隔热材料的结构优化控制研究

原题再现：

  新型隔热材料 A 具有优良的隔热特性，在航天、军工、石化、建筑、交通等高科技领域中有着广泛的应用。
  目前，由单根隔热材料 A 纤维编织成的织物，其热导率可以直接测出；但是单根隔热材料 A 纤维的热导率（本题实验环境下可假定其为定值），因其直径过小，长径比（长度与直径的比值）较大，无法直接测量。单根纤维导热性能是织物导热性能的基础，也是建立基于纤维的各种织物导热模型的基础。建立一个单根隔热材料 A 纤维的热导率与织物整体热导率的传热机理模型成为研究重点。该模型不仅能得到单根隔热材料 A 纤维的热导率，解决当前单根 A 纤维热导率无法测量的技术难题；而且在建立的单根隔热材料 A 纤维热导率与织物热导率的关系模型的基础上，调控织物的编织结构，进行优化设计，能制作出更好的满足在航天、军工、石化、建筑、交通等高科技领域需求的优异隔热性能织物。
  织物是由大量单根纤维堆叠交织在一起形成的网状结构，本题只研究平纹织物，如图 1 和图 2 所示。不同直径纤维制成的织物，其基础结构参数不同，即纤维弯曲角度、织物厚度、经密、纬密等不同，从而影响织物的导热性能。本题，假设任意单根 A 纤维的垂直切面为圆形，织物中每根纤维始终为一个有弯曲的圆柱。经纱、纬纱弯曲角度 10° < θ ≤ 26.565°。
  热导率是纤维和织物物理性质中最重要的指标之一。织物的纤维之间存在空隙，空隙里空气为静态空气，静态空气热导率 0.0296 W/(mK)。计算织物热导率时既考虑纤维之间的传热，也不能忽略空隙中空气的传热。


  我们在 25℃实验室环境下，用 Hotdisk 装置对织物进行加热和测量，Hotdisk恒定功率为 1mW，作用时间 1s，在 0.1s 时热流恰好传递到织物另一侧。实验测得 0~0.1s 之间织物位于热源一侧的温度随时间变化的数据见附件 1。
  请建立数学模型，回答下列问题:
  问题 1：假设附件 1 的温度为热源侧织物的表面温度，只考虑纤维传热和空隙间的气体传热，建立平纹织物整体热导率与单根纤维热导率之间关系的数学模型。在附件 2 的实验样品参数条件下，测得如图 2 所示的平纹织物的整体热导率为 0.033W/(mK)，请根据建立的数学模型计算出单根 A 纤维的热导率。
  问题 2：假设：1）制成织物的任意单根 A 纤维的直径在 0.3mm~0.6mm。2）织物位于热源一侧表面温度随时间的变化的数据依旧参考附件 1。3）由于温度和织物结构造成的织物整体密度和比热的变化可以忽略。请问如何选用单根 A纤维的直径及调整织物的经密、纬密、弯曲角度，使得织物的整体热导率最低。
  问题 3：如果附件 1 的温度实际是热源侧织物表面空气的温度，此时该侧就会发生对流换热，假定织物表面的对流换热系数为 50 W/(m2K)，请重新解答问题一和问题二。

整体求解过程概述(摘要)

  本文研究了平纹织物整体的热导率和单根 A 纤维的热导率之间模型，对题目所给的表面温度分布进行插值，对于纤维材料建立了整体分析和微观最小单位分析两种模型，通过建立热传导方程和边界条件并且利用两种模型的参数进行温度分布的计算，使得两种模型下的温度分布最接近则建立了整体热导率和单根A纤维热导率之间的模型。
  对于问题一，根据热力学定律在厚度方向上建立一维热传导方程，利用有限差分法对时间和厚度进行差分，建立向后隐式差分格式，利用数值积分的方法求出温度分布，求得整体分析的温度分布图如图3所示。根据平纹织物的几何形状建立微观最小单位几何模型，根据几何形状求得单根 A 纤维和该单位中空气的体积分布，根据比热和质量比例的关系求得单根 A 纤维的密度和比热的乘积。建立初始条件和边界条件，用有限差分法对方程和条件进行差分，数值积分求出温度分布。计算两个模型下各点温度的残差平方并建立单目标最优化模型，利用黄金分割法对单根 A 纤维的热导率进行搜索，求得最优情况下的温度分布如图4 所示。求得最优的残差平方和为2.2796，最后可以求得单根A纤维的热导率为0.0564。
  对于问题二，根据微观最小单位模型建立整体织物密度比热的乘积和纤维弯曲倾角、纤维直径、纤维尺寸密度之间的关系，带入空气的比热、密度和问题一计算得到的单根 A 纤维的比热和密度的乘积建立具体的关系式。建立两种整体分析模型和微观最小单位模型的热传导方程和边界条件，将方程和条件进行差分，将残差平方和和整体织物的热导率作为目标函数，建立双目标优化模型求解最小的整体导热率，建立约束条件，利用分层序列的思想先把残差平方和作为目标，再将整体织物的热导率作为目标，对纤维弯曲倾角、纤维直径、纤维尺寸密度进行遍历搜索。求得最优的整体织物热导率为0.0317。
  对于问题三，对两问建立的热传导方程更改厚度为0时的边界条件，利用牛顿冷却定律建立织物表层和空气接触面的边界条件，利用问题一和问题二的模型建立部分参数的计算方法，利用有限差分法将边界条件离散化，建立向后隐式差分格式，利用数值积分的方法对热传导方程进行求解，新问题一求解得到的最优残差平方和为47.6978，单根纤维A的热导率为0.8335。新问题二求解得到的最优整体织物热导率为0.0306。

模型假设：

  1、为了方便计算，将圆柱体的边界看成圆弧和直线的简易几何组合。
  2、由于对于编织物来说，深度方向的热传导应远大于横向的热传导，故忽略横向热传导。
  3、由于空气密度，比热在温度变化较小的范围内变化较小，故直接将空气的密度和比热看成常数。

问题分析：

  本文是对于隔热材料进行考虑，即要求使用热力学相关知识来解决材料问题，在已知附件一的温度下考虑织物的热导率以及单根纤维的热导率之间关系。
  问题一
  对于问题一，在已知热源侧织物的表面温度要求建立平纹织物整体热导率与单根纤维热导率之间关系的数学模型[1]。首先，根据题意热传导发生在空气之间、空气与织物之间以及织物内部单根纤维之间，所以要想求得整体织物热导率与单根纤维热导率之间关系的数学模型，我们需要考虑纤维传热和空隙间气体传热。在有上述分析后，由于静态空气热导率，经向热导率以及纬向热导率大小不一样，所以需要进行考虑，但是无论哪种材料，对于本文的作用效果都是为了传热，所以我们可以对于本文织物空间中的温度进行考虑。我们可以分别对于整体织物进行热传递考虑，以及对于局部纤维之间进行热传递考虑，再分别观察两种方程下温度分布，将问题转化为优化类问题，即要对单根纤维热导率进行遍历，使得局部考虑下对比整体考虑下，织物空间温度分布大小差值最小。温度差值最小情况下，得出对应织物整体热扩散率与单根纤维热扩散率之间的关系。在有热扩散率之间的关系后，由于空间既有空气分布又有纤维材料分布，导致密度分布不均匀，所以要想求得热导率之间的关系，还需要得出密度关系，再根据比热的定义式，计算出比热，最终得出本题所要求的热导率之间的关系模型。得到关系模型后，再结合附件2 的样品参数以及整体热导率，得出最后要求的单根纤维的热导率。
  问题二
  对于问题二，仍然将问题视作最优化问题，即对于优化问题最后的目标是整体热导率最低，局部考虑热传递作用和整体考虑热传递作用结合第一问的整体热导率与单根纤维热导率之间的数学模型。问题中纤维的直径以及织物的经密、纬密、弯曲角度，这些因素会对密度以及比热造成影响[2]，从而热扩散率换算成热导率会受到影响。结合热传导方程以及文章的约束条件对密度和比热进行考虑。
  问题三
  对于问题三进行分析。问题三对于前两问，附件1的温度变成热源侧织物表面空气的温度，即此时对于热传导方程中边界条件以及初始条件发生了变化，所以需要对于前两问的热传递模型进行改变，其余分析步骤几乎一致。

模型的建立与求解整体论文缩略图

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可

部分程序代码：

1. clear   
2. clc   
3. load('u1.mat')   
4. load('ua.mat')   
5. as=0.663*10^-6;%s:整体 a:空气 t:纤维   
6. ks=0.033;   
7. cs=0.05*10^6;   
8. aw=800;   
9. aj=600;   
10. d=0.6/1000;   
11. ow=26.54*pi/180;   
12. oj=19.8*pi/180;   
13. pa=1.1691;   
14. ca=1012;   
15. ka=0.0296;   
16. B=100000000;   
17. % ps=ks/(as*cs);   
18. % ps=1.5*10^3;   
19. % ms=2*d*ps/(aw*aj);   
20. vw=((1/aj-3*d*ow)/cos(ow)+3*d*ow)*pi*d^2/4;   
21. vj=((1/aw-3*d*oj)/cos(oj)+3*d*oj)*pi*d^2/4;   
22. vs=2*d/(aw*aj);   
23. ma=(vs-vw-vj)*pa;   
24. % pt=(ms-ma)/(vw+vj);   
25. % mt=pt*(vw+vj);   
26. % ct=ms*(cs-(ca*ma)/ms)/mt;   
27. aa=ka/(ca*pa);    
28. ptct=(ks*vs-as*ca*ma)/(vj+vw)/as;   
29. load('datatu.mat')   
30. fai=datatu(:,2);   
31. uc=zeros(11,641);   
32. u=zeros(11,641);   
33. U=zeros(11,641);   
34. deltaz=1.2/1000/10;   
35. deltat=0.00015625;   
36. t1=0:0.02:0.1;   
37. t2=0:0.01:0.1;   
38. t3=0:0.005:0.1;   
39. t4=0:0.0025:0.1;   
40. t5=0:0.00125:0.1;   
41. t6=0:0.000625:0.1;  
42. t7=0:0.0003125:0.1;   
43. t=0:0.00015625:0.1;   
44. z=0:deltaz:(1.2/1000);   
45. interfai=interp1(t1,fai,[0.01,0.03,0.05,0.07,0.09],'spline');   
46. fai1=[];   
47. for i=1:5   
48.       fai1=[fai1,fai(i),interfai(i),];   
49. end   
50. fai1=[fai1,fai(6)];   
51. interfai1=interp1(t2,fai1,[0.005 0.015,0.025,0.035,0.045,0.055 0.065 0
 .075 0.085 0.095],'spline');   
52. fai2=[];   
53. for i=1:10   
54.       fai2=[fai2,fai1(i),interfai1(i),];   
55. end   
56. fai2=[fai2,fai1(11)];   
57. interfai2=interp1(t3,fai2,[0.0025 0.0075,0.0125,0.0175,0.0225,0.0275 
0.0325 0.0375 0.0425 0.0475 0.0525 0.0575 0.0625 0.0675 0.0725 0.0775 0.08
 25 0.0875 0.0925 0.0975],'spline');   
58. fai3=[];   
59. for i=1:20   
60.       fai3=[fai3,fai2(i),interfai2(i),];   
61. end   
62. fai3=[fai3,fai2(21)];   
63. interfai3=interp1(t4,fai3,0.00125:0.0025:0.1,'spline');   
64. fai4=[];   
65. for i=1:40   
66.       fai4=[fai4,fai3(i),interfai3(i),];   
67. end   
68. fai4=[fai4,fai3(41)];   
69. interfai4=interp1(t5,fai4,0.000625:0.00125:0.1,'spline');   
70. fai5=[];   
71. for i=1:80   
72.       fai5=[fai5,fai4(i),interfai4(i),];   
73. end   
74. fai5=[fai5,fai4(81)];   
75. interfai5=interp1(t6,fai5,0.0003125:0.000625:0.1,'spline');   
76. fai6=[];   
77. for i=1:160   
78.       fai6=[fai6,fai5(i),interfai5(i),];   
79. end   
80. fai6=[fai6,fai5(161)];   
81. interfai6=interp1(t7,fai6,0.00015625:0.0003125:0.1,'spline');   
82. fai7=[]; 
83. for i=1:320   
84.       fai7=[fai7,fai6(i),interfai6(i),];   
85. end   
86. fai7=[fai7,fai6(321)];   
87. uc(:,1)=25+273;   
88. uc(1,:)=fai7'+273*ones(641,1);   
89. v=zeros(11,641);   
90. for l=0:0.01:2   
91.     kt=l;   
92.     v=zeros(11,641);   
93.     at=kt/ptct;   
94.     ra=aa*deltat/(deltaz^2);   
95.     rt=at*deltat/(deltaz^2);   
96. for tt=1:6   
97.     for i=1:640   
98.         if tt==1||tt==2||tt==6||tt==5    %介质全为纤维的情况   
99.             uc(2:10,i+1)=(1-2*rt)*uc(2:10,i)+rt*(uc(1:9,i)+uc(3:11,i));   
100.             uc(11,i+1)=uc(10,i+1);   
101.         else   
102.             dastart=tt+1;daend=dastart+5;   
103.             %%空气内部的计算过程   
104.             uc(2:dastart,i+1)=(1
2*ra)*uc(2:dastart,i)+ra*(uc(1:dastart-1,i)+uc(3:dastart+1,i));   
105. %           %%边界条件   
106. %           u(dastart+1,i+1)=ka/kt*(2*u(dastart,i+1)-u(dastart
1,i+1));   
107.             %%材料内部计算过程   
108.             uc(dastart+2:daend,i+1)=(1
2*rt)*uc(dastart+2:daend,i)+rt*(uc(dastart+1:daend
1,i)+uc(dastart+3:daend+1,i));   
109.             %%空气内部的计算过程   
110.             uc(daend+1:end-1,i+1)=(1-2*ra)*uc(daend+1:end
1,i)+ra*(uc(daend:end-2,i)+uc(daend+2:end,i));   
111.             %%最后一项   
112.             uc(end,i+1)=uc(end-1,i+1);   
113.         end   
114.     end   
115.        v=v+uc;   
116. end   
117.    u=v/6;   
118. %    u=0.84*uc+0.16*ua;   
119.    U=(u-u1).^2;   
120.    b=sum(sum(U));   
121.     if B>b  
122.         B=b;   
123.         A=at;   
124.         K=kt;   
125.         g=u;   
126.     end   
127. end   
128. disp(['寻找到的最佳位置为：',num2str(K)])   
129. disp(['寻找到的最佳位置为:',num2str(A)])   
130. g=g-273*ones(11,641);   
131. mesh(t,z,g)   
132. title('微观分析的温度分布')   
133. xlabel('时间')   
134. ylabel('厚度')   
135. zlabel('温度')

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可

原文地址：https://blog.csdn.net/weixin_43292788/article/details/144179303

免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！