某二階偏微分方程計(jì)算方法對(duì)比分析

賈濤 牛華東 解景濤 盧瑞軍 蘇茂輝

摘要：CAE在汽車和發(fā)動(dòng)機(jī)設(shè)計(jì)中發(fā)揮越來越重要的作用，然而CAE的本質(zhì)是求解各種常微分方程組或偏微分方程組;常微分方程又稱動(dòng)力系統(tǒng)方程，通常用來求解動(dòng)力系統(tǒng)問題。偏微分方程則廣泛應(yīng)用于電磁學(xué)、流體力學(xué)、結(jié)構(gòu)力學(xué)等多個(gè)領(lǐng)域;對(duì)于偏微分方程通常很難求得其解析解，需要借助數(shù)值計(jì)算來獲取其方程的近似解。拉普拉斯方程（Laplace）是一種橢圓形二階偏微分方程，并且可以求取其解析解。本文通過解析法以及數(shù)值法對(duì)拉普拉斯方程求解，并對(duì)比不同求解方法的效率和精度;結(jié)論顯示解析法雖然精度較高，但是需要很大的計(jì)算量，并且大多數(shù)偏微分方程沒有解析解。因此，在汽車和發(fā)動(dòng)機(jī)等工程應(yīng)用中應(yīng)該根據(jù)精度需求選擇最優(yōu)的途徑求解偏微分方程問題。

關(guān)鍵詞：Laplace方程;偏微分方程;數(shù)值計(jì)算;CAE

1? 理論分析

拉普拉斯方程的數(shù)學(xué)表達(dá)式如下，

2? 問題描述

研究某矩形二維傳熱問題，已知傳熱問題的導(dǎo)熱微分方程式如下[2]：

3? 解析解

求解偏微分方程的核心在于邊界條件，不同的邊界條件甚至?xí)霈F(xiàn)完全不同的解;根據(jù)該問題描述，通過邊界條件將上述方程簡(jiǎn)化：

由于公式（5）比較繁瑣，并且涉及到n（0～∞）求和計(jì)算;方程組規(guī)模小則計(jì)算量小，可以快速得到計(jì)算結(jié)果;方程組規(guī)模大則計(jì)算量很大，通常需要很久才能完成計(jì)算。本例中算法1取n=110，對(duì)5000個(gè)數(shù)據(jù)編程計(jì)算，則耗時(shí)近1小時(shí);而采用算法2則只需要數(shù)秒時(shí)間。

算法1：

clc;

clear;

sum=zeros（100，50）;

syms x;

for i=1：100;

for j=1：50;

for k=1：110;

Dn=int（100*sin（k*pi*x/50），0，50）;

An=1/（25*sinh（2*pi*k））;

Cn=sinh（k*pi*（100-i）/50）;

Bn=sin（k*pi*j/50）;

Un=Dn*An*Cn*Bn;

sum（i，j）=sum（i，j）+Un;

end

end

end;

算法2：

clc;

clear;

sum=zeros（100，50）;

for i=1：100;

for j=1：50;

for k=1：110;

Dn=5000*（1-cos（k*pi））/（k*pi）;

An=1/（25*sinh（2*pi*k））;

Cn=sinh（k*pi*（100-i）/50）;

Bn=sin（k*pi*j/50）;

Un=Dn*An*Cn*Bn;

sum（i，j）=sum（i，j）+Un;

end

end;

end;

4? 數(shù)值解法

通過數(shù)值方法求解偏微分方程是常用的求解方法，如下通過有限元方法以及EXCEL迭代、軟件編程等不同方法進(jìn)行數(shù)值計(jì)算，對(duì)比不同計(jì)算方法的精度和計(jì)算效率。

4.1 有限元方法

借助有限元軟件，建立平面模型并搭建穩(wěn)態(tài)熱力學(xué)模型，生成5151個(gè)節(jié)點(diǎn)、10000個(gè)單元;設(shè)置邊界條件如上所述，計(jì)算得到的溫度場(chǎng)分布如圖2所示：該方法更偏向于工程的應(yīng)用。在求解過程中，模型搭建過程花費(fèi)較多時(shí)間;但是求解計(jì)算效率較高，普通電腦計(jì)算只需要20sCPU時(shí)間。

4.2 Excel迭代方法

設(shè)置步長(zhǎng)dx=dy=1，采用五點(diǎn)差分法求解;如圖4所示，使用EXCEL建立數(shù)值計(jì)算表格，啟用迭代計(jì)算，設(shè)置誤差以及最大迭代次數(shù)進(jìn)行求解;該方法計(jì)算效率較高，只需要1～2s CPU時(shí)間;得到計(jì)算結(jié)果如圖5[3]。

4.3 軟件編程計(jì)算

借助軟件編寫五點(diǎn)差分程序，使用高斯賽德爾迭代法，同樣可以計(jì)算出最終的結(jié)果;編程耗時(shí)較長(zhǎng)，計(jì)算耗時(shí)60s左右的CPU時(shí)間[4]。

計(jì)算程序：

clear

clc

u=1：1：50;

v=1：1：50;

[x，y]=meshgrid（u，v）;%通過坐標(biāo)軸點(diǎn)在平面畫網(wǎng)格

p=[x，y];

plot（p）;

spy（p）;

for i=2：50

for j=2：100

p（i，j）=（j-1）+98*（i-2）;

end;

end;

for i=1;

for j=1：100

p（i，j）=0;

end;

end;

for i=50;

for j=1：100

p（i，j）=0;

end;

end;

for j=1;

for i=1：50

p（i，j）=0;

end;

end;

for j=100;

for i=1：50

p（i，j）=0;

end;

end;

PP=delsq（p）;

load（'my.dat'）;

b=spconvert（my）;

N=length（b）;? ?%Gauss – Seidel迭代法

u=inv（PP）*b;

u=zeros（N，1）;%迭代初始值

D=diag（diag（PP））;

E=-tril（PP，-1）;%下三角

F=-triu（PP，1）;%上三角

B=inv（D-E）*F;g=inv（D-E）*b;

eps=0.000001;

for k=1：10000 %最大迭代次數(shù)

y=B*u+g;

if norm（u-y）<eps

break;

end

u=y;

end

s=zeros（50，100）;

for i=1：48

for j=1：98;

s（i+1，j+1）=u（98*（i-1）+j，1）;

end

end

for i=1：50

s（i，1）=100;

end

subplot（1，1，1）;

mesh（s）;

5? 誤差分析

對(duì)比計(jì)算結(jié)果第二列（x=1，u（x，y））數(shù)據(jù)并繪制圖表;如圖8所示，數(shù)值解和解析解的四條曲線重合在一起，說明這些數(shù)值大小相同，整體趨勢(shì)一致。圖9誤差分析顯示，EXCEL和編程計(jì)算誤差重合在一起，誤差在0～0.01%之間;有限元計(jì)算在中間部分誤差接近0，說明其計(jì)算誤差相對(duì)較小;三種數(shù)值計(jì)算方法誤差在工程可以接受的范圍內(nèi)。

6? 總結(jié)

本文通過不同方法求解偏微分方程，對(duì)偏微分方程的求解思路有了深入的理解;各種解法中，效率最高的是使用EXCEL迭代求解，操作簡(jiǎn)單并且計(jì)算速度快、精度較高。借助軟件編程的方法求解效率相對(duì)較低，并且需要對(duì)程序調(diào)試，精度方面和EXCEL求解差別不大。借助有限元求解則需要搭建有限元模型，并設(shè)置邊界條件，計(jì)算效率同樣很高，精度方面可以滿足工程需要。而解析法求解拉普拉斯方程雖然精度高，但是對(duì)算法要求較高。因此，對(duì)于偏微分方程工程方面應(yīng)用，在滿足精度的前提下，解析解不一定是最優(yōu)的方法;具體需要根據(jù)實(shí)際情況選擇相應(yīng)的求解方法。

參考文獻(xiàn)：

[1]陸金甫，關(guān)治.偏微分方程數(shù)值解法[M].清華大學(xué)出版社，2004.

[2]楊世銘，陶文銓.傳熱學(xué)[M].四版.高等教育出版社，2006.

[3]賈新民，嚴(yán)文.有限差分法求解拉普拉斯方程[J].昌吉學(xué)院學(xué)報(bào)，2009.

[4]劉大衛(wèi)，高明.正方形環(huán)域Laplace方程的簡(jiǎn)明數(shù)值解法[J].沈陽師范大學(xué)學(xué)報(bào)，2006.