第5章偏微分方程数值解

第5章偏微分方程数值解5.1问题的提出5.2基本离散化公式5.3几种常见方程的离散化计算5.4吸附床传热传质模型中偏微分方程求解实例目录5.1问题的提出包含有偏导数的微分方程称为偏微分方程。从实际问题中归纳出来的常用偏微分方程可分为三大类：波动方程、热传导方程和调和方程。对于它们特殊的定解条件，有一些解决的解析方法，而且要求方程是线性的、常系数的。但是在实际中碰到的问题却往往要复杂得多，尤其在化工和化学模拟计算中，不仅偏微分方程的形式无一定标准，且边界条件五花八门，方程中的系数随工况改变而改变，想利用解析求解是不可能的。另一方面实际问题的要求不一定需要严格的精确解，只要求达到一定精度，所以就可借助于差分方法来求偏微分方程的数值解。在第4章里，我们介绍了一个套管式换热器稳态的传热问题。如果我们考虑一个动态的传热过程，且不忽略纵向的热传导，就可以得到以下的偏微分方程：ltultCtTCrKtPWP22)(25.15.45.35.2总目录5.1问题的提出上面方程中变量的含义如下：通过求解上面的偏微分方程，就可以得到传热管各点温度随时间的变化值，从而确定达到传热平衡所需的时间，为实验测量提供依据。想求解上述方程，就必须首先学会偏微分方程的求解方法，下面我们首先介绍如何对偏微分方程进行离散化的工作，然后再对各类不同的偏微分方程进行求解，我们一般只给出离散化的基本公式及计算方法，对离散化公式的具体推导工作一般不作详细介绍，对这方面感兴趣的读者可自行参考有关数值计算的书籍。ssKmJmrmKgKKgJCsmuKTlKtPW，时间，，流体导热系数，内套管半径，套管内流体的密度套管内流体的比热套管内流体的速度套管的管壁温度置流体在套管内所处的位套管内某一点的温度/,/,/,,/,,,m,,,,3,5.15.45.35.2总目录5.2基本离散化公式在偏微分方程中，自变量都在两个或两个以上，应变量随两个或两个以上的自变量变化而变化。在化工或化学动态模拟方程中，常常有一个自变量是时间，其它的自变量为空间位置。如果只考虑一维空间，则只有两个自变量；如果考虑两维空间，则有3个自变量。一般我们将自变量在时间和空间以一定的间隔进行离散化，则应变量就变成了这些离散变量的函数，以3维空间为例，我们将离散化的应变量表示成，它所表示的真正含义如下:有了以上的定义，对于一阶偏导我们可以利用第四章的欧拉公式直接得出向前欧拉公式：对于时间偏导而言，有时我们常常采用向后欧拉公式，时间的向后欧拉公式如下：zkzyjyxixtntnkjizyxtuu,,,,,),,,(xuuzuyuuyuxuuxutuutunkjinkjizkzyjyxixtntnkjinkjizkzyjyxixtntnkjinkjizkzyjyxixtntnkjinkjizkzyjyxixtnt,,1,,,,,,,,1,,,,,,,,1,,,,,1,,,,,tuutunkjinkjizkzyjyxixtnt,,1,,,,,)1(5.15.45.35.2总目录5.2基本离散化公式这样在以后的计算中，得到的是隐式的计算公式，需通过求解线性方程组才能求解。具体的计算过程我们在下面会针对具体的偏微分方程进行讲解。对于二阶偏导，我们可以通过对泰勒展开式处理技术得到下面离散化计算公式：有了以上的离散化公式，就可以进行偏微分方程的数值求解工作。当然，在具体求解时，还会碰到不同的问题，需要区别对待，同时在利用计算机编程计算时也会碰到困难，这些问题我们会通过具体的例子加以说明。21,,,,1,,,,,222,,1,,,1,,,,222,,1,,,,1,,,221,,,,1,,,,,22)(2)(2)(22zuuuzuyuuuyuxuuuxutuuutunkjinkjinkjizkzyjyxixtntnkjinkjinkjizkzyjyxixtntnkjinkjinkjizkzyjyxixtntnkjinkjinkjizkzyjyxixtnt5.15.45.35.2总目录5.3几种常见方程的离散化计算1、波动方程其中：为初值条件为边值条件当该波动方程只提初值条件时，称此方程为波动方程的初值问题，二者均提时，称为波动方程的混合问题。对于初值问题，是已知t=0时，u与依赖于x的函数形式，求解不同位置，不同时刻的u值。而u是定义在的二元函数，即上半平面的函数。对于混合问题除初值外，还有边值。是已知初值及x=0及x=l时u依赖于t的函数，求解不同位置x，不同时刻的u值。此时u是定义在的带形区域上的二元函数。如图可以看出初值问题和混合问题的定义域。)(),()(),(),(2100022222tutuxtuxutxfxuatulxxtttuxt,0)(),(00xtuxutt)(),(210tututxxlxt0,05.15.45.35.2总目录5.3几种常见方程的离散化计算根据5.2节提供的公式，将上面波动方程离散化，得到：(5-1)将式（5-1）进行处理，把（n+1）时刻的变量留在右边，其余放在左边得到：（5-2）同时将边界条件和初始条件也离散化，得到：(5-3)xt0a)初值问题tx0lb)混合问题),2,1(1)-m,1,2,(i0)(2)(22112211nxuuuatuuunininininini1122222212221)()())()(22()()(nininininiuuxtauxtauxtau)1,2,(n)(),(m),1,2,(i)(),(210010tnutnuxituuxjunmniii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）可写为)0()(),0()0(21l),,2,1,0(0miui)0(t)2,1,0(,0nuunmn1iu)1,,2,1)((01iixituuii1iu2iu2iu3iu1iu)1,,2,1)((10mixituuii1iu1iu2iu0),(txf),()(2)(22112211tnxifxuuuatuuunininininini5.15.45.35.2总目录5.3几种常见方程的离散化计算当方程是初值问题时，边界条件没有了，由于在t=0时，u与值是已知的，若需要求某的值，只要按“波及原则”多算一些初值，即可推得，所图5-2所示。为了保证差分方程的解在时收敛于原来波动方程的解，要求式（5-2）中等式右边的各项系数均大于0，即：化简得：而且，可以证明，只要初始条件，边界条件满足一定的光滑性要求，且满足收敛关系式时，差分格式是稳定的。tuniu0,0tx0)()(22222xta1xta图5-2x(b)(a)0xJx5.15.45.35.2总目录5.3几种常见方程的离散化计算例5.1:用数值法求解下面偏微分方程，并写出VB程序。解：首先根据前面的知识，将所求的方程离散化，先假设以下各式：代入微分方程并化简得：（5-3）分析式（5-3）可知，如果知道了某一时刻的各点t，（j=0,1,2….100），就可以求下一时刻的各点温度值。有了以上各式，上面的微分方程就可以求解了。其实这个微分方程，是在不考虑流体本身热传导时的套管传热微分方程。由计算结果可知，当计算的时间序列进行到７２时，传热过程已达到稳态，各点上的温度已不随时间的增加而改变。如果改变套管长度或传热系数，则达到稳态的时间亦会改变。0,1030,30,1503)(2100xnjWWxtxttTxttTt1.0,01.011xxttxttttnjnjnjnjnjnjWnjttTt113.068.002.05.15.45.35.2总目录5.3几种常见方程的离散化计算2、一维流动传热传导方程的混合问题与波动方程的情形类似，用差商近似代替偏商，可以得到差分方程，以其解作为流动传热传导方程的近似解。一维流动传热传导方程的混合问题：上面的偏微分方程其实就是在5.1节中提出的偏微分方程，利用5.2节中的离散化公式进行离散化，得到其离散化公式：t)((t)μut)(xul)x((x),ut)l,x(tufxubxuatuxlxt000000),(100222)0,1,2,()(),2,1,0(0)21(i)(),()(21010121121ntnunxuu,m,,xiutnxifxuubxuuuatuunnmnmininininininini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）niniuxtbxtatnxiftu1221))((),(niniuxtauxtbxta12222)())(21(0iu001mmuu1iuxniut))1(,(tnxi))1(,()(211121111121tnxifxuubxuuuatuuninininininini5.15.45.35.2总目录5.3几种常见方程的离散化计算从图5-3中可见要由初值及边界条件一排一排推上去是不行的，需解线性方程组,同时添上二边界条件：正好共有m+2个方程，同时有m+2个