小提琴自动演奏机器人中的齿轮系统设计与制造

项目简介

小提琴作为一种历史悠久、表现力丰富的乐器，其演奏通常需要高度的技巧和长时间的训练。然而，自19世纪末以来，随着机械和电子技术的发展，人们开始尝试开发能够自动演奏小提琴的装置，旨在模拟甚至超越人类的演奏能力。这些自动演奏小提琴的出现，不仅是工程技术上的突破，也反映了人类对音乐自动化和完美演奏的持续追求。

小提琴自动演奏机器人是一个集精密机械、电子、控制和人工智能技术于一体的复杂系统，旨在模拟人类小提琴演奏家的精湛技艺，实现音乐的自动化演奏。该项目的核心在于通过精确控制弓弦摩擦、指板按压和琴身姿态，准确再现音高、音色、节奏和情感表达，将小提琴演奏的复杂动作分解为可编程的机械运动。在整个系统中，齿轮传动系统作为核心传动部件，其设计与制造的精度、稳定性与可靠性直接影响着机器人的演奏性能和音质。

本报告将深入探讨小提琴自动演奏机器人中齿轮系统的设计、材料选择、制造工艺以及相关的仿真与校核，旨在分析小提琴演奏对齿轮系统的特殊要求（如高精度定位、低噪音运行、快速响应），对比不同加工工艺的优缺点，从而得出最适合自动演奏机器人的齿轮制造解决方案。通过对齿轮传动系统的详细分析，旨在为提升小提琴自动演奏机器人的整体性能提供理论依据和实践指导。

项目背景及意义

项目进展

我们是如何过渡到这个项目的。

在上个学期，团队深入研究了小提琴教育行业，分析了小提琴产业链各环节的市场状况与主要问题，归纳出三个核心问题：第一，音乐教育市场缺乏创新——与十年前相比，教学模式和方法基本未变，但学习成本却增长了十倍；第二，小提琴制造业虽表面标准化，实则是非标准化生产，从生产到销售缺乏有效的管理和标准化指标，导致消费者需要付出极高成本才能购得理想的琴；第三，也是最核心的问题，音乐消费的价值观正在变化，音乐就像绘画一样，最开始在墙上画，后来在布料上面、在教堂的玻璃上面，材料和载体一直都在变化，而如今，越来越少的人会单纯为音乐消费，但同时，随着流媒体的传播和基础音乐教育的普及，音乐正逐渐成为生活场景中一个必不可少的元素。

上个学期在技术层面，我们从第一性原理出发，研究了小提琴的基础发声原理，并开发了基于共振原理的”小提琴音响”样机，通过利用小提琴的琴腔，能够让使用者直观感受到，音乐是由小提琴琴腔发出来的。以及不同材质对音乐风格的独特影响。这种效果与传统扬声器（Speaker）有着本质区别，每种材质的小提琴都能形成独特的”音乐滤镜”。

基于这项技术原理，我们展开产品的设计，基于学琴者家庭日常练琴的“反馈”缺失的痛点，设计了一款能够实施反馈的“智能练琴搭档”。

在寒假期间，一次偶然的发现，我们关注庞大的“闲置乐器”市场。2000 年至 2020年，是国内音乐教育的高速增长期，乐器，承载着千千万万家庭的音乐梦想走入千家万户，堪称一场现代”文艺复兴”。作为这场教育浪潮的亲身参与者，我们观察到年轻的乐器学习者随着年龄增长和学习压力加大，往往会逐渐搁置乐器。这个被忽视的闲置现象涉及庞大的市场规模和用户群体，恰恰为我们提供了机遇，本学期的多个项目都以此为出发点。

当前的核心方向

当前我们的工作聚焦于以下三个核心方向：

1. 更好地帮助人们学琴 —— 提高学习效率，降低学习成本；

2. 打造更优质的小提琴及配件 —— 提升产品体验与性能；

3. 帮助人们创造并享受音乐 —— 让音乐成为生活中有趣的部分。

项目背景

机器人技术在当今社会已不再局限于工业生产线，而是逐渐渗透到各个领域，并在艺术领域展现出其独特的潜力和广阔的应用前景。小提琴自动演奏机器人正是这一趋势下的产物，它代表了精密机械、控制系统、人工智能等多学科交叉融合的最新成果。小提琴演奏作为一门高度复杂的艺术形式，对演奏者的生理和心理素质提出了极高的要求，需要长期的专业训练才能掌握其精髓。然而，机器人技术的发展为突破人类生理极限、实现音乐的自动化演奏提供了可能。

项目意义

从用户需求出发

“让沉睡的琴声醒来。”在中国家庭中，数百万把价值不菲的小提琴正面临被束之高阁的命运——调查显示，超过60%的家庭小提琴在购买后三年内沦为摆设。我们的自演奏小提琴项目应运而生，我们通过智能演奏再生系统，让每一把被遗忘的小提琴重新焕发生机，成为家庭美育的活跃分子。

我们的自动演奏核心能让闲置的小提琴自主发声：清晨自动演奏《晨曲》唤醒全家，晚餐时即兴伴奏营造温馨氛围，睡前播放舒缓旋律助眠。通过智能琴弦驱动技术，无需人工干预，就能使小提琴持续输出高品质音乐，解决”买来就闲置”的痛点。对于仍有学习需求的家庭，我们考虑提供双模陪伴系统：在练习时，AI可同步示范标准演奏；在休息时，系统会自动播放经典曲目进行熏陶。特别设计的成长记忆功能会记录孩子的练习历程，定期生成”音乐成长报告”，让家长直观看到投资回报。

从技术到体验，我们的项目始终围绕一个核心：让每一把被冷落的小提琴重新获得宠爱。我们不只是提供硬件，更是提供一种全新的乐器使用理念——当您的小提琴能够自主创造价值时，闲置将永远成为过去式。

其他意义

• 推动技术进步：项目涉及精密机械设计、高性能驱动控制、复杂路径规划、实时传感反馈和人工智能算法等多个前沿领域，其研发过程将极大地推动这些交叉学科的技术进步。

• 创新音乐教育：作为一种辅助教学工具，小提琴自动演奏机器人可以帮助学生更好地理解和掌握演奏技巧，普及音乐教育，激发更多人对音乐的兴趣。

• 拓展艺术表现：机器人演奏为音乐创作和表演提供了新的可能性，艺术家可以利用机器人创造出人类难以实现的音效和演奏技巧，从而拓展艺术表现的边界。

• 文化传承与保护：通过数字化记录和机器人演奏，可以有效地传承和保护濒临失传的演奏技巧和音乐文化遗产。

• 人才培养实践：该项目为多学科交叉工程实践提供了宝贵的平台，有助于培养具备复合型知识和创新能力的优秀人才。

• 成本效益与普及性：随着技术的不断成熟和成本的降低，小提琴自动演奏机器人有望走向市场普及，让更多人能够接触和享受高品质的音乐演奏。

材料选择

小提琴自动演奏机器人中关键构件的材料选择至关重要，它直接影响着机器人的性能、精度、寿命和成本。在选择材料时，需要综合考虑其力学性能(如强度、硬度、韧性、耐疲劳性)、物理性能(如密度、热膨胀系数)、化学性能(如耐腐蚀性)、加工性能以及成本等因素。以下将详细阐述齿轮、运弓机构和按弦机构的材料选择及其理由。

齿轮材料选择

齿轮作为核心传动部件，其材料选择需满足高强度、高硬度、优异的耐磨性和抗疲劳性能，以确保在频繁启停和变速条件下保持稳定的性能和精度，并有效降低噪音。根据项目需求和常见齿轮材料特性，做出以下选择：

• 合金钢（Bow Controller）（齿轮、齿条）：20CrMnTi。这种材料经过渗碳淬火回火处理后，表面硬度可达到HRC58-62，心部仍保持高韧性。其优点在于强度高、耐疲劳性能优异，能够承受较大的载荷和冲击，适用于主传动系统，确保齿轮的高精度、高强度和长寿命。同时，其良好的耐磨性有助于延长齿轮的使用寿命。
• 工程塑料（Fingerboard traversal）（齿轮）：POM(聚甲醛)。工程塑料齿轮具有低噪音、自润滑、轻量化和耐腐蚀的优点。在对噪音要求较高或载荷较小的辅助传动系统中，可以考虑使用POM齿轮，以降低整体噪音并减轻重量。然而，其承载能力相对较低，不适用于主传动或高载荷场合。

对于小提琴自动演奏机器人这种对精度和稳定性要求极高的设备，主传动系统中的齿轮必须采用高强度合金钢20CrMnTi，以确保其在长期高频运转条件下具备出色的耐疲劳特性和精度保持能力。而对于一些次级传动或对噪音控制有更高要求的部位，可以辅助性地采用工程塑料齿轮，以平衡性能与成本，并进一步降低噪音。此外，精确的热处理(如淬火和回火)对于高强度合金钢齿轮至关重要，它能使其获得理想的表面硬度和耐疲劳性能。此外，经过查阅资料，对于噪音敏感的部件，还可以考虑在金属齿轮表面施加DLC(类金刚石)涂层，以提高耐磨性和降低摩擦噪音。

运弓（Bow Controller）构件材料选择

运弓机构需要实现轻量化、高刚度、多自由度控制以及平稳性与稳定性，以模拟人类运弓的精细动作。因此，材料选择应侧重于比强度、比刚度、阻尼性能和加工性。

• 碳纤维复合材料（主体部分）：碳纤维复合材料具有极高的比强度和比刚度，同时具备良好的阻尼性能，能够有效抑制振动。将其应用于运弓机构主体，可以显著减轻运弓机构的重量，从而减小运动惯性，提高系统的响应速度和控制精度。其优异的力学性能确保了运弓机构在高速往复运动中的稳定性。

• 铝合金 **(连接件、标准件)**：主要采用 7075、6061。铝合金具有重量轻、易于CNC加工的特点，成本相对较低。适用于运弓机构的连接件和支撑结构，可以在保证一定强度的前提下，平衡整体重量和制造成本。

运弓作为直接与琴弦接触并进行高速往复运动的关键部件，其轻量化和高刚度是实现高精度控制和快速响应的基础。碳纤维复合材料能够完美满足这些要求，使其成为运弓主体的理想选择。而铝合金则可以作为辅助材料，用于制造对刚度要求稍低但需要精密加工的连接部件，以实现整体性能和成本的优化。

按弦（Fingerboard traversal）构件材料选择

按弦机构需要实现高精度定位、快速响应和可调按压力，并模拟多指操作。因此，其材料选择需具备高强度、耐腐蚀、表面光洁度高和低摩擦等特性。

• 不锈钢 **(电机传动齿轮)**：例如304、316L。不锈钢具有优异的耐腐蚀性和高强度，表面易于抛光，能够获得较高的光洁度，从而降低摩擦系数。将其应用于按弦杆，可以保证按弦动作的长期稳定性和精度，并能有效抵抗琴板处手部汗液等可能造成的腐蚀。

• 工程塑料 **(模组滑轨)**：PEEK (聚醚醚酮)、UHMW-PE(超高分子量聚乙烯)。 PEEK 具有自润滑和耐磨特性，UHMW-PE则具有极低的摩擦系数。这些材料适用于按弦机构中的导套、滑块等部件，可以进一步降低摩擦阻力，提高机构的灵敏度和响应速度。

按弦杆需要频繁、精确地按压琴弦，因此材料的强度、耐磨性和表面光洁度至关重要。不锈钢能够提供所需的力学性能和表面特性，确保按弦动作的准确性和稳定性。而高性能工程塑料则可以作为辅助材料，用于制造需要低摩擦和自润滑特性的部件，以优化按弦机构的整体性能。

材料类型对比总结

综合来看，小提琴自动演奏机器人设计过程中，材料选择是一个多目标优化的过程，需要在性能、成本、加工性等方面进行权衡。通过合理搭配不同材料，可以最大限度地发挥各材料的优势，从而实现机器人的高性能、高精度和长寿命运行。

关键构件设计与建模

通过机械结构代替演奏者的左右手运动，实现小提琴的自动演奏功能，且演奏能够达到学习两年的演奏者水平。

将整琴分为运弓和按压两个模块，根据小提琴的四弦特性，又可将各模块分为单弦上的小模组。利用模块化设计，首先实现单弦的按压和拉弓功能，再将单个模组衍生至整把小提琴的自动演奏。

本节将详细阐述小提琴自动演奏机器人中关键构件的设计理念、详细尺寸参数及其设计理由，并展示其三维模型。

齿轮设计与建模

设计理念与理由

• 高精度传动：为确保运弓和按弦动作的精确性，齿轮传动系统必须达到极高的精度。因此，设计中采用ISO5级或更高精度的齿轮，以最大限度地减小传动误差，保证音准的准确性和一致性。
• 低噪音运行：齿轮啮合产生的噪音会直接影响演奏音质。为降低噪音，设计中优化了齿轮齿形，采用修形齿轮，并选择合适的模数和齿数，以减少啮合冲击和滑动摩擦。斜齿轮相较于直齿轮，具有渐进啮合的特点，能够显著降低冲击和振动，从而有效降低噪音。
• 紧凑型设计：考虑到机器人内部空间有限，设计中倾向于采用行星齿轮传动或谐波齿轮传动等紧凑型结构，以在有限空间内实现大传动比，满足机器人小型化和集成化的需求。
• 长寿命可靠性：频繁的演奏操作要求齿轮具备出色的耐疲劳特性。通过选用优质材料并结合适当的热处理工艺，确保齿轮具有足够的强度和耐磨性，以满足长期稳定运行的需求。

初代齿轮基本参数与计算

根据项目提供的参数，我们以一个模数(m)为 1mm，齿数(Z)为 60 的齿轮为例，进行主要尺寸参数的计算：

• 模数(m)：1mm

• 齿数(Z)：60

• 齿根半径：0.5mm

• 厚度：5mm

• 孔径：8mm

基于这些基本参数，可以计算出齿轮的其他关键尺寸：

• 分度圆直径(d)：d=m*z=1mm*60=60mm

• 齿顶高(ha)：ha=1mm

• 齿根高(hf)：hf=1.25*m=1.25*1=1.25mm

• 齿顶圆直径(da)：da=d+2*ha=60+2*1=62mm

• 齿根圆直径(df)：df=d-2*hf=60-2*1.25=57.5mm

参数选择说明：

• 齿数 60：选择 60齿的齿轮有助于实现紧凑型设计，并在保证传动比的同时，提供精细的运动控制能力。

• 齿根半径0.5mm：较小的齿根半径有助于减小齿根处的应力集中，从而提高齿轮的疲劳强度，延长使用寿命。

• 厚度 5mm：5mm的厚度旨在保证齿轮在承受载荷时的承载能力和整体刚度，防止变形。

• 孔径 8mm：8mm的孔径设计用于与电机轴进行匹配，并考虑采用键槽连接方式，以确保可靠的动力传输。

斜齿轮齿条设计需求与优化目标

针对小提琴自动演奏机器人对噪音和运动平稳性的高要求，斜齿轮齿条的设计尤为关键。直齿轮在高速运行时容易产生较大的冲击和噪音，而斜齿轮则能有效改善这一问题。

设计需求：

• 直齿轮冲击振动大，噪音峰值突出，不适用于本机器人系统。

• 速度需大于 120mm/s。

• 噪音需控制在 35dB 以下。

• 阻力需小于 5N。

• 要求往复运动具有高稳定性。

设计思路与优化目标:

• 螺旋角：选择15°～25°的螺旋角，以实现渐进啮合，从而降低啮合冲击力，减少噪音和振动。设计值设定为20°，优化目标是在轴向力与噪音之间取得最佳平衡。

• 齿向修形：通过对齿向进行修形，可以有效补偿偏载，使齿轮在实际运行中受力更加均匀，提高传动平稳性。

• 齿廓抛物线修型：对齿廓进行抛物线修型，可以消除啮合干涉，进一步优化啮合性能，降低噪音。

参数与设计值及优化目标：

通过上述设计理念和参数选择，旨在构建一个高精度、低噪音、长寿命的齿轮传动系统，为小提琴自动演奏机器人提供稳定可靠的动力输出。

运弓（Bow Controller）构件设计与建模

运弓构件是模拟人类运弓动作的关键部件，其设计需兼顾轻量化、高刚度、多自由度控制和运动平稳性。

初代拉弓模组存在“环形弓毛材料选择难”、“整体质量大”、“多弦拓展性不足”的问题，根据上述问题，对拉弓模组进行迭代，决定仍然将琴弓保留。

设计需求与理念

• 轻量化与高刚度：采用碳纤维复合材料和优化结构设计，旨在减小运弓的惯性，从而提高其响应速度和运动精度。轻量化设计有助于降低驱动系统的负担，提高动态性能。

• 多自由度控制：运弓设计具备三个自由度：往复运动（模拟弓速）、压力控制（模拟弓压）和角度调整（模拟弓毛与琴弦的接触角度）。这使得机器人能够模拟人类运弓的复杂技巧，实现对音色和音量的精细控制。

• 平稳性与稳定性：为确保运弓运动平稳无颤动，设计中采用高精度导轨和直线电机或丝杠传动。这有助于消除运动过程中的抖动，保证音质的纯净性。

• 模块化设计：弓毛部分采用可拆卸更换的模块化设计，方便日常维护和根据不同音色需求进行快速更换。

结构设计参数与材料选取

详细尺寸参数：

• 运弓总长：300-400 mm，通过小提琴尺寸和演奏范围确定。

• 弓压范围：0-5 N，可通过弹簧调节，以模拟不同力度下的弓压。

• 导轨行程：200-300 mm 可调，确保运弓在琴弦上的有效运动范围。

• 响应频率：≥100 Hz，保证运弓能够快速响应音乐节奏的变化。

材料配置：

• 运弓主体：碳纤维复合材料，采用工字型或箱型截面设计，以最大化刚度重量比。

• 连接件：7075 铝合金，通过 CNC 精密加工，确保连接精度和强度。

• 导轨系统：高精度直线导轨，重复定位精度可达±0.01mm，保证运弓运动的精确性和稳定性。

设计理由：

• 轻量化设计：显著提高运弓的动态响应性能，减少运动惯性，使得机器人能够更灵活地进行运弓操作。

• 多自由度：使机器人能够模拟人类运弓的复杂技巧，实现对音色、音量和节奏的精细控制，提升演奏表现力。

• 高精度导轨：确保运弓运动的平稳性和重复定位精度，对于音质的稳定性和一致性至关重要。

• 模块化结构：便于弓毛部分的维护和更换，同时也为未来的功能扩展提供了便利 。

按弦（Fingerboard traversal）构件设计与建模

按弦机构是模拟人类手指按弦动作的关键部件，其设计要求高精度定位、快速响应、可调按压力和多指模拟。

从技术角度而言，机器人需要能够在连续域内随时间调节音高。上图描绘了小提琴的侧视图。Ls代表弦长，即从琴桥到琴枕的距离。红色箭头描绘了一个示例手指位置，演奏大三度音。L”是指手指位置与琴桥的距离。指板穿越系统应能在琴弦上的任意位置停止琴弦，从而实现L”
＜L#。通过改变L”，我们可以改变琴弦长度，进而改变音高。简化起见，假设等间距（实际情况稍有偏差）音阶，音高位置是2 的 12 次方根的函数。对于给定的空弦调音，可以计算出琴桥距离L”在任何品格位置 n 的距离。
$$
L_p=L_s(1-2^{\frac{-n}{12}})
$$

设计需求与理念

设计需求

• 高精度定位：按弦点在指板上的精确位置直接决定音高。因此，设计中采用高精度步进电机或伺服电机驱动，确保按弦位置的准确性，将定位误差控制在0.1mm 以内，以保证音准的准确性和一致性 。

• 快速响应：按弦动作需要快速完成，以适应乐曲节奏的快速变化。通过优化机构传动链和减小运动部件质量，实现按弦动作的快速响应。

• 可调按压力：模拟人类手指按压力度，避免过度按压导致琴弦损伤或按压不足导致音色不佳。通过力传感器反馈和闭环控制，实现按压力的精确调节 。

• 多指模拟：能够同时或快速连续按压多个琴弦，以演奏和弦或快速音阶。设计中实现四弦独立控制，分别对应G 弦（196Hz）、D 弦（294Hz）、A 弦（440Hz）和 E弦（659Hz）的音高 。

设计理念

• 采用齿轮齿条传动（耐磨、结构简单）

• 步进电机驱动（开环控制简单、位置精准、保持转矩）

• 滚珠式直线导轨（精度高、噪音大）

结构设计参数与材料选取

结构参数：

• 按弦杆数量：4 根（试验阶段为 1 根），每根按弦杆独立控制，对应小提琴的四根琴弦 。

• 按弦杆直径：3-5 mm，旨在平衡强度与空间限制，确保按弦杆的刚性 。

• 按压力范围：0.5-2N，可调节的按压力度，模拟人手力度，以保护琴弦并优化音色 。

• 按弦杆行程：5-10 mm，覆盖多把位按压琴弦所需的有效行程 。

• 定位精度：±0.1 mm，确保按弦位置的精确性，保证音高准确 。

• 响应时间：≤10 ms，保证按弦动作的快速性，适应各种演奏技巧 。

材料配置：

• 按弦杆：不锈钢 304/316L，以确保低摩擦和长期稳定性 。

• 导向套：PEEK工程塑料，利用其自润滑特性，减小按弦杆运动时的摩擦阻力 。

• 驱动系统：伺服电机，构成闭环控制系统，实现按弦动作的精确控制 。

• 按压末端：硅藻泥，模拟人手指的触感，保护琴弦以及琴面板，更好还原按压效果

参数设计需求：

• 高精度定位：确保音高准确和节奏流畅，是实现高质量演奏的基础 。

• 可调按压力：保护琴弦，并能够根据音乐表现需求优化音色 。

• 独立控制：实现和弦和复杂音阶的演奏，极大地丰富了机器人的演奏能力。

• 快速响应：使机器人能够适应各种演奏技巧的要求，如快速换把和颤音等。

装配模型

小提琴自动演奏机器人的整体装配模型展示了各个关键构件如何协同工作，以实现复杂的演奏功能。该模型包括指板遍历机构和琴弓控制机构等主要部件，它们共同构成了机器人的核心运动系统。

工作原理流程简述

1. 乐谱解析：机器人首先接收 MIDI 文件，并从中解析出音高、音长、节奏等音乐信息 。

2. 运动规划：根据解析出的音乐信息，系统生成运弓和按弦机构的详细运动轨迹和控制指令 。

3. 运弓驱动：伺服电机驱动运弓进行往复运动，并精确控制弓压和弓速，模拟人类运弓的力度和速度变化。

4. 按弦控制：独立的驱动系统控制各按弦杆，精确按压指板上的琴弦，以产生准确的音高 。

5. 传感器反馈：系统实时监测运弓和按弦机构的运动状态，并通过传感器反馈数据进行闭环控制和校正，确保运动的精确性 。

6. 音色优化：结合人工智能算法，系统实时优化演奏参数，以提升音色表现力，使其更接近人类演奏 。

控制系统简述

功能特点与意义

功能特点：

• 自动化演奏：实现小提琴音乐的完全自动化演奏。

• 高精度音高：通过精确的按弦控制，确保音高准确无误。

• 音色表现力：通过弓压、弓速和角度的精细控制，模拟人类演奏的音色变化。

• 节奏准确性：精确的运动规划和控制确保演奏节奏的准确性。

• 可编程性：用户可以通过编程控制机器人演奏不同的乐曲和风格。

• 实时反馈：传感器反馈系统保证了机器人运动的精确性和稳定性。

技术意义：

• 技术集成验证：装配模型是验证机器人各部件协调性和系统整体可行性的重要基础。

• 直观展示：通过爆炸视图等形式，可以直观地展示机器人内部结构和各部件的工作原理，有助于理解其复杂性。

• 性能评估基础：为后续的仿真分析和性能校核提供了基础模型，有助于在实际制造前发现并解决潜在问题 。

• 创新应用：该项目是机器人技术在艺术领域创新应用的典型范例，展示了机器人技术在非传统领域的巨大潜力。

工艺链选择

关键构件的工艺链选择对于小提琴自动演奏机器人的性能、精度、成本和生产效率具有决定性影响。本节将详细阐述齿轮、运弓构件和按弦机构的完整工艺链，并说明选择这些工艺的理由，强调多种工艺的结合应用。

齿轮工艺链选择

齿轮作为核心传动部件，其制造工艺链必须确保高精度、高强度、优异的耐磨性和低噪音。以下是所用齿轮完整工艺链及其理由：

毛坯制备与预处理

• 锻造制坯：对于高强度齿轮，通常采用热模锻工艺制备齿轮毛坯。锻造能够细化晶粒，改善材料的组织结构，从而显著提高齿轮的抗冲击韧性和疲劳强度。这对于承受频繁启停和变速载荷的齿轮至关重要。

  • 理由：相较于铸造，锻造能够消除铸造缺陷，使金属内部组织更致密，力学性能更优异。虽然铸造适用于大型或结构复杂的齿轮毛坯，但对于小提琴机器人所需的高精度、高强度齿轮，锻造是更优选择 。

• 正火处理：锻造后的齿轮毛坯需要进行正火处理。正火是将钢加热到临界温度以上，保温后在空气中冷却的热处理工艺。其目的是消除锻造过程中产生的内应力，细化晶粒，改善钢的切削性能，为后续的加工和热处理做好组织准备。

  • 理由：正火处理能够均匀化组织，降低硬度，提高塑性，从而便于后续的车削加工，并减少加工硬化现象。

• 车削加工：利用数控车床对正火后的齿轮毛坯进行粗车削。这一步骤旨在高效去除多余材料，形成齿轮的基本几何形状，并提供精确的基准面，为后续的齿形加工奠定基础。

  • 理由：数控车削具有高效率和高精度的特点，能够确保齿轮毛坯的尺寸精度和表面质量，为后续的精密加工提供良好的起点。

齿形加工与热处理

• 滚齿加工：滚齿是制造渐开线齿轮最常用且高效的展成法加工方法。它通过滚刀与齿轮毛坯的相对运动，连续切削出齿形。滚齿加工具有高效率、高精度和加工范围广的优点，适用于大批量生产。

  • 理由：滚齿能够保证齿形精度，对于小提琴机器人齿轮所需的ISO5级精度，滚齿是实现这一目标的关键步骤。同时，其高效性有助于控制生产成本。

• 插齿加工（可选）：对于内齿轮或对齿形精度有更高要求的半精加工，可以采用插齿加工。插齿机通过插齿刀的往复运动和旋转运动，切削出齿形，能够获得更高的齿形精度和表面质量。

  • 理由：在某些特定情况下，如内齿轮或需要更高齿形精度的齿轮，插齿可以作为滚齿的补充或替代工艺，以满足更严苛的设计要求。

• 渗碳淬火回火：这是齿轮制造中至关重要的热处理环节。渗碳是将碳原子渗入齿轮表面，使其表面碳含量增加，然后进行淬火和低温回火。处理后，齿轮表面可获得HRC58-62的高硬度，而心部仍保持高韧性。这显著提高了齿轮的耐磨性和疲劳强度，使其能够承受高载荷和频繁的啮合。

  • 理由：小提琴机器人齿轮需要承受高频次的载荷，渗碳淬火回火能够赋予齿轮“外硬内韧”的特性，确保其在长期运行中的可靠性和寿命。这是提高齿轮性能的关键步骤。

精加工与后处理

• 磨齿精加工：热处理后，齿轮会产生一定的变形。磨齿精加工是消除热处理变形、进一步提高齿形精度（达到ISO5 级甚至更高）的关键步骤。磨齿能够显著降低传动噪音，提高传动的平稳性 。

  • 理由：磨齿是获得高精度齿轮的必要手段，对于小提琴机器人对噪音和运动平稳性的高要求，磨齿能够有效提升齿轮的性能。

• 珩磨/抛光 (可选)：为了进一步改善齿轮表面粗糙度，获得镜面光洁度，可以进行珩磨或抛光处理。这有助于降低齿面摩擦系数，减少噪音，并提高齿轮的传动效率。

  • 理由：在对噪音和传动效率有极致要求的场合，珩磨/抛光能够提供更优异的表面质量，进一步提升齿轮的性能。

• 质量检测：在齿轮制造的各个阶段，特别是最终阶段，需要进行严格的质量检测。检测内容包括齿形误差、齿向误差、表面粗糙度、硬度、尺寸精度等，以确保齿轮满足设计要求和精度等级。

  • 理由： 质量检测是保证产品质量的最后一道防线，确保每个齿轮都符合小提琴自动演奏机器人的严苛要求。

• 表面处理：为了进一步提高齿轮的耐磨性、耐腐蚀性和降低摩擦，可以进行额外的表面处理。例如，对于斜齿轮齿条，可以采用三氯甲烷、氯化烷、冰醋酸的混合溶剂进行微溶解抛光，并施加 PTFE（聚四氟乙烯）润滑喷雾涂层，以降低阻力并减少噪音 。

  • 理由：这些表面处理能够优化齿轮的摩擦学性能，延长使用寿命，并进一步降低运行噪音，对于小提琴机器人这种对噪音敏感的应用至关重要。

齿轮工艺链总结

运弓构件工艺链选择

运弓构件工艺

运弓构件主要采用碳纤维复合材料制造，其工艺链需要确保轻量化、高刚度和优异的力学性能：

• 预浸料铺层与热压罐固化：这是碳纤维复合材料成型的核心工艺。首先，将碳纤维预浸料（预先浸渍树脂的纤维布）按照设计要求进行精确铺层，形成所需的形状。然后，将铺层后的构件放入热压罐中，在高温高压条件下进行固化。这种工艺能够精确控制纤维含量和树脂分布，消除气泡和孔隙，从而获得高性能的复合材料构件。

  • 理由：热压罐固化能够确保复合材料的致密性和力学性能，是制造高强度、高刚度运弓的关键。

• CNC铣削加工：固化成型后的碳纤维运弓需要进行精密加工，以达到最终的尺寸和形状精度。采用高精度CNC铣床进行精加工，并使用专用的金刚石刀具和优化切削参数，以避免碳纤维材料在加工过程中出现分层、毛刺等缺陷。

  • 理由：CNC铣削能够实现复杂形状的精密加工，确保运弓的几何精度和表面质量，满足机器人对运动精度的要求。

• 表面处理：完成机械加工后，运弓可以进行喷漆或涂层处理，以提高其表面耐磨性、耐腐蚀性，并改善外观美观度。

  • 理由：表面处理不仅能保护运弓，还能提升其视觉效果，使其更符合产品设计要求。

运弓构件工艺链总结

按弦机构工艺链选择

按弦机构工艺

按弦机构中的按弦杆主要采用不锈钢制造，其工艺链需要确保高精度、高表面光洁度和优异的耐腐蚀性：

• 精密棒料下料：选用高品质不锈钢精密棒料作为原材料，通过自动送料和精确下料设备，确保原材料的尺寸精度和一致性。

  • 理由：精确的原材料尺寸是后续精密加工的基础，有助于减少材料浪费和加工时间。

• 数控车削与铣削：利用数控车床对不锈钢棒料进行车削，形成按弦杆的圆柱外形、台阶和螺纹等特征。对于键槽等非旋转对称特征，则通过数控铣削完成。数控加工能够确保高精度和高效率。

  • 理由：数控车铣复合加工能够一次装夹完成多道工序，提高加工效率和精度，确保按弦杆的几何尺寸符合设计要求。

• 热处理：根据不锈钢的具体牌号和性能要求，可以进行固溶处理或淬火回火处理。固溶处理旨在消除加工应力，提高材料的塑性和韧性；淬火回火则可以进一步优化材料的力学性能和耐腐蚀性。

  • 理由：热处理能够优化按弦杆的内部组织，提高其强度、硬度和耐腐蚀性，确保其在长期使用中的稳定性和可靠性。

• 精密研磨**/**抛光：为了获得极高的表面光洁度，按弦杆需要进行精密研磨或抛光处理。可以使用无心磨床或专用抛光设备，使按弦杆表面达到镜面效果。

  • 理由：极高的表面光洁度能够显著降低按弦杆与导向套之间的摩擦阻力，提高按弦动作的灵敏度和响应速度，同时也能提升产品的外观质量。

按弦机构工艺链总结：

通过上述详细的工艺链选择和组合，可以确保小提琴自动演奏机器人中各个关键构件的制造质量，从而为机器人的高性能、高精度和长寿命运行提供坚实的基础。

模型仿真及校核

有限元分析(FEA)

有限元分析是一种数值计算方法，用于求解复杂结构在载荷作用下的应力、应变、变形和振动等问题。在小提琴自动演奏机器人中，有限元分析主要用于关键构件的强度、刚度和疲劳寿命评估。

• 应力应变分析：对齿轮、弓臂和按弦杆等关键构件进行应力应变分析，可以识别高应力集中区域，评估构件在最大载荷下的安全裕度，确保其强度满足设计要求，避免发生塑性变形或断裂。

• 变形分析：计算构件在载荷作用下的变形量，验证其刚度是否满足要求。对于高精度机器人，微小的变形都可能导致定位误差，因此严格控制变形量至关重要。

• 疲劳寿命预测：小提琴自动演奏机器人需要长时间、高频率地运行，因此关键部件的疲劳寿命是重要的考量因素。有限元分析可以预测齿轮、轴承等易疲劳部件的疲劳寿命，指导材料选择、结构优化和维护计划，从而延长机器人的使用寿命。

• 仿真软件：常用的有限元分析软件包括 ANSYS、Abaqus 和 SolidWorks Simulation等。这些软件提供了强大的网格划分、材料模型和载荷施加功能，能够进行精确的结构分析，本次仿真使用Fusion 360 内置的 ANSYS插件进行仿真。

仿真结果

以下是齿轮与齿条啮合的仿真分析结果图示例，这些图表直观地展示了齿轮系统在工作状态下的性能表现，是验证齿轮系统性能和安全性的重要依据。

通过上述全面的模型仿真和校核，可以确保小提琴自动演奏机器人的齿轮系统及其他关键构件在设计阶段就达到高标准，为最终产品的成功开发奠定坚实基础。

可持续性分析

小提琴自动演奏机器人项目不仅在技术和艺术领域具有创新意义，其可持续性也是一个值得深入探讨的方面。可持续性分析主要关注项目在环境和经济两个维度上的可行性，以评估其长期发展潜力和社会价值。

环境可行性分析

环境可行性分析旨在评估小提琴自动演奏机器人在其生命周期内对环境可能产生的影响，并探讨如何通过设计和制造优化来降低这些影响。

材料选择与资源消耗

• 积极影响：项目在材料选择上，如齿轮材料推荐合金钢（20CrMnTi）和工程塑料（POM），弓臂采用碳纤维复合材料，按弦机构采用不锈钢。这些材料的选择在保证性能的同时，也考虑了其可回收性和耐用性。例如，合金钢和不锈钢是可回收材料，有助于减少原生资源消耗。工程塑料如POM 在某些应用中可替代金属，降低能耗和材料消耗。
• 潜在挑战与优化：碳纤维复合材料的回收目前仍面临挑战，其生产过程也相对能耗较高。未来可以探索更环保的碳纤维生产技术，或开发易于回收的复合材料。此外，在制造过程中应尽量减少废料产生，并对废料进行分类回收。

能源消耗

• 运行能耗：机器人运行过程中，驱动电机（伺服电机、步进电机、BLDC电机）是主要的能耗来源。通过优化控制算法，提高电机效率，以及采用低功耗的传感器和控制系统，可以有效降低机器人的运行能耗。

• 制造能耗：齿轮的锻造、热处理、机加工以及碳纤维复合材料的固化等工艺都需要消耗大量能源。选择能效更高的制造设备，优化工艺流程，例如采用精密锻造减少后续机加工量，可以降低整体制造能耗。

噪音与振动

• 积极影响：项目设计中特别强调了低噪音运行，例如采用斜齿轮、优化齿形、精密磨齿和表面处理（PTFE 涂层）等措施，旨在将噪音控制在 35dB以下。这不仅提升了演奏音质，也降低了机器人运行对周围环境的噪音污染，尤其是在家庭或教育场所使用时，其环境友好性更高。

• 潜在挑战与优化：尽管采取了降噪措施，但机械运动仍会产生一定噪音。未来可以进一步研究主动降噪技术，或采用更先进的静音材料和结构设计。

废弃物管理

• 设计优化：在设计阶段就考虑产品的可拆卸性和模块化，便于报废时对不同材料进行分类回收。例如，弓毛部分的模块化设计便于更换和回收。

• 生命周期评估(LCA)：进行全面的生命周期评估，从原材料获取、制造、使用到报废回收的全过程进行环境影响分析，识别环境热点，并制定相应的改进策略。

经济可行性分析

经济可行性分析旨在评估小提琴自动演奏机器人项目的经济效益、市场潜力、成本控制和商业模式，以确保其可持续发展。

市场潜力与需求

• 艺术与教育市场：小提琴自动演奏机器人面向艺术表演、音乐教育和个人娱乐等多个市场。在音乐教育领域，它可以作为辅助教学工具，降低学习门槛，激发学习兴趣。在艺术表演领域，可以拓展新的表演形式和创作可能性。随着人们对高品质音乐体验和个性化娱乐需求的增长，市场潜力巨大。

• 技术发展趋势：机器人技术、人工智能和精密制造的不断进步，将进一步降低生产成本，提高产品性能，从而扩大市场接受度。

成本控制与效益

• 研发成本：项目涉及多学科交叉，研发投入相对较高。然而，通过模块化设计、标准化部件和高效的仿真校核，可以有效缩短研发周期，降低试错成本。

• 制造成本：齿轮的精密加工、碳纤维复合材料的成型等工艺成本较高。但随着制造技术的成熟和规模化生产，单位成本有望降低。例如，3D打印技术在小批量生产或原型制造中具有成本优势，未来可探索其在某些部件上的应用以降低成本。

• 运营与维护成本：高耐用性材料和精密设计有助于降低机器人的故障率和维护成本。模块化设计也便于部件更换和维修，进一步降低了长期运营成本。

商业模式与盈利

• 产品销售：直接销售小提琴自动演奏机器人给个人用户、音乐学院、演出团体等。

• 服务模式：提供机器人租赁、定制化演奏服务、音乐教育课程等。

• 技术授权：将核心技术授权给其他制造商，扩大技术影响力并获取收益。

• 内容生态：开发配套的乐谱库、演奏技巧库和 AI优化算法，形成内容生态，增加用户粘性。

综上所述，小提琴自动演奏机器人项目在环境和经济方面都具备良好的可行性。通过在设计、制造和运营全生命周期中融入可持续发展理念，并积极探索多元化的商业模式，该项目有望实现长期的社会、环境和经济效益。

总结

小提琴自动演奏机器人项目是精密机械、电子、控制和人工智能等多学科深度融合的典范，旨在模拟人类小提琴演奏家的精湛技艺，实现音乐的自动化演奏。本报告围绕该项目，详细阐述了其核心传动部件——齿轮系统的设计与制造，并深入探讨了相关材料选择、关键构件设计、工艺链选择、模型仿真及校核，以及项目的可持续性。

项目背景与意义部分强调了机器人技术在艺术领域的独特潜力，以及小提琴演奏对人类生理和心理素质的极高要求。该项目不仅推动了精密机械、控制系统和人工智能等前沿技术的发展，更在创新音乐教育、拓展艺术表现形式、传承文化遗产以及培养复合型人才方面具有深远意义。其发展趋势预示着更高精度、智能化、多模态融合和模块化定制的未来。

在材料选择方面，报告针对齿轮、弓臂构件和按弦机构等关键部件，推荐了合金钢（如20CrMnTi）、碳纤维复合材料和不锈钢（如304/316L）等高性能材料，并详细阐述了选择理由。这些材料的选择旨在满足机器人对高强度、高刚度、轻量化、低噪音和耐腐蚀等严苛性能要求，并通过材料对比分析，展现了在性能与成本之间进行权衡的综合考量。

关键构件设计与建模部分，以齿轮为例，详细介绍了其基本参数计算、设计理念（高精度传动、低噪音运行、紧凑型设计、长寿命可靠性）和斜齿轮齿条的优化目标。同时，弓臂构件和按弦机构的设计也充分考虑了轻量化、多自由度控制、高精度定位和快速响应等特性。装配模型部分则清晰地展示了机器人的工作原理流程、功能特点和技术意义，凸显了各部件的协同作用。

工艺链选择是确保机器人性能的关键环节。报告详细阐述了齿轮的完整工艺链，包括锻造制坯、正火、车削、滚齿、渗碳淬火回火、磨齿精加工以及珩磨/抛光和表面处理等，并强调了每一步工艺对齿轮性能提升的重要性。弓臂构件的碳纤维复合材料制造工艺（预浸料铺层与热压罐固化、CNC铣削、表面处理）和按弦机构的精密制造工艺（精密棒料下料、数控车铣、热处理、精密研磨/抛光）也得到了详细介绍，展现了多工艺结合以实现高性能部件制造的策略。

模型仿真及校核部分，强调了运动学仿真、动力学仿真和有限元分析在设计验证、安全评估和优化指导中的重要作用。通过这些仿真手段，可以在物理样机制造前发现并解决潜在问题，确保齿轮系统及其他关键构件在设计阶段就达到高标准，从而提高设计的可靠性和效率。

最后，可持续性分析从环境和经济两个维度评估了项目的可行性。在环境方面，通过材料选择、能源消耗优化和噪音控制等措施，努力降低对环境的影响。在经济方面，项目具有广阔的市场潜力，并通过成本控制、多元化商业模式和人才培养，确保其长期可持续发展。

总而言之，小提琴自动演奏机器人项目不仅是工程技术与艺术的完美结合，更是一个系统性、多学科交叉的复杂工程实践。通过对齿轮系统及其他关键构件的精细设计、材料优选、先进制造工艺的应用以及全面的仿真校核，该项目为实现高精度、高性能的自动化音乐演奏奠定了坚实基础，展现了未来机器人技术在非传统领域的巨大应用前景和价值。

参考文献

1. 黄信行, 黄淑芳, 李文鸿, 邱燕玉, 游凯程, 温扬圣. 小提琴机器人－六轴机器手臂的整合运用. 机械工业网 , 2012(4): 124-135.[https://www.automan.tw/tw/magazine-period/54/982]

2. 抗疲劳齿轮的研究进展：提升机械可靠性的新方法. 温岭市螺伞齿轮厂.[http://www.lsgears.com/NewsDetail.aspx?ID=78]

3. 齿 轮 制 造 101 ： 齿 轮 生 产 过 程 指 南 . RapidDirect.[https://www.rapiddirect.com/zh-CN/blog/gear-manufacturing/]

4. 某车型汽车差速器结构设计及其齿轮的性能研究. hanspub.org, 2024.[https://pdf.hanspub.org/mos2024135_102571870.pdf]

5. 协作机器人结构设计及齿轮传动系统仿真. hanspub.org, 2024.[https://www.hanspub.org/journal/paperinformation?paperid=79757]

6. 粉末锻造齿轮材料的组织与性能研究. 粉末冶金技术, 2020.[https://pmt.ustb.edu.cn/cn/article/id/fmyjjs202002006]

7. 齿 轮 加 工 工 艺 流 程 分 析 . ResearchGate, 2025.[https://www.researchgate.net/publication/384875475_chilunjiagonggongyiliuchengfenxi]

8. 运用力量控制技术于提升小提琴机器人演奏音色之探讨. 华艺线上图书馆.[https://www.airitilibrary.com/Article/Detail/U0017-2211202315312089]

9. 工业机器人领域谐波齿轮传动应用及运动精度分析. 中国学术期刊网络出版总库.[https://cpfd.cnki.com.cn/Article/CPFDTOTAL-KDIE201710001126.htm]

10. 齿轮用什么材料，齿轮常用材料有哪几种. 知乎专栏, 2023.[https://zhuanlan.zhihu.com/p/666265403]

1 研究背景简介

1.1 研究背景与意义

香烟烟雾中含有数千种化学物质，其中包括一氧化碳、尼古丁、苯等多种对人体有害的成分。这些物质在室内封闭环境中容易积聚，对人体健康构成显著威胁，尤其是二手烟和三手烟的危害已引起广泛关注。报告显示，2018年中国不吸烟者的二手烟暴露率为68.1%，其中家庭和工作场所是二手烟暴露的主要场所，这表明室内烟雾污染已成为影响公众健康的重要问题。因此，我们希望通过研究烟雾在室内环境中的扩散规律，为烟雾的控制、空气污染防治或室内空气质量管理提供参考，并帮助优化通风设计或制定健康防护措施。

烟雾的扩散过程受到室内温湿度、气流模式、房间结构及家具布置等多种因素的影响，具有较高的复杂性和随机性。深入研究烟雾扩散的传输机制，可以帮助揭示其污染范围与浓度分布，为优化室内空气质量管理提供理论支持。但在估计和预测吸入暴露风险时，由于与生物伦理学和动物保护相关的潜在限制，涉及人类志愿者和其他用于毒理学研究的哺乳动物替代模型的体内研究存在局限性。出于这样的考虑，我们希望通过理论分析建模与数值模拟的方式，得到烟雾扩散的一般规律，从而为提出可能的防控手段提供理论依据。

通过实验和仿真研究，能够验证不同通风方案对烟雾扩散的控制效果，为改善居住环境提供技术支持；利用数值模拟，可以精确再现不同环境下烟雾的扩散过程，从而更直观地评估各类干预措施的有效性。这种方法不仅降低了研究成本，还能为制定更加科学、精确的室内空气污染防控策略提供强有力的支持。

1.2 研究目标

本研究旨在探讨室内封闭环境下烟雾的扩散特性，基于扩散方程等物理原理与数学物理方法对烟雾的扩散情况进行解析求解与数值模拟，重点分析烟雾的流动路径、扩散速率和浓度分布，以优化空气管理和烟雾控制策略。

在具体的研究过程中，我们主要采用特定尺寸的密封容器来模拟室内封闭环境以便于实验开展与现象观察；同时，我们的研究主要关注扩散过程中烟雾的扩散效果、扩散速度以及不同时刻的浓度分布这三个要素，并通过定量计算测量与定性观察现象相结合的方式，得到烟雾扩散浓度空间分布随时间变化的的一般数学物理规律，从而为设计室内烟雾浓度控制系统提供理论依据。

2 理论分析求解

2.1 研究对象

在实际研究过程中，为方便实验的开展，我们选择了粒径分布在0.1~10微米范围内的植物甘油（丙三醇）气雾作为研究对象，这也是电子烟烟雾的主要成分之一。之所以使用植物甘油气雾来代替香烟烟雾进行研究，主要出于以下考量：

• 安全性：植物甘油为食品级材料，无健康与环境污染风险

• 可控性：植物甘油气雾生成过程可控，可进行定量化分析

• 易获取性：价格低廉，易于获取

• 可重复性：物质纯度高，化学性质稳定

具体而言，植物甘油气雾具有以下几条显著的理化特性：

• 无色无臭，有甜味，极易溶于水

• 密度（25℃）：1.26 g/cm^3

• 分子量： 92.09 g/mol

• 主要粒径属于典型气溶胶粒径范围

• 挥发性极低，以小颗粒悬浮于空气中

在研究过程中进行这样的材料替换，具有如下的合理性：

• 植物甘油气雾与香烟烟雾的气溶胶形成机制相同（蒸发-冷凝），二者在扩散和湍流行为上具有共性

• 二者粒径分布相似，具备相似的扩散动力学特征

• 二者在扩散过程中的视觉效果接近，便于追踪扩散行为

在本研究的理论分析部分，重点关注植物甘油气雾的浓度（即单位体积中该种物质颗粒的质量）ρ(x, y, z, t)这一物理量在研究区域（长为x_0 = 40cm，宽为y_0 = 50cm，高为z_0 = 60cm的长方体规则区域）空间范围内的分布情况及其随时间的变化情况。

2.2 数学物理方程的导出

针对烟雾扩散这一常见现象，有现成的实验规律——菲克定律（又称扩散定律）可以直接使用，其基本公式如下：

$$
\overrightarrow{q} = - D\nabla\rho
$$
其中：

• 式中负号代表扩散转移的方向（浓度减小的方向）和浓度梯度（浓度增大的方向）相反；
• D：扩散系数，与物质的种类以及环境温度有关；
• ∇ρ：浓度梯度，描述浓度不均匀分布的程度；
• 扩散流强度，即单位时间通过单位（横截）面积的质量，描述扩散运动的强度：

$$
\overrightarrow{q} = \frac{dm_{g}}{dVdt}
$$
在本研究中，针对常温环境（25℃）下植物甘油气雾在空气中的扩散系数，可进一步通过Chapman-Enskog方程
$$
D_{AB} = \frac{0.00143 \times T^{1.75} \times \sqrt{\frac{1}{M_{A}} + \frac{1}{M_{B}}}}{P \times \sigma_{AB}^{2} \times \Omega_{D}}
$$
进行估算，代入空气的平均分子量28.97 g/mol和估计的分子直径0.5nm可计算得到此时的扩散系数约为：
$$
D = 0.0962{cm}^{2}/s
$$
由于我们的研究区域是一个长方体的三维规则区域，故在研究坐标系的选择上，以长方体区域的底端顶点为原点，长、宽、高延伸方向作为$x、y、z$正方向建立空间直角坐标系，原点在长方体上的对角顶点坐标为$(x_{0},\ y_{0},z_{0})$。针对三个不同方向，还可以进一步地给出菲克定律的分量形式：

$$
q_{x} = - D\frac{\partial\rho}{\partial x},\ \ q_{y} = - D\frac{\partial\rho}{\partial y},\ \ q_{z} = - D\frac{\partial\rho}{\partial z}
$$
接下来将基于菲克定律和质量守恒定律给出三维扩散方程，即本研究涉及的数学物理方程：

对于空间中任一点(x, y, z)附近的无穷小空间而言：该空间内浓度变化取决于穿过其表面的扩散流。

考虑单位时间内x方向的扩散流：

左表面流入流量
$$
\left. \ q_{x} \right|{x}dydz
$$
右表面流出流量
$$
\left. \ q{x} \right|{x + dx}dydz
$$
故在单位时间内x方向净流入流量为
$$
\left. \ q{x} \right|{x}dydz - \left. \ q{x} \right|_{x + dx}dydz
$$

由于取无穷小空间，dx足够小，于是有
$$
\frac{\partial q_{x}}{\partial x} = \frac{\left. \ q_{x} \right|{x + dx} - \left. \ q{x} \right|{x}}{dx} 即 \left. \ q{x} \right|{x + dx} - \left. \ q{x} \right|{x} = \frac{\partial q{x}}{\partial x}dx
$$
代入化简后可以得到单位时间内x方向净流入流量为
$$
D\frac{\partial^{2}\rho}{\partial x^{2}}dxdydz
$$
同理可得：单位时间内y方向净流入流量为
$$
D\frac{\partial^{2}\rho}{\partial y^{2}}dxdydz
$$
单位时间内z方向净流入流量为
$$
D\frac{\partial^{2}\rho}{\partial z^{2}}dxdydz
$$
在本研究中，出于解析求解可行性的考虑，在理论分析部分暂时仅考虑常压下有源无汇的理想场景（即仅在某固定坐标点(x_0, y_0, z_0)处有一恒定扩散源，其强度（单位时间内单位体积产生的气体质量）为F(x,y, z,t)），而在后续数值求解与仿真模拟时才会把负压源视作汇的产生原因从而进一步补充对流的相关条件（对流项）。

在此条件下，基于质量守恒定律，单位时间在该无穷小空间内增加的气体质量等于单位时间内净流入的气体质量与扩散源产生的气体质量，其中扩散源产生的气体质量F(x,y, z,t)仅在(x_0, y_0, z_0)处有恒定值F_0而在其余空间各处均为0。于是列出并化简得到如下三维扩散方程（相比通用控制方程暂时不考虑对流项）：

$$
\frac{\partial\rho}{\partial t} - D\left( \frac{\partial^{2}\rho}{\partial x^{2}} + \frac{\partial^{2}\rho}{\partial y^{2}} + \frac{\partial^{2}\rho}{\partial z^{2}} \right) = \rho_{t} - D\nabla^{2}\rho = 0, 0 < x < x_{0}, 0 < y < y_{0}, 0 < z < z_{0}
$$
其中：

• 浓度的时间增长率（关于时间t的偏导数），即单位时间该无穷小空间内增加的气体质量:
  $$
  \frac{\partial\rho}{\partial t}
  $$

• 扩散系数D仅与物质种类和温度有关，若研究某种特定气体在理想的均匀介质（温度在空间各处保持一致且不随时间变化）条件下的扩散，该系数值为常数；反之则其也应表示为关于空间坐标与时间的函数D(x, y, z, t)（在确定气体类型时仅与温度的分布函数T(x, y, z, t)有关）。

2.3 初始条件与边界条件

该问题的初始条件是显然的：研究区域空间内在初始状态（t = 0）没有该种烟雾气体分布，且在某固定点处有一恒定扩散源，始终以恒定的扩散强度（单位时间内单位体积产生的气体质量）产生气体并向区域空间内扩散。将这样的初始条件写成数学表达式的形式：

$$
\rho(x,y,\ z,0) = F_{0}\delta(x - \frac{x_{0}}{2})\delta(y - \ \frac{y_{0}}{2})\delta(z - 0)
$$
式中F_0表示在固定坐标端点
$$
\left( \frac{x_{0}}{2},\ \frac{y_{0}}{2},\ 0 \right)
$$
处的恒定扩散源的扩散强度。结合实验时采用的实验装置，参考其技术文档可知，其值约为5g/cm^3。

关于边界条件，针对我们所研究的扩散问题，通常采用Dirichlet边界条件：
$$
\frac{\partial\rho}{\partial n} = 0
$$
显然这是一个齐次边界条件，也是一个典型的第二类边界条件。为方便后续求解处理，针对我们所建立的空间直角坐标系，可以将其进一步写成分量形式：

$$
\left{ \begin{array}{r}
\left. \ \frac{\partial\rho}{\partial x} \right|{x = 0} = \left. \ \frac{\partial\rho}{\partial x} \right|{x = x_{0}} = 0 \
\left. \ \frac{\partial\rho}{\partial y} \right|{y = 0} = \left. \ \frac{\partial\rho}{\partial y} \right|{y = y_{0}} = 0 \
\left. \ \frac{\partial\rho}{\partial z} \right|{z = 0} = \left. \ \frac{\partial\rho}{\partial z} \right|{z = z_{0}} = 0
\end{array} \right.
$$

2.4 数学物理方程求解

针对2.2节提出的数学物理方程，可以通过分离变数法，结合2.3节给出的边界条件与初始条件，得到本征值问题的解和本征值的取值，并将多个本征值问题的解依次相乘得到本征解，再对所有可能的本征值求和，即可得到级数解。更进一步的，通过2.3节给出的初始条件，还可以对于级数解中的常数系数进行进一步确定，从而得到该数学物理方程完整的级数解解析表达式。

2.4.1 变量分离

对于扩散的烟雾而言，虽然其浓度会在空间中的不同位置随时间而变化，但时间与空间这两者在浓度瞬态变化的过程中实际上可以看作是相互独立的，即浓度随时间变化的速率ρ_t与浓度在空间分布中的变化速率∇ρ无关。

基于这样的假设，可以对烟雾浓度函数ρ(x, y, z, t)做如下的变量分离：

$$
\rho(x, y, z, t) = v(x,y,z) \bullet T(t)
$$

• 空间部分v(x,y,z)：描述烟雾浓度在三维空间中的分布，可以认为它与时间无关，或者说在某个时刻，密度的空间分布是静态的；

• 时间部分T(t)：描述烟雾密度随时间的变化，通常假设它是一个仅依赖于时间的函数，反映了扩散过程中气体密度随时间的衰减或增长。

在该式中，自变数x, y, z只出现于v之中，自变数t只出现于T之中，烟雾浓度的一般表示式具有分离变数的形式。

更进一步的，由于在空间中x、y、 z方向上的烟雾浓度变化情况也彼此独立，还可以对于空间部分v(x,y,z)做进一步的变量分离，即：

$$
v(x,y,z) = X\left( x) \bullet Y(y) \bullet Z(z \right)
$$
其中，X(x)为烟雾浓度在空间中x方向的分布情况，Y(y)为烟雾浓度在空间中y方向的分布情况，Z(z)为烟雾浓度在空间中z方向的分布情况。

2.4.2 常微分方程与本征值问题

显然三维扩散方程
$$
\rho_{t} - D\nabla^{2}\rho = 0
$$
为齐次方程，对于该齐次方程，可以将分离变量后的烟雾浓度代入三维扩散方程中，化简得到如下式子：

$$
T_{t}(t)v(x,y,z) = DT(t)\nabla^{2}v(x,y,z)，即\frac{T_{t}(t)}{DT(t)} = \frac{\nabla^{2}v(x,y,z)}{v(x,y,z)} = \frac{X^{‘’}(x)}{X(x)} = \frac{Y^{‘’}(y)}{Y(y)} = \frac{Z^{‘’}(z)}{Z(z)}
$$
其中，左边均为关于时间t的函数，与坐标位置x, y, z无关；右边则是关于坐标位置x, y, z的函数，与时间t无关。两边相等是显然不可能的，除非两边实际上是同一个常数，不妨将其记为-λ，于是有：

$$
\frac{T_{t}}{DT} = \frac{X^{‘’}}{X} = \frac{Y^{‘’}}{Y} = \frac{Z^{‘’}}{Z} = - \lambda
$$
将该式关于空间部分X(x)、Y(y)、Z(z)与时间部分T(t)分别分离，可以得到关于X、Y、Z以及关于T的一系列常微分方程：

$$
\left{ \begin{array}{r}
T^{‘} + \lambda DT = 0 \
\nabla^{2}v + \lambda v = 0
\end{array} \right.\ 即 \left{ \begin{array}{r}
T^{‘} + \lambda DT = 0 \
X^{‘’} + \lambda X = 0 \
Y^{‘’} + \lambda Y = 0 \
Z^{‘’} + \lambda Z = 0
\end{array} \right.
$$
考虑2.3节给出的齐次边界条件即Dirichlet边界条件：边界处
$$
\frac{\partial\rho}{\partial n} = 0
$$
有：

$$
\left{ \begin{array}{r}
X^{‘}(0) = X^{‘}\left( x_{0} \right) = 0 \
Y^{‘}(0) = Y^{‘}\left( y_{0} \right) = 0 \
Z^{‘}(0) = Z^{‘}\left( z_{0} \right) = 0
\end{array} \right.\
$$
以上常微分方程与齐次边界条件共同构成了本征值问题。

2.4.3 本征值问题求解

首先，对于第一个方程
$$
\mathbf{T}^{‘} + \mathbf{\lambda DT} = \mathbf{0}
$$
而言，该方程为一阶常微分方程，可直接进行求解，结合初始条件
$$
\rho(x,y,\ z,0) = F_{0}\delta(x - \frac{x_{0}}{2})\delta(y - \ \frac{y_{0}}{2})\delta(z - 0)
$$
可得到：

$$
\mathbf{T}\left( \mathbf{t} \right) = \mathbf{F}_{\mathbf{0}}(\mathbf{1} - \mathbb{e}^{- \mathbf{\lambda Dt}})（C为某常数）
$$
其物理意义上的合理性如下：

• 初始时刻t = 0：
  $$
  T(0) = F_{0}\left( 1 - \mathbb{e}^{- 0} \right) = 0
  $$
  即空间内在初始状态下（t = 0）没有该种烟雾气体分布，与初始条件吻合；

• 稳定时刻t->∞：
  $$
  T(\infty) = F_{0}\left( 1 - \mathbb{e}^{- \infty} \right) = F_{0}
  $$
  即浓度达到稳定值F_0，与扩散物理过程实际现象一致；

• 时间项指数因子-λDt决定了扩散速度，其中D为扩散系数，越大扩散越快，λ为空间模态本征值，模态越高衰减越快，这也符合实际物理扩散过程。

而针对后面三个方程，与对应的齐次边界条件形成了三组本征值问题，由于其形式完全一致，下面将以其中一组本征值问题
$$
\left{ \begin{array}{r}
X^{‘’} + \lambda X = 0 \
X^{‘}(0) = X^{‘}\left( x_{0} \right) = 0
\end{array} \right.\
$$
为例，利用幂级数法求得其本征解：

选定点x_1 = 0，显然常微分方程
$$
X^{‘’} + \lambda X = 0
$$
的系数函数在该点领域内是解析的，因此该点为方程
$$
X^{‘’} + \lambda X = 0
$$
的常点。针对常点x_1，可以给出X(x)在该点邻域上的泰勒级数形式：

$$
X(x) = \sum_{k = 0}^{\infty}{a_{k}{(x - x_{1})}^{k}}
$$
将其代入常微分方程可得：
$$
\sum_{k = 0}^{\infty}{\lbrack(k + 2)(k + 1)a_{k + 2} + \lambda}a_{k}\rbrack{(x - x_{1})}^{k} = 0
$$
于是有：
$$
(k + 2)(k + 1)a_{k + 2} + {\lambda a}{k} = 0,\ k = 0,1,2\ldots
$$
可进一步得到系数递推公式：
$$
a{k + 2} = \frac{- \lambda}{(k + 2)(k + 1)}a_{k}
$$
推得：
$$
\left{ \begin{array}{r}
a_{2k} = \frac{\lambda^{k}{( - 1)}^{k}}{(2k)!}a_{0} \
a_{2k + 1} = \frac{\lambda^{k}{( - 1)}^{k}}{(2k + 1)!}a_{1}
\end{array} \right.\ ,\ k = 1,2,3\ldots
$$
因此可得到常微分方程
$$
X^{‘’} + \lambda X = 0
$$
的解：

$$
X(x) = a_{0}X_{0}(x) + a_{1}X_{1}(x)
$$
其中：

$$
X_{0}(x) = 1 + \frac{- \lambda}{2!}x^{2} + \frac{\lambda^{2}}{4!}x^{4} + \ldots + \frac{\lambda^{k}( - 1)^{k}}{(2k)!}x^{2k} + \ldots = \sum_{k = 0}^{\infty}{\frac{\lambda^{k}( - 1)^{k}}{(2k)!}x^{2k}} = \cos(\sqrt{\lambda}x)
$$

$$
X_{1}(x) = x + \frac{- \lambda}{3!}x^{3} + \frac{\lambda^{2}}{5!}x^{5} + \ldots + \frac{\lambda^{k}( - 1)^{k}}{(2k + 1)!}x^{2k + 1} + \ldots = \sum_{k = 0}^{\infty}{\frac{\lambda^{k}( - 1)^{k}}{(2k + 1)!}x^{2k + 1}} = \frac{\sin(\sqrt{\lambda}x)}{\sqrt{\lambda}}
$$

此时的收敛半径为
$$
R = \lim_{n \rightarrow \infty}\left| \frac{a_{n}}{a_{n + 2}} \right| = \lim_{n \rightarrow \infty}\left| \frac{(n + 2)(n + 1)}{- \lambda} \right| = \infty
$$
即级数解X_0(x)$和X_1(x)在复数域上始终收敛。X_0(x)仅含x的偶次幂，为偶函数；X_1(x)仅含x的奇次幂，为奇函数。同时，式中含有sqrt（λ）项，这要求λ≥0。

于是有：
$$
X^{‘}(x) = - a_{0}\sqrt{\lambda}\sin(\sqrt{\lambda}x) + a_{1}\cos(\sqrt{\lambda}x)
$$
将齐次边界条件
$$
X^{‘}(0) = X^{‘}\left( x_{0} \right) = 0
$$
代入可得:

$$
\left{ \begin{array}{r}
a_{1} = 0 \a_{0}\sqrt{\lambda}\sin(\sqrt{\lambda}x_{0}) + a_{1}\cos(\sqrt{\lambda}x_{0}) = 0
\end{array} \right.\
$$
显然有a_1= 0，但对于a_0而言，如果a_0也为0，这样的解是没有意义的；而要使得a_0不为零，就要求
$$
\lambda = 0或\sin\left( \sqrt{\lambda}x_{0} \right) = 0
$$
恒成立，这意味着
$$
\sqrt{\lambda}x_{0} = k\pi
$$
由于
$$
x_{0} > 0、\sqrt{\lambda} \geq 0
$$
于是k的可取值为
$$
k = 0,1,2,\ldots
$$
在这样的条件下，由于x_0与π均为定值，本征值λ则有常数k唯一确定，此时不妨将常数k视为新的本征值。于是可得到该本征值问题的本征解：

$$
X(x) = a_{0}\cos\left( \sqrt{\lambda}x \right) = a_{0x}\cos\left( \frac{k\pi}{x_{0}}x \right)，k = 0,1,2,\ldots
$$
与此类似的，另外两个本征值问题的本征解为：

$$
Y(y) = a_{0y}\cos\left( \frac{m\pi}{y_{0}}y \right)，m = 0,1,2,\ldots
$$

$$
Z(z) = a_{0z}\cos\left( \frac{n\pi}{z_{0}}z \right)，n = 0,1,2,\ldots
$$

值得注意的是，三个本征值问题中出现了三个新的本征值k、m、n，他们分别与时间部分涉及到的本征值λ有定量关系，但是当他们取不同值时，原本征值λ也会有不同的取值，即本质上有三个不同的本征值k、m、n取代了原先的本征值λ。

基于以上本征值问题的解和本征值的取值，可以将其各部分相乘得到本征解，再对本征解按照所有可能的本征值求和，即可得到烟雾浓度的级数解：

$$
\rho(x,\ y,\ z,\ t) = \sum_{k = 0}^{\infty}{\sum_{m = 0}^{\infty}{\sum_{n = 0}^{\infty}{a_{0x}a_{0y}a_{0z}F_{0}\cos\left( \frac{k\pi}{x_{0}}x \right)\cos\left( \frac{m\pi}{y_{0}}y \right)\cos\left( \frac{n\pi}{z_{0}}z \right)}}}(1 - \mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t})
$$
由于式中a_0x、a_0y、a_0z、F_0均为常数（可由初始条件
$$
\rho(x,y,\ z,0) = F_{0}\delta(x - \frac{x_{0}}{2})\delta(y - \ \frac{y_{0}}{2})\delta(z - 0)
$$
进一步确定），不妨将其替换为某一常数
$$
A = a_{0x}a_{0y}a_{0z}F_{0}
$$
以便于后续计算，即此时可得到级数解为：
$$
\rho(x,\ y,\ z,\ t) = \sum_{k = 0}^{\infty}{\sum_{m = 0}^{\infty}{\sum_{n = 0}^{\infty}{A\cos\left( \frac{k\pi}{x_{0}}x \right)\cos\left( \frac{m\pi}{y_{0}}y \right)\cos\left( \frac{n\pi}{z_{0}}z \right)}}}(1 - \mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t})
$$

2.4.4 级数解常数系数的确定

需要指出，函数
$$
\varnothing_{kmn} = \cos\left( \frac{k\pi}{x_{0}}x \right)\cos\left( \frac{m\pi}{y_{0}}y \right)\cos\left( \frac{n\pi}{z_{0}}z \right)
$$
是定义在长方体区域
$$
x \in \left\lbrack 0,x_{0} \right\rbrack,\ \ y \in \left\lbrack 0,y_{0} \right\rbrack,z \in \lbrack 0,z_{0}\rbrack
$$
上的一组正交基函数，其满足如下正交性与归一化条件：

$$
\int_{0}^{z_{0}}{\int_{0}^{y_{0}}{\int_{0}^{x_{0}}{\varnothing_{kmn}\varnothing_{k^{‘}m^{‘}n^{‘}}dxdydz}}} = \left{ \begin{array}{r}
x_{0}y_{0}z_{0},\ (k,m,n) = (k^{‘},m^{‘},n^{‘}) \
0,(k,m,n) \neq (k^{‘},m^{‘},n^{‘})
\end{array} \right.\
$$
对于级数解
$$
\rho(x,\ y,\ z,\ t) = \sum_{k = 0}^{\infty}{\sum_{m = 0}^{\infty}{\sum_{n = 0}^{\infty}{A\cos\left( \frac{k\pi}{x_{0}}x \right)\cos\left( \frac{m\pi}{y_{0}}y \right)\cos\left( \frac{n\pi}{z_{0}}z \right)}}}\mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t} = \sum_{k = 0}^{\infty}{\sum_{m = 0}^{\infty}{\sum_{n = 0}^{\infty}{A\varnothing_{kmn}}}}{(1 - \mathbb{e}}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t})
$$
而言，将两侧同时与Φ_k’m’n’相乘，并在区域内积分：

$$
\int_{0}^{z_{0}}{\int_{0}^{y_{0}}{\int_{0}^{x_{0}}{\rho(x,\ y,\ z,\ t)\varnothing_{k^{‘}m^{‘}n^{‘}}dxdydz}}} = \int_{0}^{z_{0}}{\int_{0}^{y_{0}}{\int_{0}^{x_{0}}{(\sum_{k = 0}^{\infty}{\sum_{m = 0}^{\infty}{\sum_{n = 0}^{\infty}{A\varnothing_{kmn}}}}\mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t})\varnothing_{k^{‘}m^{‘}n^{‘}}dxdydz}}}
$$
利用正交性条件，上述式子右边仅保留(k,m,n) = (k’, m’, n’)的项，其余项积分为0：

$$
\int_{0}^{z_{0}}{\int_{0}^{y_{0}}{\int_{0}^{x_{0}}{\rho(x,\ y,\ z,\ t)\varnothing_{kmn}dxdydz}}} = A(1 - \mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t})x_{0}y_{0}z_{0}
$$
代入初始条件
$$
\rho(x,y,\ z,0) = F_{0}\delta(x - \frac{x_{0}}{2})\delta(y - \ \frac{y_{0}}{2})\delta(z - 0)
$$
即可在t = 0处解出系数
$$
A = \frac{\int_{0}^{z_{0}}{\int_{0}^{y_{0}}{\int_{0}^{x_{0}}{\rho(x,\ y,\ z,\ t)\varnothing_{kmn}dxdydz}}}}{\left( 1 - \mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t} \right)x_{0}y_{0}z_{0}}
$$
由于分子分母在t = 0处均为0，故利用洛必达法则上下对时间$t$求导得：

$$
\rho_{t} = T^{‘} = \mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)Dt}（忽略系数）
$$

$$
A = \frac{\left. \rho_{t} \right|{t = 0}\int{0}^{z_{0}}{\int_{0}^{y_{0}}{\int_{0}^{x_{0}}{\rho(x, y, z, 0)\varnothing_{kmn}dxdydz}}}}{\left( {\left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}\mathbb{e}}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t} \right)x_{0}y_{0}z_{0}}
= \frac{F_{0}}{D\left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)\pi^{2}x_{0}y_{0}z_{0}}\varnothing_{kmn}\left( \frac{x_{0}}{2},\frac{y_{0}}{2},0 \right)
= \frac{F_{0}}{D\left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)\pi^{2}x_{0}y_{0}z_{0}}\cos\left( \frac{k\pi}{x_{0}}\frac{x_{0}}{2} \right)\cos\left( \frac{m\pi}{y_{0}}\frac{y_{0}}{2} \right)\cos\left( \frac{n\pi}{z_{0}}0 \right) = \frac{F_{0}}{D\left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)\pi^{2}x_{0}y_{0}z_{0}}\cos\left( \frac{k\pi}{2} \right)\cos\left( \frac{m\pi}{2} \right)，其值随k,m,n的变化而变化
$$

最后，将其代入2.4.3节解得的级数解，可以最终得到植物甘油气雾浓度的解析表达式为：

$$
\rho(x,\ y,\ z,\ t) = \frac{F_{0}}{Dx_{0}y_{0}z_{0}}\sum_{k = 0}^{\infty}{\sum_{m = 0}^{\infty}{\sum_{n = 0}^{\infty}{\frac{1}{\left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)\pi^{2}}\cos\left( \frac{k\pi}{2} \right)\cos\left( \frac{m\pi}{2} \right)\cos\left( \frac{k\pi}{x_{0}}x \right)\cos\left( \frac{m\pi}{y_{0}}y \right)\cos\left( \frac{n\pi}{z_{0}}z \right)}}}(1 - \mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t})
$$

2.5 基于MATLAB的级数解可视化

基于上述求解过程得到的烟雾浓度解析表达式，编写了对应的MATLAB程序，将不同时刻烟雾浓度在空间内各点处的分布情况进行计算与可视化。值得注意的是，由于级数解本质是无限项求和，而在MATLAB程序实际计算中需要指定求和项数，在此综合考虑计算精度与算力需求，选取$k,m,n$三个本征值各10项进行计算与可视化；除此之外，为与实际实验相对应，计算的扩散时间设置为30秒，同时为直观观察扩散过程中烟雾浓度随时间的变化，设置计算步长为1秒。

代入相应的数值并运行程序可以得到如下的可视化结果（代码详见附件1）：

可以看到，随着时间的推移，烟雾的可见范围越来越大（设定可见浓度为0.3mg/cm^3），且浓度从扩散中心点向四周呈递减趋势，这与我们的生活经验高度一致；在扩散时间为0s（即扩散未开始）时，空间区域内烟雾浓度均为零，这也与我们设定的初始条件相一致。以上事实初步验证了理论推导结果的正确性。

除此之外，为与后续的实验过程相对应，还记录了特定测量点(x_0/2, 0, z_0)处的烟雾浓度随时间的变化情况，并根据换算公式
$$
ppm = \frac{C \times 10^{6}}{\rho}
$$
将原有的气雾质量浓度C（单位：g/cm^3）换算为体积比ppm（百万分之一），得到的结果如下：

可以看到，由于该处离扩散中心点较远，此处烟雾浓度大概在扩散时间20s后达到基本稳定，稳定时的烟雾浓度值约为380ppm。

3 数值仿真模拟

3.1 数值分析方法：有限元法

由于严格的解析求解对于方程的复杂度以及初始条件与边界条件等都具有一定的要求，第2章所涉及到的理论分析也仅仅是围绕简化后的烟雾扩散问题展开。相比于真实的场景，这样的理论分析忽略了很多现实环境中的复杂变量（比如对流）。为使得我们的研究更具有实用价值，我们需要对于更加泛化的数学物理模型，采用数值求解的方法去模拟真实的扩散情况。

在这一部分，我们主要考虑到，2.3节提到的初始条件中烟雾源处视为单点并采用δ函数进行数学表示的做法过于理想化。与实验相对应的，我们采用的烟雾发生器的烟雾出口也不仅仅是一个单点，而是一个半径约为0.6cm的圆面。因此，我们希望通过有限元这一数值求解方法，对于初始条件迭代后的问题进行数值模拟分析。

图3.1 烟雾发生器烟雾出口孔径测量

基于有限元法的基本思想，我们主要按照如下的方法步骤进行数值分析：

首先进行区域剖分，将长方体研究区域Ω划分为有限个小单元，用于构建离散有限元模型。在划分单元时，采用长方体网格剖分方式，将区域划分为n_x * n_y * n_z个小立方体。每个小立方体单元有8个顶点，因此共有N = (n_x+1) * (n_y+1) * (n_z+1)个节点。对这些节点按列优先顺序进行编号，依次为1、2、…、N，同时对于各单元按体积分块顺序编号，编号为1、2、…、E。

接下来需要选择合适的单元基函数以表示单元内的浓度分布。在此我们选取基函数Φ_ i(x,y,z)近似表示浓度：
$$
\rho(x,\ y,\ z,\ t) \approx \sum_{i = 1}^{N}{c_{i}(t)\phi_{i}(x,y,z)}（c_{i}(t)表示节点i处的浓度值）
$$
并采用线性插值基函数（梯形单元），即单元内任意一点的浓度由8个顶点的浓度值线性插值（符号根据顶点位置的坐标取正或负）：
$$
\phi_{i}(x,y,z) = \frac{1}{8}(1 \pm x)(1 \pm y)(1 \pm z),\ i = 1,2,\ldots,8
$$
值得注意的是，在单元内，基函数在对应节点取值为1，其余节点取值为0，同时基函数也满足局部支撑性，即在单元外为零。

基于划分好的单元和选取的单元基函数，可以写出对应的单元积分表达式。通过对单元内的质量矩阵、刚度矩阵和源项向量进行积分，单元积分表达式应使得方程余量最小、本质边界条件余量最小、自然边界条件余量最小。将扩散方程乘以任意测试函数Φ_j(x,y,z)后，对区域Ω积分即可得到单元积分表达式：

$$
\int_{\Omega}^{}{\phi_{j}\frac{\partial\rho}{\partial t}\mathbb{d}\Omega} = D\int_{\Omega}^{}{\phi_{j}\nabla^{2}\rho\mathbb{d}\Omega} + \int_{\Omega}^{}{\phi_{j}F\mathbb{d}\Omega}
$$
该表达式中共有三项：

• 时间导数项（差分离散）：
  $$
  \int_{\Omega}^{}{\phi_{j}\frac{\partial\rho}{\partial t}\mathbb{d}\Omega} \approx \int_{\Omega}^{}{\phi_{j}\frac{\rho^{n + 1} - \rho^{n}}{\mathrm{\Delta}t}\mathbb{d}\Omega}
  $$

• 空间扩散项：通过分部积分，将扩散项转化为：

$$
\int_{\Omega}^{}{\phi_{j}\nabla^{2}\rho\mathbb{d}\Omega} = - \int_{\Omega}^{}{\nabla\phi_{j}\nabla\rho d\Omega} + \int_{\partial\Omega}^{}{\phi_{j}\nabla\rho \bullet \overrightarrow{n}dS}
$$

在齐次边界条件下，边界项为0。

• 离散形式：代入插值函数近似：

$$
\int_{\Omega}^{}{\phi_{j}\frac{\rho^{n + 1} - \rho^{n}}{\mathrm{\Delta}t}\mathbb{d}\Omega} = D\int_{\Omega}^{}{\nabla\phi_{j}\nabla(\sum_{i = 1}^{N}{C_{i}^{n + 1}\phi_{i}})\mathbb{d}\Omega} + \int_{\Omega}^{}{\phi_{j}F\mathbb{d}\Omega}
$$

基于列出的单元积分表达式，可以对各个单元进行分析，计算其对应的质量矩阵、刚度矩阵与源项向量，从而建立有限元方程。将离散化方程写为矩阵形式：

$$
M\frac{c^{n + 1} - c^{n}}{\Delta t} + Kc^{n + 1} = F
$$
其中：

• 质量矩阵：
  $$
  M_{\mathbb{i}j} = \int_{\Omega}^{}{\phi_{i}\phi_{j}d\Omega}
  $$

• 刚度矩阵：
  $$
  K_{\mathbb{i}j} = D\int_{\Omega}^{}{\nabla\phi_{i}\nabla\phi_{j}d\Omega}
  $$

• 源项向量：
  $$
  F_{\mathbb{i}} = \int_{\Omega}^{}{F\phi_{i}d\Omega}
  $$

进而可以将单元积分表达式更新为有限元方程：

$$
(M + \Delta tk)c^{n + 1} = Mc^{n} + \Delta tF
$$
将每个单元映射到标准立方体单元（边长为2，中心为原点），可以根据基函数梯度
$$
\nabla\phi_{i} = (\frac{\partial\phi_{i}}{\partial x},\ \frac{\partial\phi_{i}}{\partial y},\ \frac{\partial\phi_{i}}{\partial z})
$$
在这些标准单元上使用高斯积分进行质量矩阵和刚度矩阵的计算：

$$
M_{\mathbb{i}j}^{(e)} = \int_{\Omega_{e}}^{}{\phi_{i}\phi_{j}d\Omega}
$$

$$
K_{\mathbb{i}j}^{(e)} = D\int_{\Omega_{e}}^{}{\nabla\phi_{i}\nabla\phi_{j}d\Omega}
$$

将源项向量
$$
F_{\mathbb{i}}^{(e)} = \int_{\Omega_{e}}^{}{F\phi_{i}d\Omega}
$$
离散化为节点的平均值后，也可以通过插值基函数对其进行积分。

完成单元分析后，可以将单元矩阵的局部节点编号1、2、…、8映射到总体矩阵的全局节点编号，并在每个单元中将其局部质量矩阵M^e、刚度矩阵K^e和源项向量F^e加入到总体矩阵M、K和向量F中，进行总体合成。

在对于总体有限元方程进行求解之前，还需要在总体矩阵方程中施加边界条件。在此我们采用消行修正法：对于齐次边界条件，将其带入总体有限元方程，并将边界上的节点浓度固定置零，同时修改总体矩阵对应行列为单位矩阵，右端项置零；对于非齐次边界条件（即初始条件），在边界节点区域（以扩散源中心（x_0/2, y_0/2, 0）为圆心，半径为0.6cm的圆面）中给定固定浓度值F_0，并对应调整右端向量F（已在单元积分中满足）：
$$
F_{i} \Leftarrow F_{i} - k_{ij}F_{0}
$$
最后，对于总体有限元线性方程

$$
(M + \Delta tk)c^{n + 1} = Mc^{n} + \Delta tF
$$
可以采用迭代方法进行求解：

• 初始条件：
  $$
  c^{0} = 0
  $$

• 时间步进：对于每个时间步，首先计算右端向量
  $$
  Mc^{n} + \Delta tF
  $$
  再采用稀疏矩阵求解器计算下一步的浓度
  $$
  c^{n + 1}
  $$

• 迭代终止：当时间达到设定终止时间T时，终止求解过程。

3.2 基于Python的有限元数值仿真

基于3.1节提出的有限元分析方法，利用Python中的MeshPy依赖库进行网格生成，并编写对应Python代码进行有限元数值仿真的计算与可视化。出于程序运行的稳定性等方面考虑，数值计算时仅选取原长方体区域的四分之一（高度不变，长宽各取一半，坐标原点即为扩散源中心）作为扩散区域。

运行程序（代码详见附录2），得到以下数值仿真结果：

可以看到，有限元仿真Python程序的运行结果能够较好地反映扩散过程中烟雾浓度分布的变化情况，基本与理论分析可视化中展现的趋势一致，在此不再赘述。相较于理论分析的可视化结果而言，有限元仿真的结果更加精细也更加丰富，特别是在截面的浓度分布上，很好地反映了扩散过程的基本特点，但在空间中的浓度分布可能受到了边界条件的过分影响，导致其明显呈现沿容器边缘扩散的趋势，而在容器内部的扩散范围较为有限。

3.3 基于COMSOL仿真软件的有限元数值仿真

尽管对于初始条件进行了一定的修正，但上述的有限元数值模拟仍然缺乏对于对流项的考虑，这可能是导致其结果受边界条件影响较大的原因之一。除此之外，上述的理论分析与数值计算也仅仅对于烟雾扩散过程中的浓度分布变化进行计算与模拟，而缺乏了烟雾扩散速度方面的计算分析，这也受到了解析求解方法的限制。因此，我们希望借助成熟的商业有限元仿真软件COMSOL进行更进一步的数值仿真，以分别模拟常压扩散与负压吸附两种情况下的植物甘油气雾颗粒扩散情况，并得到常压扩散情况下的扩散速度场与压力分布情况，帮助我们进一步完善烟雾扩散的基本模型。

为更加真实地反映实际扩散过程中的复杂物理环境，我们需要在COMSOL软件中指定相应的物理场模型：

• 湍流模型：帮助模拟湍流流动对气体扩散过程的影响，增加流体混合、热量传递和物质扩散等过程的复杂性；

• RANS模型：通过对Navier-Stokes方程进行时间平均，简化了直接求解瞬时湍流方程的复杂度，计算成本相对较低，同时可以较为准确地描述气体扩散过程中湍流对流场的影响，尤其是在稳态流动情况下；

• 低雷诺数k-ω模型：采用湍动能（k）和湍流频率（ω）作为主要变量，能够较好地捕捉近壁区和低雷诺数流动特性，在复杂的流动区域中能更好地描述气体的湍流行为，提供更准确的边界层模拟，从而提高气体扩散的模拟精度。

除此之外，我们还需要使用温度场模型来模拟扩散过程中粒子的运动情况。对于烟雾产生端，采用650~750K的的随机梯度划分以达到不同速度的效果。

为了模拟出烟雾粒子扩散过程中的运动状态，我们还需要模拟出颗粒的流动。由于实际实验中采用的植物甘油气雾本质上是一种固体烟雾颗粒，所以在仿真中，我们选取遵循一阶牛顿规则的粒子，并给予每个颗粒一定的质量（与植物甘油气雾的理化性质相一致，设置其密度为1.26g/cm^3）。

处于仿真计算的复杂度与运行时间考虑，我们仅在x = x_0/2处的垂直截面区域内（50cm*60cm）进行仿真，并设置仿真时间为6秒、时间步长为0.1秒、网格单元大小为1cm。同时与实验环境保持一致，温度设置为25℃（298.15K）。

考虑到需要对于常压扩散与负压吸附两种情形分别进行仿真，还需要对于仿真边界（壁）进行一定的设置：

• 常压扩散：将各边设置为壁，壁条件设置为满散射；

• 负压吸附：将上边界设置为出口，边界条件设置为静压力（-0.02MPa）

完成了以上的所有设置流程，可以点击”计算”运行有限元仿真，仿真结果如下所示：

可以看到，通过COMSOL软件进行有限元数值仿真，可以有效划分三角网格并准确模拟出常压扩散情况下的扩散速度场及其强度分布情况，同时仿真得到的粒子运动轨迹大体趋势也与先前的解析求解与有限元数值计算较为吻合。仿真结果成功验证了先前理论模型推导的有效性，并在扩散速度方面进行了有力补充。

4 实验现象验证

4.1 预实验：常压封闭环境下的水雾扩散

在进行理论模型的推导与演算前，为更准确地把握初始条件、边界条件以及部分环境因素的设置，并对于烟雾扩散的具体现象进行初步了解，我们首先基于手头已有的实验材料进行了定性观察的预实验。在预实验中，我们使用香薰机作为烟雾发生装置，但由于装置限制，仅采用纯净水作为烟雾发生原料以产生水雾，通过观察水雾扩散现象来初步把握扩散过程的基本特征。

{width=”1.466816491688539in”
height=”1.9565212160979877in”}

基于以上预实验装置，搭建如下的实验装置场景：

启动香薰机，开始产生水雾，并通过长方体容器底部开孔向容器内部区域空间进行烟雾扩散。设置扩散时间为10分钟，拍摄视频记录扩散实验现象。

可以看到，随着扩散时间的推移，容器内聚集的水雾越来越多，且越靠近底部、越靠近扩散源中心，水雾的浓度也就越高（从现象来看，可以认为水雾团越密集的地方浓度越高），且当扩散时间达到一定程度时，容器内的水雾逐渐达到饱和，开始向容器底部沉降，并在容器底部呈聚集态势。

预实验的定性观察帮助我们对于烟雾扩散的基本现象有了宏观上的认知，并通过实物实验的方式让我们对于研究对象以及研究区域有了更加具象的认识，这对于我们后续理论模型推导与数值计算仿真的有序推进打下了良好的基础。但定性实验只能帮助我们认识现象，要对于理论模型的研究成果进行更加严谨的实证评估，还需要进行定量化的实验与测量。

4.2 正式实验：常压封闭环境下的植物甘油气雾扩散

在完成了理论模型的推导以及解析求解与数值仿真计算之后，我们还希望通过实际定量实验的测量结果来验证理论模型的合理性。为进行定量化实验，我们将预实验中的部分实验装置进行了替换，以严格定量控制烟雾的产生过程并实时测量特定测量点位处的烟雾浓度，从而与理论计算与仿真结果进行比照。

对于烟雾发生器而言，为达到定量产生植物甘油气雾的效果，选取了Selens手持烟雾发生器作为新的实验器材，该款烟雾发生器可以产生固定强度的植物甘油气雾，查阅其产品技术文档可知，其固定产生烟雾的强度约为F_0 ≈ 5g/cm^3（与先前理论推导采用相同符号）；同时为定量测量定点(x_0/2,0,z_0)处（与2.5节一致）的烟雾浓度，选取了MQ-2烟雾浓度传感器进行测量，通过STM32F1开发板对传感器进行供电，同时将数据传输至PC端并实时在屏幕上显示测量到的烟雾浓度值（单位：ppm，与g/cm^3的换算关系在2.5节已提及）。

需要注意的是，由于该款传感器读取烟雾浓度的过程是渐变而非突变，因此在使用前需要一定的预热（初始化）时间，且初始化完成后，在还未开始烟雾扩散实验时，读取到的烟雾浓度也并不为零，因此在后续数据处理时还需要移除对应的偏移量，实测平均偏移量在180ppm左右。

基于以上实验装置，搭建如下的实验装置场景：

对于该款手持烟雾发生器而言，需要手动长按开关才可持续发生烟雾，且该产品具有安全保护机制，最多只能连续产生30s的植物甘油气雾。启动装置后即开始进行烟雾扩散实验，烟雾通过长方体容器底部开孔向容器内部区域空间不断进行扩散。拍摄视频记录扩散实验现象与扩散过程中屏幕上显示的烟雾浓度的实时变化数值及曲线。

可以发现，随着扩散时间的增加，容器内聚集的植物甘油气雾也越来越多，且扩散的过程与理论解析求解可视化以及有限元数值仿真展示的烟雾扩散趋势基本一致。同时，在常压状态下，相较于水雾扩散而言，容器内聚集的植物甘油气雾在停止发生一段时间后仍然没有明显消散，这是由于植物甘油独特的理化性质所导致的，也比水雾扩散实验的结果更加接近香烟烟雾扩散的实际情况。

将烟雾传感器读取到的烟雾浓度数值实时传输到PC端保存，利用MATLAB进行可视化并与理论解析求解计算的定点测量结果进行比照，结果如下：

可以看到，通过烟雾传感器测量得到的烟雾浓度同样在20s左右达到基本稳定，且两者达到稳定后的浓度数值基本接近（在5%的误差允许范围内，约为380ppm），这说明我们的理论模型能够较好地解释实验测量结果；但比较而言，相比起实际测量的数据，理论求解得到的浓度在达到稳定前均高于实际测量结果，这一方面可能由于理论计算时进行了部分近似假设造成了偏差，另一方面也可能是实际测量时在传感器测量与数据读取传输过程中产生了一定的延迟。

通过理论求解计算结果与实验测量观察结果的比照可以发现，我们建立的理论模型对于烟雾扩散情况的模拟效果无论从现象（可视化效果）上还是从浓度数值上都能对于实际的烟雾扩散现象进行较好的模拟，这充分证明了我们理论解析求解模型研究的有效性。

4.3 正式实验：封闭环境中负压作用下的植物甘油气雾扩散

在正式实验阶段，除了在常压条件下进行烟雾扩散实验以验证烟雾扩散理论模型的准确性之外，为了后续进一步地提出可能的烟雾吸附方案并验证其有效性，我们还在负压环境下进行了对应的植物甘油气雾扩散实验，为后续开发提供参考。

在负压源的选取上，我们采用了涵道风机来营造负压环境。该款涵道风机可以通过旋钮调节产生不同档位强度的负压，通过查阅该产品技术文档可知，我们在实验时次啊用的最低档位负压值为P_0≈-0.02MPa，在通过COMSOL有限元仿真对于负压环境下的烟雾扩散进行数值模拟时也采用了相同的参数。其余的实验装置保持不变。

基于常压实验中搭建好的实验装置，在原烟雾发生处的开口对侧开出一个相同大小（半径3cm）的圆孔，并将涵道风机进风口对准该圆孔，完成负压吸附环境下植物甘油气雾扩散实验环境的搭建：

通过STM32F1开发板驱动涵道风机（供电）并通过旋钮将其调至最低档位，并采用与常压环境下实验相同的实验步骤与实验条件，拍摄视频记录扩散实验现象与扩散过程中屏幕上显示的烟雾浓度的实时变化数值及曲线。

可以发现，随着扩散时间的增加，容器内聚集的植物甘油气雾越来越多，但没有像常压状态下一样出现大量烟雾的堆积与沉降现象，且整体扩散方向受负压源牵引作用明显，这点与COMSOL有限元仿真得到的烟雾粒子运动轨迹基本一致；同时可以注意到，同样是在停止烟雾发生30秒后（期间持续进行负压吸附），在负压吸附的持续作用下容器内聚集的植物甘油气雾迅速地被排出容器外，而不像常压环境下的扩散那样会长时间堆积在容器内难以消散，这也充分验证了负压吸附这一解决方案的可行性与有效性，为后续的开发应用提供了良好的实验基础。

将烟雾传感器读取到的烟雾浓度数值实时传输到PC端保存，利用MATLAB进行可视化并与理论解析求解计算的定点测量结果进行比照，结果如下：

可以看到，相较于常压环境中的烟雾扩散，在负压作用下测量点处的植物甘油气雾浓度达到基本稳定的速度明显变快（大概在10s左右），且稳定后测得的烟雾浓度数值（100ppm左右）相较于常压下的稳定测量结果（380ppm左右）有明显下降，这样的数值结果在实验现象的基础上更进一步地验证了负压吸附这一解决方案的可行性与有效性，同时为后续的开发应用提供了良好的实验数据基础。

5 研究结论应用

5.1 研究结论

基于以上对于植物甘油气雾扩散问题在解析求解、数值模拟以及实验测量方面的深入研究，针对1.2节提出的研究目标，可以给出如下三点主要的研究结论：

• 在扩散过程中不同时刻的浓度分布方面，针对常压下密闭环境中的烟雾扩散现象，通过理论计算得到了不同时刻空间内各点处的浓度分布，并结合实验测试结果，印证了理论分析的高精度，在20s后空间内浓度分布基本稳定；

• 在扩散过程中烟雾的扩散速度方面，通过COMSOL有限元仿真软件对常压下密闭环境中的烟雾扩散进行数值模拟，得到扩散处于基本稳定状态下空间截面内速度场的分布，为分析浓度分布变化趋势提供参考；

• 在扩散过程中的烟雾扩散效果与除烟效果方面，将常压下和负压吸附下的COMSOL有限元仿真结果与实际实验现象进行对照，可以发现理论分析、数值求解以及实际实验三者的扩散现象结果基本一致，且在引入负压对流后，稳定后的烟雾浓度整体较常压时明显降低，且在烟雾发生停止后能够迅速清除残余的烟雾，具有良好的除烟效果。

除此之外，在这三点研究结论的基础上，我们还通过理论、仿真与实验的有机结合，得到了植物甘油气雾扩散的一般规律与数学物理模型，可以通过该理论模型对于实际的烟雾扩散情况进行定量化的数值模拟分析，从而对于烟雾扩散过程进行更精准的把控与必要的干预管理，并通过多种手段对该理论模型的有效性进行了充分验证，这也为我们后续应用产品的开发打下了良好的理论与实验基础。

5.2 研究成果应用：室内烟雾浓度控制系统

基于数值分析结果与研究结论，我们开发了一套简易的闭环控制系统作为原型样机，通过实时检测监测点处的烟雾浓度，实现对于风机吸附功率的灵活动态调控，以优化用户体验。

在该系统的开发方面，采用STM32F1开发板作为系统主控，通过ADC实时读取MQ-2传感器返回的烟雾浓度数值，并通过PID闭环控制实时调控给涵道风机供电的PWM占空比，从而实现对于风机吸附功率的实施调控。

通过运行该闭环控制系统并测试其烟雾浓度控制效果，可以发现该系统对于烟雾有着较高的灵敏度，并在检测到高浓度烟雾后能够迅速调动涵道电机作出反应，快速吸附容器内的植物甘油气雾并始终维持空间内的烟雾浓度低于150ppm，说明样机能够初步满足用户的基本除烟需求。

6 课程收获感悟

通过本次数学物理方法课程的学习，以及烟雾扩散课程项目的开展，我们通过在实际项目中运用课程中学习到的理论知识，在扩散过程中烟雾浓度分布的定量化求解过程中，对于数学物理方程的导出、边界与初始条件的确定以及基于变量分离与本征值方法求解级数解的过程有了更加深刻的认识与理解，同时通过课程中介绍的数值求解有限元方法对于更加泛化的情况进行了数值模拟与讨论。在解析与数值求解的理论分析基础上，我们还利用成熟的COMSOL有限元仿真软件对于烟雾扩散的粒子运动轨迹以及速度场分布进行数值仿真以作为对我们理论模型的补充，并通过定量化的实验测量进一步验证了理论模型的有效性，为我们后续进行实际的产品开发与应用提供了良好的理论与实验基础。

在本次课程项目报告的最后，再次感谢老师的耐心指导以及提出的宝贵意见，也希望我们能够将数学物理方法教给我们的定量化分析思维带入到未来的项目研究与学习工作中，对我们未来的项目开展与人生成长起到帮助。

参考文献

[1] Gao, Naiping and Jianlei Niu. “Modeling particle dispersion and deposition in indoor environments.” Atmospheric Environment (Oxford, England : 1994) 41 (2007): 3862 - 3876.

[2] Hoegg U. R. (1972). Cigarette smoke in closed spaces. Environmental health perspectives, 2, 117–128.

[3] Kuga K, Ito K, Yoo S-J, et al. First- and second-hand smoke dispersion analysis from e-cigarettes using a computer-simulated person with a respiratory tract model. Indoor and Built Environment. 2018;27(7):898-916.

[4] Holmberg, S., & Li, Y. (1998). Modelling of the Indoor Environment – Particle Dispersion and Deposition. Indoor Air, 8, 113-122.

[5] Al-sarraf, A.A., Yassin, M.F. & Bouhamra, W. Experimental and computational study of particulate matter of secondhand smoke in indoor environment. Int. J. Environ. Sci. Technol. 12, 73–86 (2015).

[6] 陈占秀,陈冠益,王艳,等.丙三醇与1,6-己二醇混合物降温凝固过程的分子动力学模拟[J].化工学报,2013,64(07):2316-2321.

[7] 沈玉峰,孟强龙.基于物理模型的烟雾扩散模拟[J].电子技术,2011,38(03):61-63.

[8] 李强,普小云.用毛细管成像法测量液相扩散系数——等折射率薄层测量方法[J].物理学报,2013,62(09):185-191.

[9] 杜青青.室内烟雾扩散模拟技术研究与实现[D].北京工业大学,2018.

[10] 梁昆淼编.数学物理方法[M].北京：高等教育出版社,1978:593页.

附录

1 MATLAB解析求解可视化代码

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

% 参数定义
x0 = 40; % 长方体的长度
y0 = 50; % 长方体的宽度
z0 = 60; % 长方体的高度
F0 = 0.5;  % 源项强度 g/cm3
D = 0.0962;   % 扩散系数
x_s = x0 / 2; % 点源的 x 坐标
y_s = y0 / 2; % 点源的 y 坐标
z_s = 0; % 点源的 z 坐标
t_max = 30 * 40; % 最大时间
num_time_steps = 31; % 时间步数
num_terms = 10; % 级数项数
% 创建空间网格
x = linspace(0, x0, 50);
y = linspace(0, y0, 50);
z = linspace(0, z0, 50);
[X, Y, Z] = meshgrid(x, y, z);
time_steps = linspace(0, t_max, num_time_steps);   % 时间步长
rho = zeros(length(x), length(y), length(z), num_time_steps);   % 计算解的函数
% 测量点坐标 (x0/2, 0, z0)
x_measure = x0 / 2;
y_measure = 0;
z_measure = z0;
rho_measure = zeros(1, num_time_steps);   % 初始化存储测量点浓度的数组
for t_idx = 1:num_time_steps
    t = time_steps(t_idx);  % 当前时间
    % 级数求解
    sum_term = 0;
    for l = 1:num_terms
        for m = 1:num_terms
            for n = 1:num_terms
                % 计算级数各项
                lambda = (l^2 * pi^2 / x0^2) + (m^2 * pi^2 / y0^2) + (n^2 * pi^2 / z0^2);
                term = cos(l * pi * X / x0) .* cos(m * pi * Y / y0) .* cos(n * pi * Z / z0) .* ...
                       cos(l * pi * x_s / x0) .* cos(m * pi * y_s / y0) .* cos(n * pi * z_s / z0) .* ...
                       (1 - exp(-D * lambda * t)) ./ lambda;
                sum_term = sum_term + term;
            end
        end
    end
    rho(:, :, :, t_idx) = (F0 / (x0 * y0 * z0 * D)) * sum_term;   % 乘上源强度项和常数因子
    % 计算测量点处的浓度
    rho_measure(t_idx) = interp3(X, Y, Z, rho(:, :, :, t_idx), x_measure, y_measure, z_measure); 
    rho_measure(t_idx) = rho_measure(t_idx) * 10^6 / 1.26;   % 单位换算：g/cm^3 -> ppm
    disp(t_idx)
end
time_steps = time_steps / 40;
% 创建视频对象
video_file = 'fog_diffusion.avi'; 
v = VideoWriter(video_file);
v.FrameRate = 5; % 设置帧率
open(v); % 打开视频文件
% 绘制动态可视化
figure;
for t_select = 1:num_time_steps % 遍历时间步
    value = rho(:, :, :, t_select);
    % 筛选出非零值及其对应的坐标
    non_zero_indices = value >= 0.0003; % 设定阈值：可见浓度0.3mg/cm3
    X_nonzero = X(non_zero_indices);
    Y_nonzero = Y(non_zero_indices);
    Z_nonzero = Z(non_zero_indices);
    Value_nonzero = value(non_zero_indices);
    % 绘制非零值的散点图
    scatter3(X_nonzero, Y_nonzero, Z_nonzero, 20, Value_nonzero, 'filled'); % 绘制散点图
    colormap jet; % 设置颜色图
    colorbar; % 显示颜色条
    xlabel('x');
    ylabel('y');
    zlabel('z');
    title(['3D Scatter Plot of Values at t = ' num2str(time_steps(t_select))]);
    grid on;
    view(3); % 设置三维视角
    % 固定坐标轴长度
    axis([0, x0, 0, y0, 0, z0]); % 设置 x, y, z 的范围
    % 将当前帧保存到视频
    frame = getframe(gcf);
    writeVideo(v, frame);
end
close(v);   % 关闭视频文件
% 绘制测量点浓度随时间的变化
figure;
plot(time_steps, rho_measure, '-o', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('\rho at (x0/2, 0, z0) (ppm)');
title('Concentration at Measurement Point over Time');
grid on;
disp(['视频已保存为: ', video_file]);
% 将测量点浓度数据保存到txt文件
output_file = 'rho_measurement_data.txt';
fileID = fopen(output_file, 'w');
fprintf(fileID, 'Time(s)\tConcentration(ppm)\n'); % 表头
for i = 1:length(time_steps)
    fprintf(fileID, '%.2f\t%.6f\n', time_steps(i), rho_measure(i));
end
fclose(fileID);
disp(['测量点数据已保存为: ', output_file]);

2 Python有限元仿真代码

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
102
103
104
105
106
107

import numpy as np
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
import matplotlib.animation as animation

# 网格生成函数
def generate_mesh(x0, y0, z0, nx, ny, nz):
    from meshpy.tet import MeshInfo, build
    points = [(0, 0, 0), (x0, 0, 0), (x0, y0, 0), (0, y0, 0),(0, 0, z0), (x0, 0, z0), (x0, y0, z0), (0, y0, z0)]
    facets = [[0, 1, 2, 3], [4, 5, 6, 7], [0, 1, 5, 4], [1, 2, 6, 5], [2, 3, 7, 6], [3, 0, 4, 7]]
    mesh_info = MeshInfo()
    mesh_info.set_points(points)
    mesh_info.set_facets(facets)
    mesh = build(mesh_info, max_volume=(x0/nx) * (y0/ny) * (z0/nz))
    return mesh
# 全局矩阵组装
def assemble_global_matrices(mesh, D):
    num_nodes = len(mesh.points)
    M_global = lil_matrix((num_nodes, num_nodes))
    K_global = lil_matrix((num_nodes, num_nodes))
    for element in mesh.elements:
        node_ids = element
        vertices = np.array([mesh.points[i] for i in node_ids])
        B, volume = compute_B_matrix(vertices)
        K_element = D * volume * (B.T @ B)
        M_element = volume * np.eye(4)
        for i in range(4):
            for j in range(4):
                M_global[node_ids[i], node_ids[j]] += M_element[i, j]
                K_global[node_ids[i], node_ids[j]] += K_element[i, j]
    boundary_indices = find_boundary_nodes(mesh)
    apply_boundary_conditions(M_global, K_global, boundary_indices)
    return M_global, K_global
def compute_B_matrix(vertices):
    A = np.ones((4, 4))
    A[:, 1:] = vertices
    volume = abs(np.linalg.det(A)) / 6.0
    B = np.zeros((3, 4))
    for i in range(4):
        sub_matrix = np.delete(A, i, axis=0)
        coeffs = np.linalg.det(sub_matrix[:, 1:])
        B[:, i] = coeffs
    B /= (6 * volume)
    return B, volume
def find_boundary_nodes(mesh):
    boundary_indices = set()
    for facet in mesh.facets:
        boundary_indices.update(facet)
    return list(boundary_indices)
def apply_boundary_conditions(M_global, K_global, boundary_indices):
    for idx in boundary_indices:
        M_global[idx, :] = 0
        M_global[:, idx] = 0
        K_global[idx, :] = 0
        M_global[idx, idx] = 1
        K_global[idx, idx] = 1
# 时间步进求解
def solve_diffusion(mesh, M_global, K_global, F0, D, t_max, num_time_steps):
    num_nodes = len(mesh.points)
    rho = np.zeros((num_nodes, num_time_steps))
    dt = t_max / (num_time_steps - 1)
    center = np.array([x0/2, y0/2, 0])
    radius = 0.6
    for t_idx in range(1, num_time_steps):
        A = M_global + dt * K_global
        F = F0 * np.ones(num_nodes)
        for idx, point in enumerate(mesh.points):
            if np.linalg.norm(point - center) < radius:
                F[idx] += F0
        b = M_global @ rho[:, t_idx-1] + dt * F
        rho[:, t_idx] = spsolve(A, b)
    return rho
# 可视化与动画
def create_animation(mesh, rho, time_steps, output_file):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    points = np.array(mesh.points)
    def update_plot(frame_idx):
        ax.clear()
        rho_t = rho[:, frame_idx]
        norm = Normalize(vmin=0.005, vmax=0.01)
        ax.scatter(points[:, 0], points[:, 1], points[:, 2], c=rho_t, cmap='jet', norm=norm, s=10)
        ax.set_title(f"Concentration at t = {time_steps[frame_idx]:.2f}")
        ax.set_xlim(0, x0)
        ax.set_ylim(0, y0)
        ax.set_zlim(0, z0)
    ani = animation.FuncAnimation(fig, update_plot, frames=len(time_steps), interval=500)
    ani.save(output_file, writer='ffmpeg', fps=5)
plt.show()

# 主程序
if __name__ == "__main__":
    # 参数
    x0, y0, z0 = 20, 25, 60
    nx, ny, nz = 10, 10, 10
    D = 0.0962
    F0 = 0.5
    t_max = 30
    num_time_steps = 31
    mesh = generate_mesh(x0, y0, z0, nx, ny, nz)   # 网格生成
    M_global, K_global = assemble_global_matrices(mesh, D)   # 矩阵组装
    # 求解
    time_steps = np.linspace(0, t_max, num_time_steps)
    rho = solve_diffusion(mesh, M_global, K_global, F0, D, t_max, num_time_steps)
    create_animation(mesh, rho, time_steps, "diffusion_animation.mp4")   # 生成动画