2021年第十届数学建模国际赛小美赛
A题 气道阻力的评估
原题再现：
气道阻力的定义是通过肺气道产生单位气流所需的经肺压力的变化。更简单地说，它是嘴和肺泡之间的压力差，除以气流。影响气道阻力的因素是多方面的，我们需要探讨这些因素是如何影响气道总阻力的。
许多研究假设肺气道是一个光滑的相邻圆柱体系统，空气以层流方式流经该系统，其结果基于稳定流动的假设。这种假设与实际情况有许多不一致之处。例如，鼻气道阻力占总气道阻力的50%以上，但鼻腔不是圆柱形的。这种假设完全忽略了上气道的紊流。因此，预测气道的总通气量是不够的。
我们需要建立一个模型来估计总气道阻力至少由这些因素：气流速度，肺容量，以及气道各部分的横截面积。这些是一些最重要的促成因素。对于某些患者，由于手术等原因，气道的某些部分会暂时/永久改变，包括某些部分的形状或横截面积。我们的模型将用于评估患者气道的通气性能，并帮助医生设计对气道总通气性能影响最小的治疗方案。
整体求解过程概述(摘要)
气道阻力定义为产生单位气流通过肺气道所需的跨肺压力变化。气道阻力的大小间接反映了患者的气道通气性能。因此，有必要找出影响气道总阻力的因素。
在问题的第一部分，根据流体力学方程，得到了口与肺泡压差、气道长度、气道各截面截面积等因素之间的关系，进而得到气道总阻力的表达式。由于气体可能具有一定的粘度，通过计算不考虑气体粘度的伯努利方程和分别考虑气体粘度的伯努利方程，可以得到压差与其它因素的关系。计算气道各节段的气流速度和口肺泡与肺的压差。在具体的计算中，采用蒙特卡罗模拟的方法，计算出气道各截面的一些物理量。最后根据气道阻力的定义计算气道阻力。在不考虑气体粘度和考虑气体粘度的情况下，利用MATLAB绘制了两种情况的比较图。通过比较，最终得到考虑气体粘度条件下的气道阻力。
在问题的第二部分，根据呼吸系统的压力-容积曲线（拟合结果为s型曲线），压差与肺容积之间的关系可用指数函数来描述。然后，利用上一问题中的方法和MATLAB计算了气道内的平均气体流量。通过定义气道阻力，得到气道阻力与肺容量的关系曲线（其他变量在此控制不变）。
并利用MATLAB控制其他变量，绘制气道阻力与气道各截面面积的关系曲线。通过观察关系曲线，可以得出对全气道通气性能影响不大的治疗方案。
在此基础上，总结了气道阻力的影响因素和物理表现，并得到了相应的微分方程。
模型假设：
（1） 由于气体的粘度非常小，因此未考虑气体粘度对气道阻力估计值的影响。
（2） 当倾斜时，气体重力引起的倾斜角对压力的影响被忽略。
（3） 由于人气道的每一段的横截面积可以变化，因此假设人气道的表面是平滑的曲面（不一定是圆柱形的），并且气道的每一段的横截面积满足正态分布。
（4） 当气道的段数足够大时，假设每个段中的压差与每个段的长度成线性关系，即从气道的一个段到嘴的压差与它们之间的距离（深度）成线性关系。
（5） 其他阻力对气道阻力估计的影响被忽略。
（6） 假设气体密度不随气道深度变化，且气体密度p=1.293kg/m3为常数。
问题重述：
气道阻力可能受到许多不确定因素的影响，如气流速度、肺容量和每个气道段的横截面积。所以我们需要讨论这些因素是如何影响总气道阻力的。许多研究假设肺气道是一个光滑的连接柱体系统，空气以层流方式流过该系统，其结果基于稳定流动的假设。这种假设与实际情况有许多不一致之处。例如，鼻气道阻力占总气道阻力的50%以上，但鼻腔不是圆柱形的。这种假设完全忽略了上气道的紊流。因此，我们需要建立一个模型来描述影响气道阻力的因素。
1气流速度与气道阻力的关系
气道速度和气道阻力之间的关系是什么？如何利用曲线分析这些因素对气道阻力的影响？
为了解决这类问题，我们需要找到一个合适的物理模型来描述这些物理量之间的关系。由于该问题是以气体（属于流体）为基础的，因此需要用流体力学方程来描述这些物理量之间的关系·首先，通过对伯努利方程进行变形，得到了相应的微分方程形式和气道各段气速、压差的递推公式。当然，也可以得到肺泡附近的气流速度与肺泡与口部的压差和初始气流速度之间的关系。然后，利用递推公式在MATLAB中进行蒙特卡罗仿真，并进行一定次数的迭代，使气道的每个截面都能近似光滑圆柱体，从而使数值接近真实值。最后，通过MATLAB绘图分析气道阻力与其他因素的关系。
·由于以前没有考虑空气粘度对气道内其他物理量的影响，本文对伯努利方程进行了修正，采用粘性流体伯努利方程求解。重复刚才的操作，在MATLAB上绘制图形，比较考虑粘性和不考虑粘性两种情况，最后选择考虑粘性的Bernoulli方程方案求解。
2肺容量与气道阻力的关系
肺容量和气道阻力之间的关系是什么？
由于气道阻力是嘴和肺泡之间的压差与气道中气流的比率，因此问题可分为两部分：
1） 探讨肺动脉压差与肺容量的关系。
2） 通过以上两种关系间接探讨气道阻力与肺容量的关系。因此，通过寻找信息，我们找到了适合描述肺容积与压差关系的曲线：呼吸系统压力容积曲线。
3横截面积与气道阻力的关系
气道各部分的横截面积与气道阻力之间的关系是什么？
根据第一个主题的方法，通过控制其他变量不变，探索气道阻力与气道各部分横截面积之间的关系。利用MATLAB软件分析了两者之间的关系。
4如何评估患者的气道性能
通气性能通常与气道阻力有关。因此，有必要探讨通气性能与气道阻力的相关性，并通过得出的结论对通气性能进行评价。
5对全气道通气性能影响最小的方案
对于某些患者，由于手术或其他原因，气道的某些部分可能会暂时/永久改变，包括某些部分的形状或横截面积。因此，可将该问题转化为以下问题：选择与术前横截面积差值最小的合适横截面积，使该面积下的气道阻力差值小于术前。
通过分析第三个问题中气道阻力与横截面积的关系曲线，可以确定引起气道阻力微小变化的横截面积（横截面积应尽可能接近原始横截面积）。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
部分程序代码：(代码和文档not free)
function mean_p_airway=non_viscosity(dp,H,N)
v0=14.8e-6;
g=9.8;
ru=1.293;
p_airway=zeros(1,length(N));
for i=1:length(N)
v2=[v0^2,zeros(1,N(i))];
S=zeros(1,N(i));
for j=1:N(i)
S(j)=(0.0125^2*normrnd(1,rand(1)))*1e-2;
v2(j+1)=2/ru*(dp/N(i)-ru*g*H/N(i))+v2(j);
end
p_airway(i)=dp/mean(abs(sqrt(v2(2:end))).*S);
end
mean_p_airway=mean(p_airway);
end
function mean_p_airway=airway_resistance(dp,H,N0)
v0=14.8e-6;
g=9.8;
ru=1.293;
mu=0.01894;
p_airway=zeros(1,length(N0));
for N=N0
v2=[v0^2,zeros(1,N)];
S=zeros(1,N);
for i=1:N
S(i)=(0.0125^2*normrnd(1,rand(1)))*1e-2;
v2(i+1)=2/ru*(dp/N-ru*g*H/N-8*pi*mu*sqrt(v2(i))*H/(N*g*S(i)));
end
p_airway((N-N0(1))/(N0(2)-N0(1))+1)=dp/mean(abs(sqrt(v2(2:end))).*S);
end
mean_p_airway=mean(p_airway);
end
clear;clc;close all
p0=1.013e5;
v0=14.8e-6;
g=9.8;
ru=1.293;
mu=0.01894;
lN=input('Please enter the lower limit of monte Carlo simulation:');
uN=input('Please enter the maximum number of Monte Carlo simulations:');
N=lN:uN;
dp=100:300;
airway=zeros(1,length(dp));
mean_p_airway=zeros(1,length(dp));
for i=1:length(dp)
airway(i)=non_viscosity(dp(i),0.12,N);
mean_p_airway(i)=airway_resistance(dp(i),0.12,N);
end
figure(1)
plot(dp,airway,'r-p',dp,mean_p_airway,'b-*')
xlabel('differential pressure/Pa')
ylabel('airway resistance/Pa*m^{-3}*s')
legend('without the viscosity of the air','with the viscosity of the air')
figure(2)
plot(dp,airway,'r-p')
xlabel('differential pressure/Pa')
ylabel('airway resistance/Pa*m^{-3}*s')
figure(3)
plot(dp,mean_p_airway,'b-*')
xlabel('differential pressure/Pa')
ylabel('airway resistance/Pa*m^{-3}*s')
H=0.11:0.001:0.13;
airway=zeros(1,length(H));
mean_p_airway=zeros(1,length(H));
for i=1:length(H)
airway(i)=non_viscosity(200,H(i),N);
mean_p_airway(i)=airway_resistance(200,H(i),N);
end
figure(4)
plot(H,airway,'r-x',H,mean_p_airway,'-*')
xlabel('The length of the airway/m')
ylabel('airway resistance/Pa*m^{-3}*s')
legend('without the viscosity of the air','with the viscosity of the air')
figure(5)
plot(H,airway,'r-x')
xlabel('The length of the airway/m')
ylabel('airway resistance/Pa*m^{-3}*s')
figure(6)
plot(H,mean_p_airway,'-*')
xlabel('The length of the airway/m')
ylabel('airway resistance/Pa*m^{-3}*s')