第5章偏微分方程数值解

1 i 0 i 0 i −1 i
n n n u in +1 − 2u in + u in −1 2 u i +1 − 2u i + u i −1 −a = f (i∆x, n∆t ) (∆t ) 2 (∆x) 2
5.1
5.2
5.3
5.4
总目录
5.3几种常见方程的离散化计算 5.3几种常见方程的离散化计算
5.1
5.2
5.3
5.4
总目录
5.2 基本离散化公式
这样在以后的计算中，得到的是隐式的计算公式，需通过求解线性方 程组才能求解。具体的计算过程我们在下面会针对具体的偏微分方程进行 uin, +1k − 2uin, j ,k + uin, −1k 讲解。对于二阶偏导，我们可 ∂ 2 u j, j, = ∆t 以通过对泰勒展开式处理技术 ∂t 2 t =n∆t , x=i∆x, y = j∆y , z =k∆z 得到下面离散化计算公式： ∂ 2 u un − 2u n + u n
5.2
5.3
5.4
总目录
5.3几种常见方程的离散化计算 5.3几种常见方程的离散化计算
这样，由式（5-2），并结合式（5-3），就可以从n时刻的各点u值，计 算得到下一时刻的u值，这样层层递推，就可以计算出任意时刻，任意位 置的u值。而图5.2则表明了这种层层递推的计算过程，在图5.2中*表示 需求u值的点，○表示为了求x点的u 值必须已知u值的点。 需要说明的是，在应用式（5-2）进行计算时，初值与边值应当满足相容 0 ( 0) ) 性条件 ϕ (0) = µ1 (0), ϕ (l ) = µ 2 。由初值得到 u i (i = 0,1,2,⋯ , m，由边值得到(t ≥ 0) ， n n u 0 , u m ( n = 0,1,2 ⋯), 但在利用式（5-2）进行第一轮计算时，若取n=0，则发 1 现等式右边出现了 u i−，这是一个无法计算的值。这时可以利用另一个初 1 u −u 值条件 ∆t = ψ (i∆x)(i = 1,2,⋯, i − 1) 算得 u i ，这样，可在第一轮计算的时候，取 n=1，计算得到 u i2 ，由 u i2，递推得到 u i3 ，这样就可由式（5-2）一排一 排往上推，计算得到所有希望得到的u值。对于式（5-2）取n=0计算中碰 到的 u i−1 ，也可利用另一种方法进行计算，解决的办法是将另一个初值条 u −u −1 件利用向后欧拉离散化 ∆t = ψ (i∆x)(i = 1,2,⋯, m − 1) 算得 u i ，这样利用式（51 2 2），取n=0就可以得到 u i ，取n=1，u i ，和前一种处理方法一样一排一 排往上推，计算得到所有希望得到的u值。象这样可以用已知点上函数值 直接推出所有点上函数值的格式，称为显式格式。当方程非齐次 时，f ( x, t ) ≠ 0 ，式（5-1）可写为
∂x
2 t = n∆t , x =i∆x , y = j∆y , z = k∆z
=
i +1, j , k
i , j ,k
i −1, j ,k
(∆x)
2
uin, j +1,k − 2uin, j ,k + uin−1, j ,k ∂ 2u = (∆y) 2 ∂y 2 t =n∆t , x=i∆x , y = j∆y , z =k∆z uin, j ,k +1 − 2uin, j ,k + uin, j ,k −1 ∂ 2u = 2 ∂z t =n∆t , x=i∆x , y = j∆y , z =k∆z (∆z ) 2
5.4
总目录
5.3几种常见方程的离散化计算 5.3几种常见方程的离散化计算
1、 波动方程 、
∂u ∂t
其中：
u
t =0
= ϕ ( x ),
t =0
为初值条件
u
x =0
= µ1 (t ), u x =t = µ 2 (t ) 为边值条件
当该波动方程只提初值条件时，称此方程为波动方程的初值问题，二者均 提时，称为波动方程的混合问题。
(∆t ) (∆x) (n = 1,2,⋯)
将式（5-1）进行处理，把（n+1）时刻的变量留在右边，其余放 在左边得到： 2 2 2 n +1 2 ( ∆t ) n 2 ( ∆t ) n 2 ( ∆t ) ui = a u i +1 + ( 2 − 2a u in−1 + u in −1 )u i + a （5-2） 2 2 2
(∆x)
∂u ∂t
∆x
而且，可以证明，只要初始条件，边界条件满足一定的光滑性 要求，且满足收敛关系式时，差分格式是稳定的。
0 ∆x (a)
J ∆x
图5-2
x
(b)
5.1
5.2
5.3
5.4
总目录
5.3几种常见方程的离散化计算 5.3几种常见方程的离散化计算
例5.1: 用数值法求解下面偏微分方程，并写出VB程序。 ∂t ∂t
(∆x) (∆x) ( ∆x)
同时将边界条件和初始条件也离散化，得到：
u i1 − u i0 u = ϕ ( j∆x), = ψ (i∆x) ∆t n n u 0 = µ1 (n∆t ), u m = µ 2 (n∆t )
0 i
(i = 1,2,⋯ , m) (n = 1,2,⋯)
(5-3)
5.1
∂t ∂t 2K λ ∂ t = (TW − t ) + −u ∂τ rρC P ρC P ∂l 2 ∂l
5.1
5.2
5.3
5.4
总目录
5.1 问题的提出
t , 套管内某一点的温度 , K
上面方程中变量的含义如下：
l , 流体在套管内所处的位 置 , m TW , 套管的管壁温度 , K u , 套管内流体的速度 , m / s C P , 套管内流体的比热 , J / Kg ⋅ K
ρ 套管内流体的密度 , Kg / m 3
r ,内套管半径， m
λ ，流体导热系数， J / m ⋅ K ⋅ s τ ，时间， s
通过求解上面的偏微分方程，就可以得到传热管各点温度随时间 的变化值，从而确定达到传热平衡所需的时间，为实验测量提供依据。 想求解上述方程，就必须首先学会偏微分方程的求解方法，下面我们 首先介绍如何对偏微分方程进行离散化的工作，然后再对各类不同的 偏微分方程进行求解，我们一般只给出离散化的基本公式及计算方法， 对离散化公式的具体推导工作一般不作详细介绍，对这方面感兴趣的 读者可自行参考有关数值计算的书籍。
= 2(TW − t ) − 3 ∂τ ∂x 0 n TW = 150, t j = 30, t 0 = 30 0 ≤ x ≤ 1, ∂t ∂x =0
x =1
解：首先根据前面的知识，将所求的方程离散化，先假设以下各 式：∂t t nj +1 − t nj
∂τ = ∂t 代入微分方 = n n ∂x ∆x t n +1 程并化简得：j = 0.02TW + 0.68t j + 0.3t j −1（5-3） ∆τ = 0.01, ∆x = 0.1
5.1
5.2
5.3
5.4
总目录
5.2 基本离散化公式
在偏微分方程中，自变量都在两个或两个以上，应变量随两个或两 个以上的自变量变化而变化。在化工或化学动态模拟方程中，常常有一个 自变量是时间，其它的自变量为空间位置。如果只考虑一维空间，则只有 两个自变量；如果考虑两维空间，则有3个自变量。一般我们将自变量在 时间和空间以一定的间隔进行离散化，则应变量就变成了这些离散变量的 函数，以3维空间为例，我们将离散化的应变量表示成，它所表示的真正 n 含义如下 : u i , j , k = u ( t , x , y , z ) t = n ∆ t , x = i ∆ x , y = j ∆ y , z = k ∆ z 有了以上的定义，对于一阶偏 u in, +1k − u in, j ,k ∂u j, = 导我们可以利用第四章的欧拉 ∂ t t = n ∆ t , x = i ∆ x , y = j ∆ y , z = k∆ z ∆t 公式直接得出向前欧拉公式： u in+1, j , k − u in, j ,k ∂u 对于时间偏导而言，有时我们 = ∂ x t = n ∆ t , x = i ∆ x , y = j ∆ y , z = k∆ z ∆x 常常采用向后欧拉公式，时间的 u in, j +1, k − u in, j ,k ∂u 向后欧拉公式如下： = ∂ y t = n ∆ t , x = i ∆ x , y = j ∆ y , z = k∆ z ∆y u in, +,1k − u in, j ,k ∂u j = u in, j ,k +1 − u in, j ,k ∂u ∂t t =( n +1) ∆t , x =i∆x , y = j∆y , z = k∆z ∆t = ∂ z t = n ∆ t , x = i ∆ x , y = j ∆ y , z = k∆ z ∆x
有了以上的离散化公式，就可以进行偏微分方程的数值求解工作。当 然，在具体求解时，还会碰到不同的问题，需要区别对待，同时在利用计 算机编程计算时也会碰到困难，这些问题我们会通过具体的例子加以说明。
5.1
5.2
5.3
∂ 2u ∂ 2u − a 2 2 = f ( x, t ) 2 ∂x ∂t ∂u = ψ ( x) u t =0 = ϕ ( x), ∂t t =0 u x =0 = µ1 (t ), u x =l = µ 2 (t ) = ψ ( x)
5.1
5.2
5.3
t
5.4
总目录
5.3几种常见方程的离散化计算 5.3几种常见方程的离散化计算