Matlab/Simulink建模详解：一阶时变偏微分方程的求解

问题的定义

这一次日笃小编来教大家如何在simulnk里面，求解偏微分方程（Partial Differential Equation-PDE）。这次给大家讲解偏微分方程。跟ODE不同的是，PDE的变量不仅是时间的函数，而且也是空间的函数。举个简单的例子：炒菜的锅的锅壁温度，温度这个变量不仅仅是时间的参数（不同时候，锅壁同一个点的温度不一样），而且也是空间的参数（同一时刻，锅壁不同位置的温度不一样）。求解PDE的时候可以直接求解，也可以先转化成ODE，再用求解ODE的方法求解。这里，小编先将PDE转化成ODE再求解。

主要方程的定义

先给出三个相互关联的方程（前两个是PDE，第三个是ODE）：

其中Tf是hf关系式，A和K是常数

其中Te是he关系式，A和K是常数

M和K是常数

PDE转化为ODE

再将前面两个PDE转化为ODE：

这样我们就得到了三个ODE，运用显示的ODE求解方法，通过简单的欧拉展开，我们可以得到三个变量的求解显示表达式：

再将前面两个PDE转化为ODE：

这样我们就得到了三个ODE，运用显示的ODE求解方法，通过简单的欧拉展开，我们可以得到三个变量的求解显示表达式：

边界条件和起始条件的讨论和设定

这三个表达式中，只有等式左边的参数是未知数，等式右边的所有变量都是已知的。要求解这三个表达式，我们需要设定一些条件：起始条件和边界条件。所谓起始条件就是开始模拟的那一个时刻，状态量的数值，在这里我们需要给定Tw,hf,he三个变量在Z1，Z2位置对应的数值。如果模拟的持续时间很长，起始条件在一定范围内一般对整个模拟的影响不是特别大，比如计算电池充放电过程中电池内部温度的变化（一般持续10个小时以上）；但是如果模拟持续时间短，起始条件的设定就显得尤为重要，比如汽车碰撞（一般在0.5秒以内就能完成），或者安全气囊的起爆过程（一般在0.015秒内就能完成）。边界条件指的是有限元或者有限体积法将一个平面或者立体划分成很多个小单元，这些单元在物理边界的位置，状态量随时间变化的关系。比如：要计算一个炒菜锅炒菜过程中锅壁温度的分布，我们需要知道锅底火焰再各个接触面的温度随时间变化的关系；再比如要计算一个轮胎的压到一颗石子过程中的变形过程，我们需要知道这个石头在轮胎表面产生的轮廓随时间变化的关系；再比如我们要计算一个空调房各个角落的温度，我们需要知道，空调出风口的风速大小方向、温度、湿度随时间变化的关系，以及房间门口地下缝隙跟室外空气交换的速度以及室外空气温度和湿度随时间变化的关系。

下面小编设定一下上面式子的起始条件(Initial Condition - I.C)和边界条件(Boundary Condition - B.C)：

I.C:

Tw = [377.540533331317 386.026288241315 394.882816943803 403.780234260442 412.528927823382 421.020618818433 429.194964414357 437.019737325395 444.479243081580 451.567215232647 458.281835132688 464.622699931517 464.045356725446 461.548572973745 460.533446196397 459.951339710425 459.573192331990 459.311655878501 459.125265766771 458.992877081766 458.900967710798 458.850165061973 458.921507977048 459.670267128948 488.518351484883 503.038359452690 516.356161826663 528.545395038828 539.687125046220 549.862685009442]

hf = [335879.230662688 383800.108882945 430054.713708057 474539.317712384 517232.789150323 558163.101428462 597387.616256567 634981.549133264 671031.215823456 705630.227842953 738877.802385151 770878.391400899 804774.697027078 841552.709305229 880672.392971595 922052.212463106 965712.639000156 1011715.66565823 1060145.65583764 1111101.06457330 1164632.21059020 1220683.53929694 1279315.98281919 1340335.98261620 1390832.64999538 1436955.33525927 1479057.18575508 1517473.69348631 1552518.42857505 1584481.89301558]

he = [502113.781884773 508072.678267337 513859.836263859 519445.766992407 524817.943825280 529973.814790241 534916.758091738 539653.701008431 544193.726929246 548547.259765169 552725.604497504 556740.743076404 560605.289626758 564698.772122112 569140.263508662 573864.546505107 578861.774480850 584134.419639887 589689.968567105 595538.609396464 601692.232036153 608156.913219795 614925.943905113 622006.682551548 629375.754089632 635473.976895104 641043.976376731 646128.400323336 650767.764378679 654999.936793652]

B.C :

Simulink建模

上述三个方程式代表的是有限体积法中一个小单元的方程，小编把很多个小单元依次连接起来，重复求解上述方程就能够得到解。对于有限体积法，在方程求解器设计的过程中，一开始并不知道需要离散多精密才能得到理想的解的精度，所以一开始不知道要求解多少组这样的三个方程。如果把一组三个方程的求解写在一个simulink的block里面，如果要用100组的话，小编需要复制粘贴100次，不是很显示，并且连线都会连得手抽筋，小编选择用simulink里面的matlab function模块，这个模块就是把matlab的m文件与simulink对接起来，能够用m文件最大的好处就是可以写for循环，这样的话，只需要改变for循环的循环次数，就能够很容易的实现不同离散程度方程的求解。下面小编把最终的simulink文件的截图奉上：

起始条件是放在simulink的记忆block里面：

而边界条件则通过查表的方式，实时输入到matlab function里面，这里只需要事先把随时间变化的边界条件存成一个随时间变化的一维向量，保证时间向量的长度跟边界条件向量的长度保持一致就行，然后利用simulink的1D(一维)表格，输入是时间，输出是对应的边界条件（如下）：

然后就设定模拟时间,在simulink的设置里面，选定离散的一阶欧拉求解器，就可以开始模拟了。

Simulink仿真结果

下面给出三个状态变量随时间变化的曲线：