2020年华数杯数学建模A题带相变材料的低温防护服御寒仿真模拟求解全过程文档及程序

2020年华数杯数学建模

A题 带相变材料的低温防护服御寒仿真模拟

原题再现：

  在一些特定的场合，人们往往需要在极寒天气下作业，如高山高原工作、潜水员水下工作、现代化工厂的低温车间以及寒冷气候下的野外作业等。为了能使工作顺利进行，科学家们一直在研究低温防护复合材料，试图做成防护服用以保护在超低温环境下的工作者。
  某研究所研制的低温防护复合材料：三层结构，包括内层织物层、中间层功能层、外层隔热层。内层织物层主要用于舒适性。中间层是一种特殊的材料，可以产生并释放热量，用以延缓人体温度过快降低，称为相变材料。外层隔热层主要是延缓热量对外过快传递。低温防护材料主要用于短时间的低温防护，有效降低外界环境对人体的伤害。
  为了延缓人体温度过快降低，研制的复合材料的中间层要具有良好的保温性能。中间层有两个特性，特性一是厚度不能大于 0.45mm，因为中间层的硬度与厚度成正比，一旦超过 0.45mm，人体将无法伸展，也就无法工作。特性二是在高于 25℃左右（根据材料不同这个临界点会有小的变化）为液态，低于 25℃开始固化，固化时就开始放热，一直到 14.7℃左右固化完毕，将不再放热。具体数据详见附件 1。注意：附件 1 中的数据放热温度范围与上面表述温度范围有差异，以附件数据为准。
  热量传递方式有对流、辐射和传导三种。但在超低温下，对外辐射微不足道，因此一般不予考虑。内层织物与人体表面之间有空气流动，外层隔热层与外部环境之间也有空气流动。空气流速不一样，所计算出来的表面换热系数也会不一样。至于热传导能力（热导率），是材料的物理属性。附件 2 给出了三层材料的物理属性值。
  为检验这种复合材料的耐低温效果，科研者按照附件 2 提供的厚度为一名身高 1.70m，体重为 60kg 的中国实验者(消耗的衣料面积一般不超过人体表面积的1.25 倍)制作了一套耐低温服装。实验者将前往南极洲长城站在-40℃的低温下进行工作实验。一般来讲，人体表面温度低于 15℃工作将很困难，低于 10℃将有生命危险。科研者为慎重起见，希望在去南极洲实验前先做仿真模拟。 。
  请你查询和阅读相关资料，运用数学模型，回答下列问题：
  问题一：假设实验者站在南极洲长城站附近。当时环境安静，晴天无风，-40℃。请你利用传热学原理对低温防护复合材料在低温环境下传热机理进行分析，建立热量传递数学模型。并预测实验者能够在室外坚持多久？
  问题二：假设长城站外面风速为 3m/s，实验者做轻微运动。如何改进第一问的模型，并计算实验者大约多久后就必须要返回室内？
  问题三：在问题一的室外环境下，实验者静止不动，实验者穿着较重的防护服，最大承受重量为 100kg。并且每多承受 10s 的时间，最大承受重量将下降0.5kg。如果实验室在低温防护服材料原有支出的基础上能够增加 50%以内的资金，请问如何增加防护服的厚度，使得实验者在外站立时间尽可能的长?
  问题四：接着问题三，假定不追加资金，要想提高在户外活动时间，可以通过提高复合材料中间层的放热能力来实现，这是一个新的研究课题。那么，中间功能层的放热能力要比原先提高多少，才能使得实验者在室外坚持的时间不比第三问中坚持的时间少？假定放热能力的提高是在各个温度下同比例增加的。

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

  本文针对带相变材料的低温防护服的御寒效果进行研究，综合考虑多种传热方式建立热传导模型，并以此模型为基础对低温防护服的防护效果以及一定资金条件下的最佳制作厚度问题进行了研究。
  对于问题一，本文首先利用传热学原理对实验者、低温防护服以及外界空气环境组成的系统建立了一维热传导偏微分方程模型，并确定了定解的边值条件，研究了在无风条件下带相变材料的低温防护服对静止的实验者的防护效果。本文利用有限差分方法对微分方程进行转化求解，结果显示：实验者在室外能够坚持的时长为675𝑠。
  对于问题二，本文在问题一模型的基础上，考虑风速（3𝑚/𝑠）对对流换热系数的影响，由经验公式给出了计算对流换热系数时的修正项，并用修正后的模型研究实验者做轻微运动时低温防护服的防护效果。结果显示：实验者在室外能够坚持的时长为152𝑠。
  对于问题三，本文研究了在多重限制条件（资金、最大承受重量、中间层厚度）下，应该如何增加防护服的厚度，使得实验者在外站立时间尽可能长。本文对此建立了目标优化模型，以实验者的在外站立时长为优化目标，以资金数额和最大承受重量为约束条件，通过编程求解该非线性优化模型的最优解，并得到低温防护服的增厚方案。结果显示：保持内层（织物层）厚度不变，中间层（功能层）的厚度增加至0.45𝑚𝑚 ，最外层（隔热层）的厚度增加至1.2𝑚𝑚时御寒效果最优，此时能在外坚持的最长时长为872.165𝑠。
  对于问题四，本文以问题三中防护服的御寒能力为研究参照，对新复合材料中间功能层的放热能力给出了预期。本文应用二次搜索法寻找两种材料放热能力的比值𝛼，结果显示：如果中间功能层的放热能力比原先提高37.9%，用其制作的低温防护服的御寒能力将不弱于第三问中低温防护服的御寒能力。

问题分析：

  为了模拟实验者站在南极洲长城站附近时的热量散失过程，我们需要对实验者、带相变材料的低温防护服以及外界空气环境组成的系统建立一维热传导数学模型，并应用此模型预测实验者在室外能够坚持的时长，并以此评估在无风条件下带相变材料的低温防护服对静止实验者的防护效果。
  在问题一模型的基础上，我们要对有风和实验者运动的情况进行研究。在计算对流换热系数时，我们需要考虑由风速（3𝑚/𝑠）产生的修正项，并用修正后的模型研究实验者做轻微运动时低温防护服的防护效果。
  对于问题三，对非线性优化模型进行求解，在多重限制条件（资金、最大承受重量、中间层厚度）下，考虑应该如何增加防护服的厚度，使得实验者在外站立时间尽可能长。我们需要建立目标优化模型，并把实验者的在外站立时长作为优化目标，把资金上限和最大承受重量作为约束条件，通过编程求解该问题的最优解，并由此得到低温防护服的最佳增厚方案。
  对于问题四，我们需要把问题三中防护服的御寒能力作为研究参照，求解复合材料中间层的放热能力和在外活动时间的关系。由此推算出中间功能层的放热能力应比原先提高多少，才能使实验者在室外坚持的时间不比第三问中坚持的时间少。

模型假设：

  假设 1：在超低温的环境中，我们不考虑物体对外辐射的热量传递方式
  假设 2：耐低温防护服的三层材料之间的间隙可以忽略不计
  假设 3：在有热对流现象的气层中空气是均匀流动的
  假设 4：材料放热能力的提高是在各个温度下同比例增加的

模型的建立与求解

热传导偏微分方程模型

  由假设可知，耐低温服装各处的温度与热量交换情况完全一致。因此，我们沿着垂直于皮肤的方向将整个服装系统的热量交换抽象为一维情况，如下图所示：

  热传导就是如果空间某物体内各点处的温度不同，热量会从温度较高的点处向温度较低的点处流动的现象。我们先来推导一般热传导方程，再根据低温防护服每一层的情况建立相应的热传导方程、边值条件与初值条件。

有限差分法[2]求解

  对于以上微分方程组，求出解析解是十分困难且没有必要的，因此我们采用有限差分法将微分方程转换为相应的差分方程，从而递推求出数值解。

  在编程求解前，为得到所需温度点处的放热能力，我们对附件 1 中的数据做三次样条插值，使得放热能力-温度成为一条通过一系列形值点的光滑曲线。然后对该有限差分方程组使用 Python 语言进行编程求解。得到温度随时间的变化关系如下图所示：

  以体表温度降低到 10℃作为坚持极限，由上图得知：实验者静止不动，环境无风的情况下能够坚持的时长为675𝑠。
  假设长城站外面风速为 3𝑚/𝑠，其大小处于受迫对流的风速区间。根据国家标准GB8175-87，这种情况下的对流换热系数计算可用如下经验关系式计算：



  由于最外层（隔热层）图层每层制作费用为44.6067元，故最多可以增加三层最外层涂层。通过编程分别求解在防护服最外层涂层增加 1、2、3 层时，温度随时间的变化关系如下图所示：

  在问题三的背景下，如果不追加资金，要想提高实验者在户外的活动时间，我们可以通过提高复合材料中间层的放热能力来实现。对于由新材料制作的防寒服，中间层（功能层）应成立：

论文缩略图：

程序代码：

calculation.py

import numpy as np 
import pandas as pd 
from scipy.interpolate import interp1d 
thickness = [0.3, 0.4, 0.7] 
rou = [300, 552.3, 208] 
c = [5463.2, 2400, 4803.8] 
k = [0.0527, 0.06, 0.068] 
v = 20 
T_out = [37] * int(thickness[0] * v + 1) 
T_mid = [37] * int(thickness[1] * v + 1) 
T_in = [37] * int(thickness[2] * v + 1) 
T_env = -40 
T_peo = 37 
target_columns = ['time', 'T'] 
t_dict = {} 
filepath="./attachment1.xlsx" 
data = pd.read_excel(io=filepath,sheet_index=0,header=0) 
g = data.values[:, 0] 
f = data.values[:, 1] 
inter = interp1d(g, f, kind="cubic") 
g = np.arange(g.min(), g.max(), 0.001) 
f = inter(g) 
g = [round(i,3) for i in g] 
def get_heat(t, g, f): 
 if t > max(g) or t < min(g): 
 return 0 
 index =int((t-14.05)*1000) 
 s = f[index] 
 f[index] *= 1 
 return s 
delta_x = 0.05 * 0.001 
delta_t = 0.035 
T = [0] * 200000
sign = 0 
afa = 1.379 
for t in np.arange(0, 200000): 
 h1 = 2.38 * (T_out[0] - T_env) ** 0.25 
 #h1 = 14 
 h2 = 3 
 r = k[0] / (c[0]*rou[0]) * delta_t / (delta_x**2) 
 #print(r) 
 if T_out[0] > -40: 
 T_out[0] = 2*r*T_out[1] + (1 - 2*r- 2*r*delta_x*h1/k[0])*T_out[0] + 
2*r*delta_x*h1*T_env/k[0] 
 if T_out[0] < -40: 
 T_out[0] = -40 
 
 for i in range(1, int(thickness[0] * v)): 
 T_out[i] = r * (T_out[i+1] + T_out[i-1]) + (1 - 2*r) * T_out[i] 
 if T_out[i] < -40: 
 T_out[i] = -40 
 
 
 T_out[-1] = (k[0] * T_out[-2] + k[1] * T_mid[1]) / (k[0] + k[1]) 
 T_mid[0] = T_out[-1] 
 
 
 r = k[1] / (c[1]*rou[1]) * delta_t / (delta_x**2) 
 for i in range(1, int(thickness[1] * v)): 
 s = T_mid[i] 
 T_mid[i] = r * (T_mid[i+1] + T_mid[i-1]) + (1 - 2*r) * T_mid[i] - afa * 
get_heat(round(T_mid[i], 3), g, f) * 1000 * delta_t / c[1] 
 if s <= T_mid[i]: 
 T_mid[i] = s 
 if T_mid[i] < -40: 
 T_mid[i] = -40 
 
 T_mid[-1] = (k[1] * T_mid[-2] + k[2] * T_in[1]) / (k[1] + k[2]) 
 T_in[0] = T_mid[-1] 
 
 
 r = k[2] / (c[2]*rou[2]) * delta_t / (delta_x**2) 
 
 for i in range(1, int(thickness[2] * v)): 
 s = T_in[i] 
 T_in[i] = r * (T_in[i+1] + T_in[i-1]) + (1 - 2*r) * T_in[i]
  if s <= T_in[i]: 
 T_in[i] = s 
 if T_in[i] < -40: 
 T_in[i] = -40 
 
 T_in[-1] = 2*r*T_in[-2] + (1 - 2*r- 2*r*delta_x*h2/k[2])*T_in[-1] + 
2*r*delta_x*h2*T_peo/k[2] 
 print(T_in[-1]) 
 T[t] = T_in[-1] 
 if T_in[-1] <= 14 and sign == 0: 
 sign = 1 
 t2 = t 
 if T_in[-1] <= 10: 
 print("-----------hard-------------") 
 print(t2) 
 print("-----------threat-------------") 
 s = t 
 print(s) 
 break 
t = np.arange(s+1) * delta_t 
t_dict[target_columns[0]] = t[:] 
t_dict[target_columns[1]] = T[:t.size] 
t_generation = pd.DataFrame(t_dict) 
t_generation.drop_duplicates(subset=None, keep='first', inplace=True) 
t_generation.to_csv("./T5.csv", index=None) 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101

plot.py

import pandas as pd 
from matplotlib import pyplot as plt 
from matplotlib import font_manager 
my_font = font_manager.FontProperties(fname="C:\Windows\Fonts\msyh.ttc") 
plt.figure(dpi=100) 
data = pd.read_csv('./T5.csv', sep=',') 
t = data.values[:,0] 
T = data.values[:,1] 
plt.plot(t, T) 
r1 = 24190 
r2 = 24919 
plt.scatter([t[r1],t[r2]],[T[r1],T[r2]],c = 'r') 
plt.text(t[r1]-190,T[r1],(round(t[r1],3),14),color='r') 
plt.text(t[r2]-210,T[r2],(round(t[r2],3),10),color='r')
plt.xlabel("时间 t", fontproperties=my_font) 
plt.ylabel("温度 T", fontproperties=my_font)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16