1、哈爾濱工程大學航天與建筑工程學院王濱生 E-mail:wang_，,計算力學Computational Mechanics(13),本次課的主要問題： 差分原理 應力函數(shù)的差分解,基本內(nèi)容,差分原理 設(shè)有x的解析函數(shù)y=f(x)，從微分學知道函數(shù)y對x的導數(shù)為 式中 dy/dx是函數(shù)對自變量的導數(shù)，又稱微商 Dy、Dx分別稱為函數(shù)及自變量的差分 Dy/Dx為函數(shù)對自變量的差商。,差分原理,常微分方程數(shù)值解 工程上很多數(shù)學表述都可以歸結(jié)為常微分的定解問題，很多偏微分方程，也可以化為常微分方程問題來近似求解。 對于一個常微分方程 通常會有無窮個解，如 因此，要加上一個定解條件，如,差分原理,常微分

2、方程的解是一個函數(shù)，但是，計算機沒有辦法對函數(shù)進行運算，因此，常微分方程的數(shù)值解并不是求函數(shù)的近似值，而是求解函數(shù)在某些節(jié)點上的近似值。 這種方法，稱為數(shù)值離散方法。求的是未知函數(shù)在一系列離散點上的值的近似。 得到高精度方法的一個直接想法是利用Taylor展開。 假設(shè)式 中的 f(x,y) 充分光滑，將y(xi+1)在xi點作Taylor展開，若取右端不同的有限項作為y(xi+1)的近似值，就可得到計算y(xi+1)的各種不同截斷誤差的數(shù)值公式。,差分原理,例如：取前兩項可得到截斷誤差為O(x)2的公式 若取前三項，可得到截斷誤差為O(x)3的公式 類似地，若取前p+1項作為y(xi+1)的近

3、似值，便得到,差分原理,歐拉(Euler)法 尋求通過某種方法去獲得解y(x)在點xi上的近似值yi，即 為實現(xiàn)這一目標，Euler方法首先將微分算子離散化 即 推廣 這就是顯式的Euler公式，又稱差分格式。 顯然與y(xi+1)的Taylor展開式的前兩項完全相同，即局部截斷誤差為O(x)2 。 Euler公式可改寫成,差分原理,改進的歐拉法 對方程y = f(x,y)，兩邊由xi到xi+1積分，并利用梯形公式，有： 故有公式： 此即改進的歐拉法。,差分原理,同理，改進Euler公式可改寫成 局部截斷誤差為O(x)3 。 上述兩組公式在形式上共同點：都是用f(x,y)在某些點上值的線性組合

4、得出y(xi+1)的近似值yi+1，且增加計算的次數(shù)，可提高截斷誤差的階。 歐拉法：每步計算一次f(x,y)的值，為一階方法。 改進歐拉法：需計算兩次f(x,y)的值，為二階方法。,差分原理,龍格-庫塔(Runge-Kutta)法 于是可考慮用函數(shù)f(x,y)在若干點上的函數(shù)值的線性組合來構(gòu)造近似公式，構(gòu)造時要求近似公式在(xi,yi)處的Taylor展開式與解y(x)在xi處的Taylor展開式的前面幾項重合，從而使近似公式達到所需要的階數(shù)。既避免求高階導數(shù)，又提高了計算方法精度的階數(shù)。 或者說，在xi,xi+1這一步內(nèi)多計算幾個點的斜率值，然后將其進行加權(quán)平均作為平均斜率，則可構(gòu)造出更高精

5、度的計算格式。 這就是龍格-庫塔(Runge-Kutta)法的基本思想。,差分原理,一般龍格-庫塔方法的形式為 稱為p階龍格-庫塔方法。 其中ai，bij，ci為待定參數(shù)，要求上式y(tǒng)i+1在點(xi，yi)處作Taylor展開，通過相同項的系數(shù)確定參數(shù)。 一階龍格-庫塔公式，也就是歐拉法 二階龍格-庫塔公式，也就是改進的歐拉法,差分原理,差分方法是一類重要的數(shù)值解法。 尋求一系列離散節(jié)點 x1 x2 xn上的近似解y1，y2，yn，。 h=xn+1-xn稱為步長。 初值問題差分方法的特點： 步進式求解過程順著節(jié)點排列的次序一步一步地向前推進。 描述這種算法，只要給出從已知信息yn，yn-1，

6、yn-2 ，計算yn+1的遞推公式差分格式。 求解的核心消除導數(shù)，離散化方法。,差分原理,有限差分法(finite difference method) 微分方程和積分方程數(shù)值解的方法。 基本思想把彈性力學的基本方程和邊界條件(一般均為微分方程)近似地改用差分方程(代數(shù)方程)來表示，把求解微分方程的問題改換成為求解代數(shù)方程的問題，屬于數(shù)學上的近似。(有限單元法是以有限個單元的集合體來代替連續(xù)體，屬于物理上的近似)。 有限差分法的主要內(nèi)容包括： 根據(jù)問題的特點將定解區(qū)域作網(wǎng)格剖分 把原微分方程離散化為差分方程組 解差分方程組 此外為了保證計算過程的可行和計算結(jié)果的正確，還需從理論上分析差分方程組

7、的性態(tài)，包括解的唯一性、存在性和差分格式的相容性、收斂性和穩(wěn)定性,有限差分法簡介,構(gòu)造差分的方法有多種形式，目前主要采用的是泰勒級數(shù)展開方法。其基本的差分表達式主要有三種形式： 一階向前差分 一階向后差分 一階中心差分和二階中心差分 其中前兩種格式為一階計算精度，后兩種格式為二階計算精度。通過對時間和空間這幾種不同差分格式的組合，可以組合成不同的差分計算格式。 向前差分 向后差分 中心差分 上面談的是一階導數(shù)，對應的稱為一階差分。,有限差分法簡介,對一階差分再作一階差分，得到二階差分，記為D2y。 以向前差分為例，有 依此類推，任何階差分都可由其低一階的差分再作一階差分得到。 例如n階前差分為

8、,差分原理,以一維的情況為例 向前差分 一階差分 二階差分 向后差分 一階差分 二階差分 中間差分 一階差分 二階差分,差分原理,函數(shù)的差分與自變量的差分之比，即為函數(shù)對自變量的差商。 一階向前差商為 一階向后差商為 一階中心差商為 或,差分原理,二階差商多取中心式，即 對于多元函數(shù)f(x,y,)的差分與差商也可以類推。如一階向前差商為,差分原理,截斷誤差 差商與導數(shù)之間的誤差表明差商逼近導數(shù)的程度，稱為截斷誤差。由函數(shù)的Taylor展開，可以得到截斷誤差相對于自變量差分(增量)的量級，稱為用差商代替導數(shù)的精度，簡稱為差商的精度。 一階向后差商也具有一階精度。,差分原理,將f(x+Dx)和f(

9、x-Dx)的Taylor展開式相減可得 可見一階中心差商具有二階精度。 可見二階中心差商的精度也為二階精度。,差分原理,非均勻步長 在有些情況下要求自變量的增量本身是變化的，相應的差分和差商就是不等距的。 一階向后差商 一階中心差商,差分原理,彈性體上，用相隔等間距h而平行于坐標軸的兩組平行線織成正方形網(wǎng)格，Dx=Dy=h，如圖。 設(shè)f=f(x，y)為彈性體內(nèi)的某一個連續(xù)函數(shù)。該函數(shù)在平行于x軸的一根網(wǎng)線上，如在3-0-1上，它只隨x坐標的改變而變化。在鄰近節(jié)點0處，函數(shù)f可展為泰勒級數(shù)如下：,應力函數(shù)的差分解,只考慮離開節(jié)點0充分近的那些節(jié)點，即(x-x0)充分小。于是可不計(x-x0)的三

10、次及更高次冪的各項，則上式簡寫為： 在節(jié)點3，x=x0-h， 在節(jié)點1， x=x0+h，代入上式得： 聯(lián)立上兩式，解得差分公式,應力函數(shù)的差分解,同理，在網(wǎng)線4-0-2上可得到差分公式 差分公式中分母為2h的是以相隔2h的兩節(jié)點處的函數(shù)值來表示中間節(jié)點處的一階導數(shù)值，可稱為中點導數(shù)公式。 以相鄰三節(jié)點處的函數(shù)值來表示一個端點處的一階導數(shù)值，可稱為端點導數(shù)公式。 中點導數(shù)公式與端點導數(shù)公式相比，精度較高。因為前者反映了節(jié)點兩邊的函數(shù)變化，而后者卻只反映了節(jié)點一邊的函數(shù)變化。因此，盡可能應用前者，而只有在無法應用前者時才不得不應用后者。,應力函數(shù)的差分解,以上4式是基本差分公式，從而可導出其它的差

11、分公式如下：,應力函數(shù)的差分解,當不計體力時，將彈性力學平面問題歸結(jié)為在給定邊界條件下求解雙調(diào)和方程的問題。用差分法解平面問題，就應先將雙調(diào)和方程變換為差分方程，而后求解。 一旦求得彈性體全部節(jié)點的j值后，就可按應力分量差分公式(對節(jié)點0)算得彈性體各節(jié)點的應力。,應力函數(shù)的差分解,可見，用差分法解平面問題，共有兩大任務： 建立差分方程 將(D.1)(D.3)代入雙調(diào)和方程 對于彈性體邊界以內(nèi)的每一節(jié)點，都可以建立這樣一個差分方程。,應力函數(shù)的差分解,聯(lián)立求解這些線性代數(shù)方程，就能求得各內(nèi)節(jié)點處的值 一般建立和求解差分方程，在數(shù)學上不會遇到很大困難。但是，當對于邊界內(nèi)一行的(距邊界為h的)節(jié)點

12、，建立的差分方程還將涉及邊界上各節(jié)點處的j值，并包含邊界外一行的虛節(jié)點處的j值。 為了求得邊界上各節(jié)點處的j值，須要應用應力邊界條件，即：,應力函數(shù)的差分解,由右圖可見 于是，上式可改寫為 由此得,應力函數(shù)的差分解,關(guān)于邊界上任一點處的j/y，j/x值，可將上式從A點到B點對s積分得到：,應力函數(shù)的差分解,由高等數(shù)學可知 將此式也從A點到B點沿s進行積分，就得到邊界上任一點B處的j值。為此利用分部積分法，得： 整理得： 可見，已知 即可根據(jù)面力分量及導數(shù)求得,應力函數(shù)的差分解,由分析可知，把應力函數(shù)加上一個線性函數(shù)，并不影響應力。因此，可設(shè)想把應力函數(shù)加上a+bx+cy，然后調(diào)整a，b，c三個

13、數(shù)值，使得 這樣,應力函數(shù)的差分解,從圖看出，前面3式右邊的積分式分別表示A與B之間的x方向的面力之和、A與B之間的y方向的面力之和、A與B之間的面力對于B點的矩。 至此，解決了計算邊界上各節(jié)點 的值的問題。,應力函數(shù)的差分解,邊界外一行虛節(jié)點處j的值，則可用邊界上節(jié)點處的j/y或j/x值和邊界內(nèi)一行相應節(jié)點處j的值來表示。例如，對于圖中的虛節(jié)點14，因為有 所以有 當求出全部節(jié)點上的j值以后，我們就 可按應力分量的差分公式計算應力分量。,應力函數(shù)的差分解,用差分法解彈性平面問題時，可按下列步驟進行： 1.在邊界上任意選定一個節(jié)點作為基點A，取 然后由面力的矩及面力之和算出邊界上所有各節(jié)點處j

14、的值，以及所必需的一些j/y及j/x值，即垂直于邊界方向的導數(shù)值。 2.應用公式(D.6)，將邊界外一行虛節(jié)點處的j值用邊界內(nèi)的相應節(jié)點處的j值來表示。 3.對邊界內(nèi)的各節(jié)點建立差分方程(D.5)，聯(lián)立求解這些節(jié)點處的值。,應力函數(shù)的差分解,4.按照公式(D.6)，算出邊界外一行的各虛節(jié)點處的j值。 5.按照公式(D.4)計算應力的分量。,應力函數(shù)的差分解,說明 以上是針對單連體導出的結(jié)果。 如果一部分邊界是曲線的，或是不與坐標軸正交，則邊界附近將出現(xiàn)不規(guī)則的內(nèi)節(jié)點。對于這樣的節(jié)點，差分方程(D.4)必須加以修正。,應力函數(shù)的差分解,差分方程的相容性、收斂性、穩(wěn)定性 相容性 對于一個微分方程建

15、立的各種差分格式，為了有實用意義，一個基本要求是它們能夠任意逼近微分方程。 導數(shù)與其差分近似式之間存在截斷誤差。因此，差分方程的解并不是嚴格的，而是近似地滿足原來的偏微分方程。但是，當時間步長Dt和空間步長Dx都趨近于零時，差分方程的截斷誤差也趨近于零，差分方程的極限形式就是原偏微分方程。這時，認為差分方程與偏微分方程是相容的，這種相容性表示差分方程“收斂”于原偏微分方程。是收斂性的必要條件。,有限差分法,收斂性 一個差分格式是否有用，最終要看差分方程的精確解能否任意逼近微分方程的解。 差分方程的解，即當步長Dt0 、Dx0時收斂于原偏微分方程的解。,有限差分法,穩(wěn)定性 差分格式的計算過程是逐

16、步推進的，在計算第n+1步的近似值時要用到第n步的近似值 ，直到與初始值有關(guān)。 前面各步若有舍入誤差，必然影響到后面各步的值，如果誤差的影響越來越大，以致差分格式的精確解的面貌完全被掩蓋，這種格式是不穩(wěn)定的。 相反如果誤差的傳播是可以控制的，或者說若這種誤差積累保持有界，就認為格式是穩(wěn)定的，只有在這種情形，差分格式在實際計算中的近似解才可能任意逼近差分方程的精確解。,有限差分法,優(yōu)點 差分法是解微分方程邊值問題和彈性力學問題的有效方法； 差分法簡便易行； 對于某些結(jié)構(gòu)，為了更精確地分析局部的應力狀態(tài)，可以用差分法進行分析； 缺點 對于曲線邊界和斜邊界等產(chǎn)生的不等間距網(wǎng)格的處理，比較麻煩和易于出錯； 比較適用于求解二維問題或平面問題； 比較適用于等間距網(wǎng)格，對于應力變化較為劇烈時，需要采用二次網(wǎng)格進行計算。,有限差分法,有限差分法,以位勢問題為例，例如對于各向同性二維穩(wěn)態(tài)導熱方程，如果物體內(nèi)部熱源密度為0，方程可寫成如下形式： 如方程和邊界條件如下，按等步長0.25，求解各節(jié)點的溫度。方程及邊界條件如下：,算例,算