计算流体力学课程作业

计算流体力学大作业

——有限差分法解Poisson 方程

五点格式解区域内Poisson 方程

摘要：本文结合计算流体力学课上所学知识，采用数值解法中的有限差分法求解Poisson 方程（偏微分方程中椭圆型方程的一种），并用其五点格式采用高斯—塞德尔（Gauss-Seidel ）迭代求解。并比较了数值近似解与真实解，以及不同步长情况下误差的大小，得到了一定的结论。

关键词：Poisson 方程 有限差分法 五点格式

一、计算流体流体力学的特点

计算流体力学中许多问题求解最终都会变成偏微分方程的求解，而在数学上，除了几种极少数情况外，要求出它们精确解是很难的。计算机技术的发展使得这一难题的一很好地解决。

二、偏微分方程的种类

2.1、 椭圆型偏微分方程

椭圆型偏微分方程的一般形式为

(

)(,)div c u au f x t -?+= 其中：若12(,,,,)(,)n u u x x x t u x t ==，u ?为u 的梯度，则其定义为 12

,,,n u u x x x ??????=??????? 散度()div v 的定义为

12

()n div v v x x x ?????=+++ ??????

这样，()div c u ?可以更明确地表示为

1122()n n u u u div c u c c c x x x x x x ???????????????=++

+?? ? ? ???????????

????

若c 为常数，则进一步化简为 22

222212()n div c u c u c u x x x ??????=+++=? ??????

其中，?又称为Laplace 算子。这样椭圆型偏微分方程可以简单地写为

22

222212

(,)n c u au f x t x x x ?????-++++= ??????

2.2、抛物型偏微分方程

抛物型偏微分方程的一般形式为 ()(,)u d div c u au f x t t

?-?+=? 根据上面叙述，若c 为常数，则该方程可以更简单地写为

22222212

(,)n u d c u au f x t t x x x ??????-++++= ??????? 2.3、双曲型偏微分方程

双曲型偏微分方程的一般形式为

22()(,)u d div c u au f x t t

?-?+=? 若c 为常数，则可以将该方程简化为

2222222212(,)n u d c u au f x t t x x x ??????-++++= ???????

三类方程的直接的区别在于u 对t 的导数的阶次。

若对t 没有求导，可以理解为其值为常数，故称为椭圆型的。

若取u 对时间t 的一阶导数，则与u 对x 的二阶导数直接构成了抛物线关系，故称为抛物型偏微分方程。

若取u 对时间t 的二阶导数，称其为双曲型偏微分方程。

三、Poisson 方程：

泊松方程是数学中一个常见于静电学、机械工程和理论物理的偏微分方程。是从法国数学家、几何学家及物理学家泊松而得名的。

泊松方程一般可写为：

△φ=f

在这里△代表的是拉普拉斯算符(也就是哈密顿算符▽的平方)，而 f 和φ可以是在流形上的实数或复数值的方程。

拉普拉斯方程：

因此泊松方程通常写成：

在三维直角坐标系，可以写成

如果没有f，这个方程就会变成拉普拉斯方程△φ=0.

泊松方程可以用格林函数来求解；如何利用格林函数来解泊松方程可以参考screened Poisson equation。现在有很多种数值解法，数学上，泊松方程属于椭圆型方程。

四、Poisson方程的解法

泊松首先在无引力源的情况下得到泊松方程，△Φ=0（即拉普拉斯方程）；当考虑引力场时，有△Φ=f（f为引力场的质量分布）。后推广至电场磁场，以及热场分布。该方程通常用格林函数法求解，也可以分离变量法，特征线法求解。由于Poisson方程难以求得其解析解，计算机技术发展之后，数值解法成为工程实际中应用最广泛的求Poisson方程解的方法，常见的数值解法有：里兹（Ritz）法，加权余量法，有限差分法，有限元法，边界元法及有限体积法等。

4.1、里兹（Ritz）法

瑞利-里兹法（也称里兹法）是通过泛函驻值条件求未知函数的一种近似方法，是英国的瑞利于1877年在《声学理论》一书中首先采用，后由瑞士的W.

里兹于1908年作为一个有效方法提出。这一方法在许多力学、物理学、量子化学问题中得到应用。

同时它也是广泛应用于应用数学和机械工程领域的经典数值方法，它可以用来计算结构的低阶自然频率。它是直接变分法的一种，以最小势能原理为理论基础。通过选择一个试函数来逼近问题的精确解，将试函数代入某个科学问题的泛函中，然后对泛函求驻值，以确定试函数中的待定参数，从而获得问题的近似解。

4.2、加权余量法

加权余量法（Weighted residual approach），又称加权残量法，加权残余法。

当n 有限时，定解方程存在偏差（余量）。取权函数，强迫余量在某种平均意义上为零。采用使余量的加权积分为零的等效积分的“弱”形式来求得微分方程近似解的方法称为加权余量法。

加权余量法在固体力学中，是求解线性、非线性微分方程的一种有效方法，它是基于等效积分形式的近似方法，也是通用的数值计算方法．有限元法、边界元法、无网格法都是加权余量法的特殊情况，由于这三种方法各有其特点，所以都各自发展为一种独立的方法，加权余量法最早是用于流体力学，传热等科学领域，后在固体力学中得到了更大的发展。

权函数的选择

加权余量法是求解微分方程近似解的一种有效方法．显然，任何独立的完全函数都可用来作为权函数，加权余量法可分为内部法、边界法和混合法，在内部法中，又可分为：配置法（以笛拉克函数δ作为权函数），子域法，最小二乘法，矩量法，伽辽金法等。

4.3、有限差分法

有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来

代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。

在采用数值计算方法求解偏微分方程时，若将每一处导数由有限差分近似公式替代，从而把求解偏微分方程的问题转换成求解代数方程的问题，即所谓的有限差分法。

有限差分法求解偏微分方程的步骤如下：

1、区域离散化，即把所给偏微分方程的求解区域细分成由有限个格点组成的网格；

2、近似替代，即采用有限差分公式替代每一个格点的导数；

3、逼近求解。换而言之，这一过程可以看作是用一个插值多项式及其微分来代替偏微分方程的解的过程。

4.4、有限元法

有限元法（finite element method）是一种高效能、常用的计算方法。有限

元法在早期是以变分原理为基础发展起来的，所以它广泛地应用于以拉普拉斯方程和泊松方程所描述的各类物理场中（这类场与泛函的极值问题有着紧密的联系）。自从1969年以来，某些学者在流体力学中应用加权余数法中的迦辽金法(Galerkin)或最小二乘法等同样获得了有限元方程，因而有限元法可应用于以任何微分方程所描述的各类物理场中，而不再要求这类物理场和泛函的极值问题有所联系。基本思想：由解给定的泊松方程化为求解泛函的极值问题。

将连续的求解域离散为一组单元的组合体，用在每个单元内假设的近似函数来分片的表示求解域上待求的未知场函数，近似函数通常由未知场函数及其导数在单

元各节点的数值插值函数来表达。从而使一个连续的无限自由度问题变成离散的有限自由度问题。

有限元法的一般步骤：

步骤1：剖分:

将待解区域进行分割，离散成有限个元素的集合．元素（单元）的形状原则上是任意的．二维问题一般采用三角形单元或矩形单元，三维空间可采用四面体或多面体等．每个单元的顶点称为节点（或结点）．

步骤2：单元分析:

进行分片插值，即将分割单元中任意点的未知函数用该分割单元中形状函数及离散网格点上的函数值展开，即建立一个线性插值函数

步骤3：求解近似变分方程

用有限个单元将连续体离散化，通过对有限个单元作分片插值求解各种力学、物理问题的一种数值方法。有限元法把连续体离散成有限个单元：杆系结构的单元是每一个杆件；连续体的单元是各种形状（如三角形、四边形、六面体等）的单元体。每个单元的场函数是只包含有限个待定节点参量的简单场函数，这些单元场函数的集合就能近似代表整个连续体的场函数。根据能量方程或加权残量方程可建立有限个待定参量的代数方程组，求解此离散方程组就得到有限元法的数值解。有限元法已被用于求解线性和非线性问题，并建立了各种有限元模型，如协调、不协调、混合、杂交、拟协调元等。有限元法十分有效、通用性强、应用广泛，已有许多大型或专用程序系统供工程设计使用。结合计算机辅助设计技术，有限元法也被用于计算机辅助制造中。

4.5、边界元法

边界元法（boundary element method）是一种继有限元法之后发展起来的一种新数值方法，与有限元法在连续体域内划分单元的基本思想不同，边界元法是只在定义域的边界上划分单元，用满足控制方程的函数去逼近边界条件。所以边界元法与有限元相比，具有单元个数少，数据准备简单等优点。但用边界元法解非线性问题时，遇到同非线性项相对应的区域积分，这种积分在奇异点附近有强烈的奇异性，使求解遇到困难。

边界元法以定义在边界上的边界积分方程为控制方程，通过对边界分元插值离散，化为代数方程组求解。它与基于偏微分方程的区域解法相比，由于降低了问题的维数，而显著降低了自由度数，边界的离散也比区域的离散方便得多，可用较简单的单元准确地模拟边界形状，最终得到阶数较低的线性代数方程组。又由于它利用微分算子的解析的基本解作为边界积分方程的核函数，而具有解析与数值相结合的特点，通常具有较高的精度。特别是对于边界变量变化梯度较大的问题，如应力集中问题，或边界变量出现奇异性的裂纹问题，边界元法被公认为比有限元法更加精确高效。由于边界元法所利用的微分算子基本解能自动满足无限远处的条件，因而边界元法特别便于处理无限域以及半无限域问题。边界元法的主要缺点是它的应用范围以存在相应微分算子的基本解为前提，对于非均匀介质等问题难以应用，故其适用范围远不如有限元法广泛，而且通常由它建立的求解代数方程组的系数阵是非对称满阵，对解题规模产生较大限制。对一般的非线性问题，由于在方程中会出现域内积分项，从而部分抵消了边界元法只要离散边界的优点。

4.6、有限体积法

有限体积法（Finite Volume Method）又称为有限容积法、控制体积法。

将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待解的微分方程对每一个控制体积积分，便得出一组离散方程。其中的未知数是网格点上的因变量的数值。为了求出控制体积的积分，必须假定值在网格点之间的变化规律，即假设值的分段的分布的分布剖面。

从积分区域的选取方法看来

有限体积法属于加权剩余法中的子区域法；从未知解的近似方法看来，有限体积法属于采用局部近似的离散方法。简言之，子区域法属于有限体积法的基本方法。

有限体积法的基本思路易于理解，并能得出直接的物理解释。离散方程的物理意义，就是因变量在有限大小的控制体积中的守恒原理，如同微分方程表示因变量在无限小的控制体积中的守恒原理一样。有限体积法得出的离散方程，要求因变量的积分守恒对任意一组控制体积都得到满足，对整个计算区域，自然也得到满足。这是有限体积法吸引人的优点。有一些离散方法，例如有限差分法，仅当网格极其细密时，离散方程才满足积分守恒；而有限体积法即使在粗网格情况下，也显示出准确的积分守恒。

就离散方法而言有限体积法可视作有限单元法和有限差分法的中间物。有限单元法必须假定值在网格点之间的变化规律（既插值函数），并将其作为近似解。有限差分法只考虑网格点上的数值而不考虑值在网格点之间如何变化。有限体积法只寻求的结点值，这与有限差分法相类似；但有限体积法在寻求控制体积的积分时，必须假定值在网格点之间的分布，这又与有限单元法相类似。在有限体积法中，插值函数只用于计算控制体积的积分，得出离散方程之后，便可忘掉插值函数；如果需要的话，可以对微分方程中不同的项采取不同的插值函数。

五、有限差分法详细介绍

在计算流体力学中，由于有限差分法发展较早，应用技术比较成熟，其应用比较的广泛。

物理学和其他学科领域的许多问题在被分析研究之后, 往往可以归结为常微分方程或偏微分方程的求解问题。一般说来，处理一个特定的物理问题，除了需要知道它满足的数学方程外，还应当同时知道这个问题的定解条件，然后才能设计出行之有效的计算方法来求解。

有限差分法以变量离散取值后对应的函数值来近似微分方程中独立变量的连续取值。在有限差分方法中，我们放弃了微分方程中独立变量可以取连续值的特征，而关注独立变量离

散取值后对应的函数值。但是从原则上说，这种方法仍然可以达到任意满意的计算精度。因为方程的连续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值插值计算来近似得到。

这种方法是随着计算机的诞生和应用而发展起来的。其计算格式和程序的设计都比较直观和简单，因而，它的实际应用已经构成了计算数学和计算物理的重要组成部分。

有限差分法的具体操作分为两个部分：

(1)用差分代替微分方程中的微分，将连续变化的变量离散化，从而得到差分方程组的数学形式；

(2)求解差分方程组。在第一步中，我们通过所谓的网络分割法，将函数定义域分成大量相邻而不重合的子区域。通常采用的是规则的分割方式。这样可以便于计算机自动实现和减少计算的复杂性。网络线划分的交点称为节点。若与某个节点P 相邻的节点都是定义在场域内的节点，则P 点称为正则节点；反之，若节点P 有处在定义域外的相邻节点，则P 点称为非正则节点。 在第二步中，数值求解的关键就是要应用适当的计算方法，求得特定问题在所有这些节点上的离散近似值。

有限差分法的差分格式

一个函数在x 点上的一阶和二阶微商，可以近似地用它所临近的两点上的函数值的差分来表示。如对一个单变量函数f(x)，x 为定义在区间[a,b]的连续变量。以步长h=Δx 将[a,b]区间离散化，我们得到一系列节点公式

()()()j i j i j i u y x u j i y x jl y ih x l y h x ,,，,,，

，，，??===?=? 则可用差商代替微商 前插和后插格式：()()j i j i j

i j i j i j i u u h x u u u h x u ,1,,,,1,11-+-=???? ????-=???? ???? 中心差分格式： ()j i j i j

i u u h x u ,1,1,21-+-=???? ????

下面即用五点格式求解Poisson 方程的Matlab 程序，一类边界条件（即在边界上给定边界函数），采用高斯-赛德尔（Gauss-Seidel ）迭代，求解。若P 点不恰好在边界上（如下图所示），则用直接转移法定义。在x 区域上给定函数为22)y ,(y x x f +=，

则其精确解即为22)y ,(y x x f +=，边界条件亦为22)y ,(y x x f +=。用五点格式求其精确解，并输出其解析解及其与精确解的误差。

P点不恰好在边界上的节点处理

六、Matlab程序及数值解分析

% Uses the five-point formula to solve Poisson's equation

% on a square.

% Inline functions define the source function (f), the boundary

% values (ux0,ux1,u0y,u1y), and the true solution (u).

f = inline('x^2 + y^2');

ux0 = inline('0');

ux1 = inline('.5*x^2');

u0y = inline('sin(pi*y)');

u1y = inline('exp(pi)*sin(pi*y) + .5*y^2');

u = inline('exp(pi*x)*sin(pi*y) + .5*(x^2)*(y^2)');

% Ask user for problem size. (Use same no. of points in each direction.)

n = input(' Enter number of subintervals in each direction: ');

h = 1/n;

N = (n-1)^2;

A = sparse(N,N); % Store A as a sparse matrix.

F = zeros(N,1);

% Set up matrix A and right-hand side vector F.

A = -4*sparse(eye(N,N));

for j=1:n-1, % Loop over grid rows

for i=1:n-1, % Loop over points in each row

k = (j-1)*(n-1)+i; % Index of this point

if j > 1,

A(k,k-(n-1)) = 1; % coupling to pt below

end;

if j < n-1,

A(k,k+(n-1)) = 1; % coupling to pt above end;

if i > 1, A(k,k-1) = 1; end; % coupling to pt to left

if i < n-1, A(k,k+1) = 1; end; % coupling to pt to right

xi = i*h; yj = j*h;

F(k) = f(xi,yj); % right-hand side vector if j==1,

F(k) = F(k) - (1/h^2)*ux0(xi); % bdry pt below

end;

if j==n-1,

F(k) = F(k) - (1/h^2)*ux1(xi); % bdry pt above

end;

if i==1,

F(k) = F(k) - (1/h^2)*u0y(yj); % bdry pt to left

end;

if i==n-1,

F(k) = F(k) - (1/h^2)*u1y(yj); % bdry pt to right

end;

end;

end;

% Remember to multiply A by 1/h^2.

A = (1/h^2)*A;

% Solve linear system.

uapprox = A\F;

% Compare to true solution.

utrue = zeros(N,1);

for j=1:n-1,

for i=1:n-1,

k = (j-1)*(n-1)+i;

xi = i*h; yj = j*h;

utrue(k) = u(xi,yj);

end;

end;

err2 = h*norm(utrue-uapprox),

errinf = norm(utrue-uapprox,'inf')

% Plot solution and error over grid.

ugrid = reshape(utrue,n-1,n-1);

figure(1)

mesh([h:h:(n-1)*h]',[h:h:(n-1)*h]',ugrid')

title('True Solution')

errgrid = reshape(utrue-uapprox,n-1,n-1);

figure(2)

mesh([h:h:(n-1)*h]',[h:h:(n-1)*h]',errgrid')

title('Error')

运行，输入两个方向的划分节点数（程序中默认两者相等）10，

Enter number of subintervals in each direction: 10，

得到结果：误差矩阵的范数为err2 =0.0329

真实解图误差图再次运行，输入两个方向的划分节点数（程序中默认两者相等）100，Enter number of subintervals in each direction: 100，

得到结果：误差矩阵的范数为err2 =

3.3207e-004

真实解图误差图

对比两次运行结果可以发现，随着节点数的增多，数值解的误差变小，同时求解速度也变大，实际工程中应按照工程类型合理选择节点数，以提高效率。

参考文献：

《计算流体力学》李万平华中科技大学出版社 2004-10

《Matlab教程》张志涌北京航空航天大学出版社 2006年8月