第5章偏微分方程数值解

第5章

偏微分方程数值解

5.1问题的提出

5.2基本离散化公式

5.3几种常见方程的离散化计算

5.4吸附床传热传质模型中偏微分方程求解实例

目录5.1问题的提出包含有偏导数的微分方程称为偏微分方程。从实际问题中归纳出来的常用偏微分方程可分为三大类：波动方程、热传导方程和调和方程。对于它们特殊的定解条件，有一些解决的解析方法，而且要求方程是线性的、常系数的。但是在实际中碰到的问题却往往要复杂得多，尤其在化工和化学模拟计算中，不仅偏微分方程的形式无一定标准，且边界条件五花八门，方程中的系数随工况改变而改变，想利用解析求解是不可能的。另一方面实际问题的要求不一定需要严格的精确解，只要求达到一定精度，所以就可借助于差分方法来求偏微分方程的数值解。在第4章里，我们介绍了一个套管式换热器稳态的传热问题。如果我们考虑一个动态的传热过程，且不忽略纵向的热传导，就可以得到以下的偏微分方程：5.15.45.35.2总目录5.1问题的提出上面方程中变量的含义如下：通过求解上面的偏微分方程，就可以得到传热管各点温度随时间的变化值，从而确定达到传热平衡所需的时间，为实验测量提供依据。想求解上述方程，就必须首先学会偏微分方程的求解方法，下面我们首先介绍如何对偏微分方程进行离散化的工作，然后再对各类不同的偏微分方程进行求解，我们一般只给出离散化的基本公式及计算方法，对离散化公式的具体推导工作一般不作详细介绍，对这方面感兴趣的读者可自行参考有关数值计算的书籍。

5.15.45.35.2总目录5.2基本离散化公式

在偏微分方程中，自变量都在两个或两个以上，应变量随两个或两个以上的自变量变化而变化。在化工或化学动态模拟方程中，常常有一个自变量是时间，其它的自变量为空间位置。如果只考虑一维空间，则只有两个自变量；如果考虑两维空间，则有3个自变量。一般我们将自变量在时间和空间以一定的间隔进行离散化，则应变量就变成了这些离散变量的函数，以3维空间为例，我们将离散化的应变量表示成，它所表示的真正含义如下:有了以上的定义，对于一阶偏导我们可以利用第四章的欧拉公式直接得出向前欧拉公式：对于时间偏导而言，有时我们常常采用向后欧拉公式，时间的向后欧拉公式如下：5.15.45.35.2总目录5.2基本离散化公式

这样在以后的计算中，得到的是隐式的计算公式，需通过求解线性方程组才能求解。具体的计算过程我们在下面会针对具体的偏微分方程进行讲解。对于二阶偏导，我们可以通过对泰勒展开式处理技术得到下面离散化计算公式：有了以上的离散化公式，就可以进行偏微分方程的数值求解工作。当然，在具体求解时，还会碰到不同的问题，需要区别对待，同时在利用计算机编程计算时也会碰到困难，这些问题我们会通过具体的例子加以说明。

5.15.45.35.2总目录5.3几种常见方程的离散化计算

1、

波动方程其中：为初值条件为边值条件

当该波动方程只提初值条件时，称此方程为波动方程的初值问题，二者均提时，称为波动方程的混合问题。对于初值问题，是已知t=0时，u与依赖于x的函数形式，求解不同位置，不同时刻的u值。而

u是定义在的二元函数，即上半平面的函数。

对于混合问题除初值外，还有边值。是已知初值及x=0及x=l时u依赖于t的函数，求解不同位置x，不同时刻的u值。此时u是定义在的带形区域上的二元函数。如图可以看出初值问题和混合问题的定义域。

5.15.45.35.2总目录5.3几种常见方程的离散化计算

根据5.2节提供的公式，将上面波动方程离散化，得到：(5-1)将式（5-1）进行处理，把（n+1）时刻的变量留在右边，其余放在左边得到：

（5-2）

同时将边界条件和初始条件也离散化，得到：

(5-3)xt0a)初值问题tx0lb)混合问题5.15.45.35.2总目录5.3几种常见方程的离散化计算

这样，由式（5-2），并结合式（5-3），就可以从n时刻的各点u值，计算得到下一时刻的u值，这样层层递推，就可以计算出任意时刻，任意位置的u值。而图5.2则表明了这种层层递推的计算过程，在图5.2中*表示需求u值的点，○表示为了求x点的u

值必须已知u值的点。需要说明的是，在应用式（5-2）进行计算时，初值与边值应当满足相容性条件。由初值得到，由边值得到

，,但在利用式（5-2）进行第一轮计算时，若取n=0，则发现等式右边出现了，这是一个无法计算的值。这时可以利用另一个初值条件算得，这样，可在第一轮计算的时候，取n=1，计算得到，由，递推得到，这样就可由式（5-2）一排一排往上推，计算得到所有希望得到的u值。对于式（5-2）取n=0计算中碰到的，也可利用另一种方法进行计算，解决的办法是将另一个初值条件利用向后欧拉离散化算得，这样利用式（5-2），取n=0就可以得到，取n=1，，和前一种处理方法一样一排一排往上推，计算得到所有希望得到的u值。象这样可以用已知点上函数值直接推出所有点上函数值的格式，称为显式格式。当方程非齐次时，，式（5-1）可写为5.15.45.35.2总目录5.3几种常见方程的离散化计算

当方程是初值问题时，边界条件没有了，由于在t=0时，u与值是已知的，若需要求某的值，只要按“波及原则”多算一些初值，即可推得，所图5-2所示。为了保证差分方程的解在时收敛于原来波动方程的解，要求式（5-2）中等式右边的各项系数均大于0，即：化简得：

而且，可以证明，只要初始条件，边界条件满足一定的光滑性要求，且满足收敛关系式时，差分格式是稳定的。

图5-2x5.15.45.35.2总目录5.3几种常见方程的离散化计算

例5.1:

用数值法求解下面偏微分方程，并写出VB程序。解：首先根据前面的知识，将所求的方程离散化，先假设以下各式：代入微分方程并化简得：（5-3）分析式（5-3）可知，如果知道了某一时刻的各点t，（j=0,1,2….100），就可以求下一时刻的各点温度值。有了以上各式，上面的微分方程就可以求解了。其实这个微分方程，是在不考虑流体本身热传导时的套管传热微分方程。由计算结果可知，当计算的时间序列进行到７２时，传热过程已达到稳态，各点上的温度已不随时间的增加而改变。如果改变套管长度或传热系数，则达到稳态的时间亦会改变。

5.15.45.35.2总目录5.3几种常见方程的离散化计算

2、一维流动传热传导方程的混合问题与波动方程的情形类似，用差商近似代替偏商，可以得到差分方程，以其解作为流动传热传导方程的近似解。一维流动传热传导方程的混合问题：

上面的偏微分方程其实就是在5.1节中提出的偏微分方程，利用5.2节中的离散化公式进行离散化，得到其离散化公式：

5.15.45.35.2总目录5.3几种常见方程的离散化计算

将上式进行处理得到：

（5-4）利用初始条件和边界条件，可以得到零时刻各点的（i=0,1,2,…m）及

，这样就可以利用公式（5-4）计算得到，依次类推，可以得到其它时刻的各点值，所以式（5-4）也是显式格式。只要保证式（5-4）中各项系数大于零，一般情况下，式（5-4）的计算公式是稳定的，可以获得稳定的解。分析式（5-4）可以发现，当为了提高数值精度取适当小的时，最有可能小于零的系数是的系数，若要保证此项系数大于零，此时必须相应地更小，这样，计算量将大大增加，这是显式格式的缺点，为了克服此缺点，下面提出一种隐式格式。偏微分方程在点上进行离散化，且对时间的偏微分采用向后欧拉公式得到原偏微分方程的离散化公式：

（i=1,2….m）

5.15.45.35.2总目录5.3几种常见方程的离散化计算

从图5-3中可见要由初值及边界条件一排一排推上去是不行的，需解线性方程组,同时添上二边界条件：

正好共有m+2个方程，同时有m+2个变量，就能解出n+1排上各点值。（至于线性方程组的求解方法我们在第3章中已作过介绍，请读者自行参照第3章的内容）。这样，每解一个线性方程组，就可以往上推算一排点的u值，虽然引入了方程组的求解，有可能增加计算量，但由于隐式格式无条件稳定，

的取法与无关，可以少计算许多排节点上的u值，相应于显式格式来说，最终反而节省了计算量。

图5-3

5.15.45.35.2总目录5.3几种常见方程的离散化计算

例5.2请计算考虑纵向导热的套管换热器内管各点温度分布微分方程：解：首先根据前面的知识，将所求的方程离散化，先假设以下各式：代入微分方程并化简得：分析式（5-3）可知，如果知道了某一时刻的各点t，（j=0,1,2….10,11），就可以求下一时刻的各点温度值t(j=1,2….10),现在已经知道了零时刻管内各点的温度分布及入口处在任何时刻的温度，如想求下一时刻的温度值，根据上面的离散化计算公式，还需知道在j=11处的温度，这个温度可利用给定的边界条件离散化求得：有了以上各式，上面的微分方程就可以求解了。

5.15.45.35.2总目录5.3几种常见方程的离散化计算

和前面不考虑热传导的情况比较，可以发现温度有细微的变化，如果导热系数足够大，则温度的变化会更大。如导热项的系数为0.2时，其计算公式变为：

由于导热的缘故，已经加热的向前流动的流体却要向后方向进行热传导，从而降低了总体传热效率，使在相同时刻、相同位置点的温度比没有热传导时要低。

前

后

5.15.45.35.2总目录5.3几种常见方程的离散化计算

3、稳态导热/扩散方程

在化工导热及扩散过程中，没有物流的流动，仅靠导热及扩散进行热量及质量的传递。如果此时系统达到稳定状态，也就是说系统中每一个控制单元的各项性质如温度、浓度等不再随时间的改变而改变，系统中的各种性质只与其所处的位置有关，利用化工知识，我们可以得到下面二维、三维的稳态导热或扩散偏微分方程：二维：三维：下面我们主要介绍二维的求解方法，二维的稳态导热或扩散偏微分方程又称调和方程，其方程示意图见图5-4所示。常见有三种边界条件：第一类边界条件：

第二类边界条件：

第三类边界条件：

图5-45.15.45.35.2总目录5.3几种常见方程的离散化计算

在化工中碰到较多的是第一类边界条件，下面以第一类边界条件为例，说明其求解方法：

首先利用5.2中提供的离散化公式（不考虑u中的上标变量n），可得下面离散化公式：

取，经化简得：对于每一个边界内的离散点均可列出这样的五点格式。若中有边界点，用边界值代入。若靠边界很近，也可以看作边界节点，从靠它最近的边界点上的值来取代。由于此计算格式不存在时间上的递推问题，它只是不同空间位置上变量的求解问题，而已知条件仅仅知道边界上的值，这样要求边界内点的值只能通过离散化的偏微分方程来求解，幸好有多少个内节点就有多少个离散化的方程，构成了一个未知数个数与方程个数相等的稀疏方程组，既可直接求解，也可迭代求解，一般用迭代法解比较好。

5.15.45.35.2总目录5.3几种常见方程的离散化计算

下面介绍3种迭代格式：

(1)同步迭代：

(2)异步迭代：

(3)超松弛迭代：当计算范围R为矩阵区域，x方向m等分，y方向n等分，那么最佳松弛因子由数学知识可知，用这些迭代法求解上面的偏微分方程均收敛。

5.15.45.35.2总目录5.3几种常见方程的离散化计算

例5.3

：处于传热平衡状态的某保温，假设其形状为长方体，在x，y两个方向上存在热传导，且导热系数相等，已知边界温度分布如下图所示：

试列出其传热微分方程，并求出各点的温度分布（间隔以），并画出温度分布图。

解：取某一微元进行能量衡算，由于已达传热平衡状态，故可得传导入热量-传导出热量=0

1xy10(1,1)5.15.45.35.2总目录5.3几种常见方程的离散化计算

根据上图的导热表达图可得：化简得：由已知条件可知，，则有：

此偏微分方程和前面介绍的调和方程一致，可用五点格式同步迭代计算，其中需要计算的内点共有16点。也可列出16个线性方程，组成方程组利用第3章介绍的方法进行求解，本书介绍利用五点格式同步迭代计算的VB程序

[+（）]－[＋]=05.15.45.35.2总目录5.3几种常见方程的离散化计算

计算结果：

图示：

x

y00.20.40.60.81.003040.873.2127.2202.83000.284104.53134.89178.09234.133000.4138158.45183.73216.13255.653000.6192207.54225.46247.06272.343000.8246254.26263.5274.3286.663001.03003003003003003005.15.45.35.2总目录VB调用5.4吸附床传热传质模型中偏微分方程求解实例

5.4.1基本设定及假设

1.吸附器结构参数的设定

热流体

上图所示的是套筒式吸附器，该吸附器的有效长度为L，其有效内径为D,环隙宽度为δ，吸附器壁厚为δb。导热流体通过环隙将热量传入或传出吸附器，吸附质通过吸附器上端的小管进入或离开吸附器。

图5-5吸附器结构示意图

5.15.45.35.2总目录5.4吸附床传热传质模型中偏微分方程求解实例

2.吸附床外流体传热的一些基本假设

１．忽略流体在环隙宽度δ上的温度梯度；２．忽略热损失；3.忽略吸附器壁厚δb上的温度梯度，用集中参数法求取吸附器壁面温度。

３.吸附床内传热传质的一些基本假设

１.

吸附床内的吸附质气体处于气滞状态；２.

忽略蒸发器、冷凝器和吸附床之间的压力差；３.

吸附床内各计算微元内达到吸附平衡。吸附量可利用回归方程计算；４.

吸附热利用微分吸附热，随吸附量和吸附温度的改变而改变；热采用有效比热，亦随温度改变，但在计算微元内，可认为是常数；５.床层活性炭导热系数采用当量导热系数，其具体数值利用实验测量值。

5.15.45.35.2总目录5.4吸附床传热传质模型中偏微分方程求解实例

5.4.2流体传热模型的建立

在轴方向上取一环隙微元，见图５－７，作能量分析如下：１.流体通过流动流入环隙微元的能量去

其中f为流体的密度uf

为环隙的流体速度，

Sf为环隙的横截面积，Cpf

为流体的比热。2.流体通过流动流出环隙微元的能量3.流体热传导在x

处的热量导入

图５－6流体传热微元模型

5.15.45.35.2总目录5.4吸附床传热传质模型中偏微分方程求解实例

4.流体热传导在x+x处的热量导入

5.微元体传递给吸附床的热量qt

6.微元体内的能量变化率

7总能量平衡方程

其中,为流体的横截面积。

5.15.45.35.2总目录5.4吸附床传热传质模型中偏微分方程求解实例

5.4.3吸附床内吸附剂传热传质模型的建立

吸附床内发生着热量和质量的传递，但质量的传递是建立在热量传递基础上的，故只要建立热量传递方程，就可以根据平衡吸附量方程求出各处的吸附量。吸附床内的热量传递主要以热传导为主，既有经向的热传导，也有轴向的热传导，为了便于建模分析，我们选取如图5.8的吸附床微元体，具体分析如下：

x+x

r

r+r

xx

图5.7吸附床内传热传质微元体

1.轴向热量导入:5.15.45.35.2总目录5.4吸附床传热传质模型中偏微分方程求解实例

2.轴向热量导出

3.经向热量导入

4.径向热量导出

5.微元体内的能量变化率

其中,为吸附床层内的有效比热。6.总能量平衡方程

其中

5.15.45.35.2总目录5.4吸附床传热传质模型中偏微分方程求解实例

5.4.4吸附器壁面温度轴向分布方程

和前面的分析方法一样，通过微元能量平衡方程并作适当化简可得：

其中为吸附器壁面的横截面积。

5.15.45.35.2总目录5.4吸附床传热传质模型中偏微分方程求解实例

5.4.5吸附器内/外无因子化方程根据前面推导已经得到的方程，并对变量作以下无因子化处理：

通过以上的无因子化处理，可得吸附器内、外无因子化传热传质方程如下：