一種彈性動(dòng)力學(xué)邊界面計(jì)算方法
【專利摘要】本發(fā)明涉及一種彈性動(dòng)力學(xué)邊界面計(jì)算方法，針對(duì)傳統(tǒng)時(shí)域法計(jì)算量大、分析時(shí)間較長(zhǎng)時(shí)結(jié)果不穩(wěn)定的問(wèn)題，提出了一種擬初始條件法和卷積數(shù)值積分法相結(jié)合的方法來(lái)求解三維彈性動(dòng)力學(xué)問(wèn)題的時(shí)域邊界積分方程。在新的方法中，首先利用擬初始條件法將整個(gè)分析時(shí)間劃分為幾個(gè)時(shí)間段，將每個(gè)時(shí)間段末的動(dòng)態(tài)響應(yīng)作為下一段的初始條件，不必再累加之前分析步的結(jié)果對(duì)當(dāng)前分析步的影響，從而可以減少系數(shù)矩陣的存儲(chǔ)量和計(jì)算時(shí)間，降低了計(jì)算量；然后在每個(gè)時(shí)間段內(nèi)，再次劃分為更小的時(shí)間步采用卷積數(shù)值積分法逐步求解。該方法既能夠提高傳統(tǒng)時(shí)域法的穩(wěn)定性，又克服了卷積數(shù)值積分法計(jì)算量大的問(wèn)題。
【專利說(shuō)明】
一種彈性動(dòng)力學(xué)邊界面計(jì)算方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明涉及計(jì)算機(jī)和力學(xué)數(shù)值分析技術(shù)，具體涉及一種彈性動(dòng)力學(xué)邊界面計(jì)算方 法。
【背景技術(shù)】
[0002] 隨著計(jì)算機(jī)技術(shù)的迅速發(fā)展，數(shù)值計(jì)算方法使得我們能夠?qū)?fù)雜的工程問(wèn)題進(jìn)行 仿真計(jì)算，從而得到其數(shù)值解。除有限元法以外，邊界元法也是常用的數(shù)值計(jì)算方法之一， 特別是在處理無(wú)限域波動(dòng)、裂紋擴(kuò)展以及應(yīng)力集中等問(wèn)題時(shí)邊界元法具有天然的優(yōu)勢(shì)。然 而，無(wú)論是邊界元法還是有限元法，都需要單獨(dú)生成一個(gè)離散化的網(wǎng)格近似模型才能進(jìn)行 分析，存在著幾何模型和分析模型不一致的問(wèn)題。
[0003] 基于邊界積分方程和計(jì)算機(jī)圖形學(xué)，參考文獻(xiàn)《基于邊界面法的完整實(shí)體應(yīng)力分 析理論與應(yīng)用》(作者:張見(jiàn)明，期刊:計(jì)算機(jī)輔助工程，2010,19(3) :5-10)提出了一種邊界 面法。邊界面法具備邊界元法的所有優(yōu)點(diǎn)，同時(shí)又克服了幾何模型和分析模型不一致的問(wèn) 題。它直接利用幾何造型系統(tǒng)中的邊界表征數(shù)據(jù)結(jié)構(gòu)對(duì)模型進(jìn)行離散和分析，能夠精確計(jì) 算分析過(guò)程中所有的幾何數(shù)據(jù)，有效避免了幾何誤差的產(chǎn)生，是一種適應(yīng)設(shè)計(jì)/分析一體化 發(fā)展趨勢(shì)的新型數(shù)值分析方法。
[0004] 根據(jù)對(duì)時(shí)間處理方式的不同，彈性動(dòng)力學(xué)問(wèn)題的邊界積分方程可以用以下三種方 法求解:時(shí)域法、頻域法(Laplace變換法或Four ier變換法)和雙重互易法。以下是彈性動(dòng)力 學(xué)問(wèn)題的時(shí)域邊界積分方程：
[0005] 對(duì)于均勻、各向同性的線彈性材料而言，彈性動(dòng)力學(xué)問(wèn)題的位移運(yùn)動(dòng)方程(亦稱為 控制方程)如下：
[0006]
(1)
[0007]式中i，j = l，2,3;Ui代表位移分量，Uj,ij(q，t)為t時(shí)亥Ijj方向位移的二階導(dǎo)數(shù)，uMj (q，t)為t時(shí)刻i方向位移的二階導(dǎo)數(shù)代表單位質(zhì)量的體積力;P是材料密度A是位移對(duì) 時(shí)間的二階導(dǎo)數(shù)，也就是加諫庠:λ和G是拉梅常數(shù)，表達(dá)式如下：
[0_
(2)
[0009] 分別表示材料的壓縮性和剪切模量，E為彈性模量，V為泊松比。
[0010] 若雎件體可弄示為加圖1所示的區(qū)域V，方程（1)滿足邊界條件：
[0011] (3)
[0012]
[0013] (4)
[0014] 其中，q表示彈性體上任意一點(diǎn)，可V)和MV)分別表示邊界P和於上已知的位移 和面力，?
分別為初始位移和初始速度。
[0015] 控制方程(1)對(duì)應(yīng)的時(shí)域邊界積分方程為：
[0016]
[0017] 其中》:(P，qU-T)是位移基本解，表示τ瞬時(shí)作用在ρ點(diǎn)i方向的單位集中力脈沖所引 起的q點(diǎn)t瞬時(shí)j方向的位移分量為是面力基本解； 是速度基本解，為位移基本對(duì)時(shí)間的一次導(dǎo)數(shù)。當(dāng)源點(diǎn)P位于光滑邊界時(shí)(^(ρ) = δ^/2,位 于域內(nèi)時(shí)Cij(p) = 3ij。
[0018] 相對(duì)而言，時(shí)域法是最為自然直接且應(yīng)用范圍最廣的一種方法，它可以直接得到 時(shí)間相關(guān)的動(dòng)態(tài)響應(yīng)。特別是處理爆炸、沖擊等瞬態(tài)問(wèn)題、接觸等非線性問(wèn)題以及幾何體形 狀隨時(shí)間變化的動(dòng)邊界問(wèn)題時(shí)，時(shí)域法比其他方法更加有利。然而，當(dāng)分析時(shí)間較長(zhǎng)時(shí)，時(shí) 域法可能會(huì)出現(xiàn)結(jié)果不穩(wěn)定的現(xiàn)象。Schanz曾將卷積數(shù)值積分法應(yīng)用到時(shí)域邊界積分方程 中，有效提高了計(jì)算結(jié)果的穩(wěn)定性。該方法采用有限項(xiàng)求和近似計(jì)算邊界積分方程中基本 解和節(jié)點(diǎn)值之間的卷積。它采用頻域基本解，同時(shí)又可以直接得到真正的時(shí)域解，不需要將 邊界條件進(jìn)行變換，也不需要將結(jié)果進(jìn)行反變換，這是一種介于時(shí)域法和頻域法之間的一 種數(shù)值計(jì)算方法。它繼承了傳統(tǒng)時(shí)域法的優(yōu)點(diǎn)，同時(shí)又能夠提高時(shí)域法的穩(wěn)定性。但是也存 在著一個(gè)問(wèn)題，由于卷積數(shù)值積分法在計(jì)算系數(shù)矩陣時(shí)是對(duì)復(fù)數(shù)進(jìn)行操作，又要進(jìn)行傅里 葉變換，所以其計(jì)算量和存儲(chǔ)量都比傳統(tǒng)時(shí)域法大很多。

【發(fā)明內(nèi)容】

[0019] 本發(fā)明提供了一種彈性動(dòng)力學(xué)邊界面計(jì)算方法，以解決現(xiàn)有技術(shù)中的彈性動(dòng)力學(xué) 邊界面計(jì)算方法計(jì)算穩(wěn)定性差且計(jì)算量、存儲(chǔ)量較大的缺陷。
[0020] 為解決上述問(wèn)題，本發(fā)明的彈性動(dòng)力學(xué)邊界面計(jì)算方法包括如下步驟：
[0021 ] 1)將所求彈性動(dòng)力學(xué)問(wèn)題的整個(gè)分析時(shí)間[0，t ]劃分為M個(gè)時(shí)間段，將對(duì)每個(gè)時(shí)間 段進(jìn)行彈性動(dòng)力學(xué)邊界積分方程的求解過(guò)程稱為一個(gè)子分析，采用擬初始條件法處理；
[0022] 2)將每一個(gè)子分析時(shí)間段[t^U]進(jìn)一步劃分為N個(gè)分析步;根據(jù)時(shí)域邊界積分 方程，利用卷積數(shù)值積分法求出第一個(gè)子分析的每一個(gè)分析步的動(dòng)態(tài)響應(yīng)及N組權(quán)系數(shù)矩 陣；
[0023] 3)將第一個(gè)子分析的最后一個(gè)分析步計(jì)算出來(lái)的動(dòng)態(tài)響應(yīng)中的節(jié)點(diǎn)位移和速度 作為第二個(gè)子分析的虛擬初始位移和速度，調(diào)用第一個(gè)子分析的N組權(quán)系數(shù)矩陣，根據(jù)時(shí)域 邊界積分方程，利用卷積數(shù)值積分法求出第二個(gè)子分析的每一個(gè)分析步的動(dòng)態(tài)響應(yīng)，以此 類推，最終求出M個(gè)子分析的每一個(gè)分析步的動(dòng)態(tài)響應(yīng)，實(shí)現(xiàn)彈性動(dòng)力學(xué)邊界面計(jì)算。
[0024] 在利用卷積數(shù)值積分法進(jìn)行計(jì)算前，將當(dāng)前子分析的節(jié)點(diǎn)位移和速度分別轉(zhuǎn)換為 單位質(zhì)量初始位移的等效虛擬力和單位質(zhì)量初始速度的虛擬力，利用單位質(zhì)量初始位移的 等效虛擬力和單位質(zhì)量初始速度的虛擬力，將時(shí)域邊界積分方程轉(zhuǎn)換為初始位移加上基于 位移增量的邊界積分方程。
[0025] 所述時(shí)域邊界積分方程為：
[0026]
[0027] 其中，<(P，q;z-『)是位移基本解，表示τ瞬時(shí)作用在P點(diǎn)i方向的單位集中力脈沖所引起 的q點(diǎn)t瞬時(shí)j方向的位移分量為< O是面力基本解；〖〇似-〇) = M(P，叫_〇)?， 是速度基本解，為位移基本解對(duì)時(shí)間的一次導(dǎo)數(shù)。
[0028] 所怵擬系教的i+笪公忒加下，
[0029]
[0030]
[0031] 其中，4
為線性多步法中的二次向后差分公式， 為誤差控制參數(shù)，計(jì)算公式為 分別為拉普拉斯變換后的位移基本解和面力基本解。
[0032] 所述基于位移增量的邊界積分方程為：
[0033] ?=()
[0034] 其中
為權(quán)系數(shù)；APj(q，tm-i+aAt) 為邊界S上任意一點(diǎn)的面力增量，Auj(q，tm-i+aA t)為邊界S上任意一點(diǎn)的位移增量，bj(q, tm-1+a At)為單位質(zhì)量的體積力，〇)為單位質(zhì)量初始位移的等效虛擬力，廣(q)為單位 質(zhì)量初始速度的等效虛擬力，P是材料密度。
[0035]所述單位質(zhì)量初始位移的等效虛擬力的公式為：
[0036]
[0037] 分別表示材料的壓縮性和剪 切模量，E為彈性模量，V為泊松比，H( t-Ο)為Heavi s ide函數(shù)，P是材料密度，Uj, ij (q，0)為初 始時(shí)亥?方向位移的二階導(dǎo)數(shù)，m,ij(q，0)為初始時(shí)刻i方向位移的二階導(dǎo)數(shù)。
[0038] 所述單位質(zhì)量初始速度的虛擬力的公式為：= vIl(tI)沖)，其中，δ(?)是Dirac delta函數(shù)，r〉'(q)為q點(diǎn)在i方向的初始速度。
[0039] 本發(fā)明的有益效果:針對(duì)傳統(tǒng)時(shí)域法計(jì)算量大、分析時(shí)間較長(zhǎng)時(shí)結(jié)果不穩(wěn)定的問(wèn) 題，本發(fā)明提出了一種擬初始條件法和卷積數(shù)值積分法相結(jié)合的方法來(lái)求解三維彈性動(dòng)力 學(xué)問(wèn)題的時(shí)域邊界積分方程。在新的方法中，首先利用擬初始條件法將整個(gè)分析時(shí)間劃分 為幾個(gè)時(shí)間段，將每個(gè)時(shí)間段末的動(dòng)態(tài)響應(yīng)作為下一段的初始條件，不必再累加之前分析 步的結(jié)果對(duì)當(dāng)前分析步的影響，從而可以減少系數(shù)矩陣的存儲(chǔ)量和計(jì)算時(shí)間，降低了計(jì)算 量;然后在每個(gè)時(shí)間段內(nèi)，再次劃分為更小的時(shí)間步采用卷積數(shù)值積分法逐步求解，該方法 既能夠提高傳統(tǒng)時(shí)域法的穩(wěn)定性，又克服了卷積數(shù)值積分法計(jì)算量大的問(wèn)題。
[0040] 無(wú)論是擬初始條件法引入的虛擬初始條件，還是工程問(wèn)題中實(shí)際存在的初始位移 和初始速度，都在邊界積分方程中涉及到體積分的計(jì)算。由于彈性動(dòng)力學(xué)問(wèn)題的時(shí)域基本 解含有Dirac Delta函數(shù)，而且邊界積分方程中的初始條件項(xiàng)不含有卷積，無(wú)法使用卷積數(shù) 值積分法，所以導(dǎo)致初始條件項(xiàng)的數(shù)值體積分變得非常復(fù)雜。針對(duì)該問(wèn)題，本發(fā)明提出了一 種適用于三維線彈性動(dòng)力學(xué)問(wèn)題的虛擬力法，基于線性系統(tǒng)的疊加原理以及沖量定理將初 始條件轉(zhuǎn)化為等效的體積力來(lái)處理，原控制方程被轉(zhuǎn)換為一個(gè)初始條件完全為零的控制方 程，從而避免了復(fù)雜的幾何求交操作以及球面積分的產(chǎn)生，大大簡(jiǎn)化了初始條件項(xiàng)的計(jì)算。
【附圖說(shuō)明】
[0041 ]圖1為彈性體邊界條件示意圖；
[0042]圖2為彈性動(dòng)力學(xué)邊界面計(jì)算方法流程圖；
[0043] 圖3為懸臂桿幾何模型和邊界條件示意圖；
[0044] 圖4(a)為外力作用下懸臂桿自由端長(zhǎng)時(shí)間的位移響應(yīng)結(jié)果對(duì)比整體曲線圖；
[0045] 圖4(b)為外力作用下懸臂桿自由端長(zhǎng)時(shí)間的位移響應(yīng)結(jié)果對(duì)比局部放大曲線圖； [0046]圖5(a)為外力作用下懸臂桿固定端長(zhǎng)時(shí)間的面力響應(yīng)結(jié)果對(duì)比整體曲線圖；
[0047]圖5(b)為外力作用下懸臂桿固定端長(zhǎng)時(shí)間的面力響應(yīng)結(jié)果對(duì)比局部放大曲線圖；
[0048]圖6為計(jì)算時(shí)間隨分析步數(shù)的變化對(duì)比圖。
【具體實(shí)施方式】
[0049]下面結(jié)合附圖，對(duì)本發(fā)明的技術(shù)方案作進(jìn)一步詳細(xì)說(shuō)明。
[0050]如圖2所示，本實(shí)施例的彈性動(dòng)力學(xué)邊界面計(jì)算方法的步驟如下：
[0051 ] 1)將所求彈性動(dòng)力學(xué)問(wèn)題的整個(gè)分析時(shí)間[0，t ]劃分為M個(gè)時(shí)間段，將對(duì)每個(gè)時(shí)間 段進(jìn)行彈性動(dòng)力學(xué)邊界積分方程的求解過(guò)程稱為一個(gè)子分析；
[0052] 2)將每一個(gè)子分析時(shí)間段[t^U]進(jìn)一步劃分為N個(gè)分析步;根據(jù)時(shí)域邊界積分 方程，利用卷積數(shù)值積分法求出第一個(gè)子分析的每一個(gè)分析步的動(dòng)態(tài)響應(yīng)及N組權(quán)系數(shù)矩 陣；
[0053] 3)將第一個(gè)子分析的最后一個(gè)分析步計(jì)算出來(lái)的動(dòng)態(tài)響應(yīng)中的節(jié)點(diǎn)位移和速度 作為第二個(gè)子分析的虛擬初始位移和速度，調(diào)用第一個(gè)子分析的N組權(quán)系數(shù)矩陣，根據(jù)時(shí)域 邊界積分方程，利用卷積數(shù)值積分法求出第二個(gè)子分析的每一個(gè)分析步的動(dòng)態(tài)響應(yīng)，以此 類推，最終求出M個(gè)子分析的每一個(gè)分析步的動(dòng)態(tài)響應(yīng)，實(shí)現(xiàn)彈性動(dòng)力學(xué)邊界面計(jì)算。
[0054] 下面詳細(xì)闡述上述步驟：
[0055] 首先，利用擬初始條件法，將整個(gè)分析時(shí)間[0，t]劃分為M個(gè)時(shí)間段，稱每個(gè)時(shí)間段 的分析為子分析。在每個(gè)子分析內(nèi)，將上一個(gè)子分析的計(jì)算結(jié)果作為當(dāng)前子分析的初始條 件，那么當(dāng)前子分析相應(yīng)的時(shí)域邊界積分方程為：
[0056]
[0057] 和速度。
[0058]然后將每一個(gè)子分析的時(shí)間段[t^U]進(jìn)一步劃分為N個(gè)更小的分析步，時(shí)間步 長(zhǎng)設(shè)為△ t，用卷積數(shù)值積分法求出每一個(gè)子分析步的動(dòng)態(tài)響應(yīng)，時(shí)間1?散后的方程為：
[0059:
[0060：
[0061：
[0062:
[0063]其中，γ (Z) = 1.5-2z+0.5z2為線性多步法中的二次向后差分公式；當(dāng)L = N+1以及 紀(jì)=石時(shí)，權(quán)系數(shù)ω n-a的誤差量級(jí)為〇(')（1〇3 SidO111)，《(p，q.i)和Κ(Μ，?)分別是拉普拉 斯變換后的位移基本解和面力基本解。其中，權(quán)系數(shù)ω "1的值只與分析步之間的時(shí)間差n-a 的大小有關(guān)，與第幾個(gè)子分析m無(wú)關(guān)，因此只需在第一個(gè)子分析中將N組系數(shù)矩陣計(jì)算出來(lái) 并存儲(chǔ)，之后的所有子分析計(jì)算時(shí)只要調(diào)用即可，不需要重新計(jì)算。系數(shù)矩陣的存儲(chǔ)量由原 來(lái)的M X N組降為N組，相應(yīng)的計(jì)算時(shí)間也能大幅減少。
[0064] 方程(7)等號(hào)右邊的后兩項(xiàng)代表的是在上一個(gè)子分析的最后一步計(jì)算出來(lái)的節(jié)點(diǎn) 位移和速度，我們將其作為當(dāng)前子分析的虛擬初始位移和初始速度。這兩項(xiàng)體積分不包括 時(shí)間卷積，因此不能使用卷積數(shù)值積分法進(jìn)行計(jì)算。為此，我們提出了一種適用于三維彈性 動(dòng)力學(xué)問(wèn)題的虛擬力法，將初始條件轉(zhuǎn)換為等效的體積力處理，具體過(guò)程為：
[0065] 首先進(jìn)行初始位移的等效。對(duì)于線彈性系統(tǒng)，可以設(shè)控制方程（1)的解等于初始位 移加上位移增量，如
[0066] Ui(q,t)=Ui(q,0)+Aui(q,t) (10)
[0067] 將公式（10)代入方程(1)可得：
[0068]
[0069] 定義單位質(zhì)量初始位移的等效虛擬力為：
[0070]
(12)
[0071]之所以定義為單位質(zhì)量的虛擬力，是為了與控制方程中的體積力單位保持一致。 這里H(t-〇)為Heaviside函數(shù)，只有當(dāng)t>0的時(shí)候才有值。當(dāng)t>0時(shí)，可以將上式簡(jiǎn)寫(xiě)為：
[0072]
Π3)
[0073] 若要求出域內(nèi)任意一點(diǎn)初始位移虛擬力fWq)的值，可以通過(guò)以下兩種方法：
[0074] (1)當(dāng)m(q，0)的函數(shù)表達(dá)式已知時(shí)，可以直接對(duì)m(q，0)進(jìn)行求導(dǎo)，然后求出相應(yīng) 的fiu°(q)的值。
[0075] (2)當(dāng)m(q，0)無(wú)法用具體的函數(shù)表達(dá)時(shí)，也就是說(shuō)初始位移的分布是任意的、沒(méi) 有規(guī)律的，那么可以通過(guò)求解下列Navier方程來(lái)得到fWq)的值。
[0076] 上述方程為彈性靜力學(xué)控制方 程，與傳統(tǒng)的控制方程求解不一樣，在這個(gè)方程中，初始位移m(q，0)為已知量初始位移，而 體積力fWd)為需要求解的未知變量。
[0077] 接下來(lái)進(jìn)行初始速度的等效。由動(dòng)量定理"物體動(dòng)量的增量等于它所受合外力的 沖量"可知，從t = (T到t = 0時(shí)刻物體動(dòng)量的變化可以等效為單位質(zhì)量的沖擊力fV^qA)作 用在物體上，即
[0078]
(15}
[0079] 其中，衫(!,『) =·()。也就是說(shuō)物體一開(kāi)始處于靜止?fàn)顟B(tài)，直到?jīng)_擊力作用后才產(chǎn)生了 初始速度。由于作用時(shí)間很短，物體的位移還沒(méi)來(lái)得及變化。將初始速度~他〇) = ^如)代入 方程(15)可以得到單位質(zhì)量初始速度的虛擬力為
[0080]
(16)
[0081] 其中δ(?)是Dirac delta函數(shù)，代表./TdO除了t = 0時(shí)其他時(shí)刻值都為零。如果我 們將沖擊力的作用平均到第一個(gè)時(shí)間步[0, At]內(nèi)，就可以得到不含時(shí)間變量的初始速度 虛擬力^
[0082] (17)
[0083]最后，將初始位移和初始速度的等效虛擬力（12)和（16)代入原始控制方程(1)，可 以得到一個(gè)新的關(guān)于位移增量的控制方程：
[0089]原始控制方程（1)的解就等于新的控制方程（18)的解△ Ui(q，t)加上初始位移Ui
[0084]
[0085]
[0086]
[0087]
[0088] (q，〇) ο
[0090] 將上述虛擬力法應(yīng)用于方程(7)，則方程（7)的解等于初始位移uKpAh)加上下 面關(guān)于位移增量的邊界積分方程的解：
[0091]
[0092]
[0093]
[0094]
[0095] 本發(fā)明使用邊界面法進(jìn)行空間變量的離散和積分，邊界面法是直接基于邊界表征 實(shí)體造型數(shù)據(jù)結(jié)構(gòu)實(shí)現(xiàn)其網(wǎng)格離散及單元積分的，實(shí)體的邊界曲面都以參數(shù)形式表達(dá)。在 其實(shí)現(xiàn)中，無(wú)論是邊界的數(shù)值積分還是場(chǎng)變量的插值都是在邊界曲面的二維參數(shù)空間里進(jìn) 行。該方法首先在曲面的二維參數(shù)空間里將每個(gè)參數(shù)曲面離散成若干個(gè)參數(shù)曲面單元，這 些單元同樣都是定義在二維參數(shù)空間，而非三維物理空間，也就相當(dāng)于CAD造型系統(tǒng)中的分 片曲面，從而曲面原始的幾何信息都能夠保留下來(lái)。積分過(guò)程使用到的幾何數(shù)據(jù)，比如高斯 積分點(diǎn)的三維坐標(biāo)、雅可比、外法線向量都是直接由參數(shù)曲面計(jì)算獲得，而不是通過(guò)分段多 項(xiàng)式插值近似，從而避免了幾何誤差。
[0096] 將物體表面離散為BN個(gè)表面節(jié)點(diǎn)，DN個(gè)域內(nèi)節(jié)點(diǎn)。對(duì)于任意一個(gè)源點(diǎn)qk，（ k = 1， 2,--,BN或者k = l .2. ···.),空間富散后方趕（21 )可弄示為加下鉭陳方趕，
[0097]
[0098] 其中m=1，2,…，Μ，η = 1，2，···，Ν.系數(shù)矩陣只需要在第一個(gè)子分析中求出 即可，只要離散的分析步長(zhǎng)不變，就不需要重新計(jì)算系數(shù)矩陣。求出當(dāng)前子分析中每一個(gè)分 析步每個(gè)節(jié)點(diǎn)的位移增量'H1 +?△〇后，加上初始位移Uj (q，tm-i)的值就是每個(gè)節(jié)點(diǎn)真 正的位移響應(yīng)。當(dāng)一個(gè)子分析求解結(jié)束后，還需要求解出最后一個(gè)分析步所有域內(nèi)節(jié)點(diǎn)的 位移和速度，為下一個(gè)子分析做準(zhǔn)備。初始位移的計(jì)算直接求解域內(nèi)邊界積分方程即可，在 域內(nèi)邊界積分方程中Cij (p) = Sij(當(dāng)i = j時(shí)，Sij = 1，否則Sij = 〇)，源點(diǎn)p取遍所有域內(nèi)節(jié)點(diǎn)。 節(jié)點(diǎn)速度的計(jì)算需要用到最后一步以及前一步的節(jié)點(diǎn)位移，設(shè)在每一個(gè)分析步內(nèi)位移呈線
f生變化，抓Λ晶臼一來(lái)的書(shū)占迪IiPff曾公才加下
[0099] (25)
[0100] 重復(fù)以上過(guò)程計(jì)算下一個(gè)子分析，直到M個(gè)子分析全部結(jié)束。加入擬初始條件法之 后，系數(shù)矩陣的存儲(chǔ)量降為原來(lái)的1/M，雖然引入了體積分，但是計(jì)算時(shí)間仍然比之前要少， 劃分的子分析數(shù)越多，內(nèi)存需求和計(jì)算時(shí)間就越短。但是子分析數(shù)也不能過(guò)多，否則會(huì)影響 計(jì)算精度。
[0101]下面通過(guò)一個(gè)算例驗(yàn)證上述方法的正確性，對(duì)擬初始條件法的計(jì)算效率和計(jì)算精 度進(jìn)行研究。所使用微機(jī)的處理器為3.4GHZ的Intel(R)C〇re(TM)i7-2600 CPU，內(nèi)存空間為 8GB〇
[0102] 計(jì)算懸臂桿的受迫振動(dòng)，懸臂桿幾何模型及尺寸如圖3所示，長(zhǎng)L = 8m，寬和高均為 2m。左端面固定，右端面沿X軸承受一 Heavi side類型的均布拉力p = 1000N/m2。材料參數(shù)為： 彈性模量E = I. I e5N/m2，密度P = 2. Okg/m3，泊松比υ = 〇. 〇。本算例暫不考慮重力的影響，之 所以取零泊松比，是為了將數(shù)值結(jié)果和解析解進(jìn)行比較，在程序中仍然將該問(wèn)題當(dāng)作三維 問(wèn)題來(lái)計(jì)算。此工況下，懸臂桿上任意一點(diǎn)任意時(shí)刻的位移響應(yīng)解析解可以用下面的公式 計(jì)算：
[0103]
[0104]
[0105]
[0106]
[0107] 在這個(gè)算例中，采用線性四邊形單元和線性六面體單元離散模型，網(wǎng)格的單元尺 寸均為lm。加入擬初始條件法之前，分析步數(shù)如果超過(guò)400步，就會(huì)出現(xiàn)計(jì)算機(jī)內(nèi)存不足的 情況。加入擬初始條件法之后，這里將分析時(shí)長(zhǎng)延伸到1.424s，時(shí)間步長(zhǎng)仍然為0.00089s， 總的分析步數(shù)為1600步，下面來(lái)驗(yàn)證本發(fā)明方法的性能：
[0108] 首先將整個(gè)分析過(guò)程劃分為M個(gè)子分析，在每個(gè)子分析中，分析步個(gè)數(shù)為1600/M。 為了考察不同M的取值對(duì)結(jié)果的影響，分別取M=IO和M=20的數(shù)值結(jié)果與解析解進(jìn)行對(duì)比。 當(dāng)然，我們也考察了很多其它取值，計(jì)算結(jié)果相差并不特別明顯，為了圖片能夠顯示的清 楚，只取其中的兩組值來(lái)觀察。
[0109] 同時(shí)，我們也實(shí)現(xiàn)了Schanz的Reformulated CQM(改寫(xiě)的卷積數(shù)值積分法），在 Reformulated CQM中，只需開(kāi)辟一組復(fù)數(shù)系數(shù)矩陣的存儲(chǔ)空間，但是每一步需要重新計(jì)算 一次單元積分，未知量的系數(shù)矩陣也需要重新求逆。這個(gè)方法和頻域法更加相似，相對(duì)來(lái)說(shuō) Reformulated CQM的計(jì)算量比傳統(tǒng)的卷積數(shù)值積分法要小，但是卻喪失了時(shí)域法的主要優(yōu) 勢(shì)。為了對(duì)比加入擬初始條件法前后卷積數(shù)值積分法的計(jì)算結(jié)果，我們將Reformulated CQM的結(jié)果也加入到對(duì)比圖片中。
[0110]圖4顯示了 1.424s內(nèi)懸臂桿自由端位移隨時(shí)間的變化，圖中曲線1代表一維桿縱向 振動(dòng)的解析解，曲線2代表Reformulated CQM的數(shù)值解，曲線3、4代表卷積數(shù)值積分法和擬 初始條件法相結(jié)合的數(shù)值解（convolution quadrature method and pseudo-initial condition method，縮寫(xiě)為CQM&PICM)?？梢钥吹?，完全理想的工況下解析解一直保持同樣 的振幅和頻率，三種數(shù)值解都是穩(wěn)定的，和解析解吻合較好，但是均有不同程度的衰減。為 了更清楚的觀察數(shù)值解之間的差別，我們將前1/6和后1/6時(shí)間段的結(jié)果放大，圖4(b)為圖4 (a)的部分截?cái)喾糯髨D，可以看到:在分析的一開(kāi)始，三種數(shù)值解的精度都較高，隨著分析時(shí) 間的增加，結(jié)果誤差越來(lái)越大;相對(duì)而言，Reformulated CQM的數(shù)值解精度高于后兩種數(shù)值 解，也就是說(shuō)擬初始條件法的加入會(huì)對(duì)卷積數(shù)值積分法的結(jié)果精度有所影響，因?yàn)樵谔摂M 力法的應(yīng)用中會(huì)引入額外的誤差;子分析個(gè)數(shù)M=IO的數(shù)值解比M = 20的數(shù)值解精度高，可 見(jiàn)劃分的子分析個(gè)數(shù)越少，對(duì)卷積數(shù)值積分法的結(jié)果精度影響就越小。同樣的現(xiàn)象在固定 端面力的動(dòng)態(tài)響應(yīng)結(jié)果中也可以觀察得到，如圖5所示。
[0111]接下來(lái)對(duì)比三種數(shù)值解法的計(jì)算時(shí)間。本算例中統(tǒng)計(jì)的均為直接計(jì)算的時(shí)間，并 未加入快速算法。圖6顯示了三種數(shù)值解法的計(jì)算時(shí)間隨分析步數(shù)的變化對(duì)比， Reformulated CQM隨分析步數(shù)的增加呈線性增加趨勢(shì)，當(dāng)分析步數(shù)繼續(xù)增加時(shí)，該方法的 計(jì)算時(shí)間將會(huì)非?？捎^;擬初始條件法中M=IO的計(jì)算時(shí)間為M = 20的兩倍，因?yàn)榍罢咝枰?計(jì)算并存儲(chǔ)的系數(shù)矩陣為后者的兩倍。當(dāng)總的分析步數(shù)較少時(shí)，本算例中擬初始條件法的 計(jì)算時(shí)間比Reformulated CQM的計(jì)算時(shí)間長(zhǎng)，這是因?yàn)樵谶@個(gè)例子中，初始條件和體積力 均為零，Reformulated CQM的邊界積分方程中不含有體積分項(xiàng)，而擬初始條件法的加入需 要增加體積分的計(jì)算，增加了系數(shù)矩陣的計(jì)算時(shí)間。但是隨著分析步數(shù)的增加，擬初始條件 法的計(jì)算時(shí)間增加非常緩慢，可見(jiàn)該方法能夠顯著減少卷積數(shù)值積分法的計(jì)算時(shí)間，提高 彈性動(dòng)力學(xué)分析的計(jì)算效率。
[0112]表1對(duì)比列出了三種數(shù)值解法中系數(shù)矩陣所需的內(nèi)存空間。因?yàn)镽eformulated CQM只需要存儲(chǔ)一組系數(shù)矩陣，所以其所需的內(nèi)存空間最小，但是該方法不能算是真正的時(shí) 域法，因?yàn)樾枰獙?duì)邊界條件和計(jì)算結(jié)果進(jìn)行傅里葉變換，只能用于邊界條件隨時(shí)間變化已 知的問(wèn)題。
[0113]表1系數(shù)矩陣所需的存儲(chǔ)空間對(duì)比
[0115]從以上算例分析可知，當(dāng)分析步數(shù)較多時(shí)，擬初始條件法能夠顯著降低傳統(tǒng)時(shí)域 法的計(jì)算時(shí)間和系數(shù)矩陣的存儲(chǔ)量，與卷積數(shù)值積分法相結(jié)合能夠分析較長(zhǎng)時(shí)間的振動(dòng)響 應(yīng)。在擬初始條件法中劃分的子分析數(shù)越多，計(jì)算效率越高，但是結(jié)果精度會(huì)有所降低。所 以，我們應(yīng)該根據(jù)實(shí)際情況，如分析時(shí)長(zhǎng)、節(jié)點(diǎn)個(gè)數(shù)、計(jì)算機(jī)能力等合理使用擬初始條件法， 劃分適當(dāng)?shù)淖臃治鰯?shù)來(lái)平衡計(jì)算時(shí)間和計(jì)算精度之間的關(guān)系。
【主權(quán)項(xiàng)】
1. 一種彈性動(dòng)力學(xué)邊界面計(jì)算方法，其特征在于，該方法包括如下步驟： 1) 將所求彈性動(dòng)力學(xué)問(wèn)題的整個(gè)分析時(shí)間[〇，t]劃分為M個(gè)時(shí)間段，將對(duì)每個(gè)時(shí)間段進(jìn) 行彈性動(dòng)力學(xué)邊界積分方程的求解過(guò)程稱為一個(gè)子分析，采用擬初始條件法處理； 2) 將每一個(gè)子分析時(shí)間段[t^U]進(jìn)一步劃分為N個(gè)分析步;根據(jù)時(shí)域邊界積分方程， 利用卷積數(shù)值積分法求出第一個(gè)子分析的每一個(gè)分析步的動(dòng)態(tài)響應(yīng)及N組權(quán)系數(shù)矩陣； 3) 將第一個(gè)子分析的最后一個(gè)分析步計(jì)算出來(lái)的動(dòng)態(tài)響應(yīng)中的節(jié)點(diǎn)位移和速度作為 第二個(gè)子分析的虛擬初始位移和速度，調(diào)用第一個(gè)子分析的N組權(quán)系數(shù)矩陣，根據(jù)時(shí)域邊界 積分方程，利用卷積數(shù)值積分法求出第二個(gè)子分析的每一個(gè)分析步的動(dòng)態(tài)響應(yīng)，以此類推， 最終求出M個(gè)子分析的每一個(gè)分析步的動(dòng)態(tài)響應(yīng)，實(shí)現(xiàn)彈性動(dòng)力學(xué)邊界面計(jì)算。2. 根據(jù)權(quán)利要求1所述彈性動(dòng)力學(xué)邊界面計(jì)算方法，其特征在于，在利用卷積數(shù)值積分 法進(jìn)行計(jì)算前，將當(dāng)前子分析的節(jié)點(diǎn)位移和速度分別轉(zhuǎn)換為單位質(zhì)量初始位移的等效虛擬 力和單位質(zhì)量初始速度的虛擬力，利用單位質(zhì)量初始位移的等效虛擬力和單位質(zhì)量初始速 度的虛擬力，將時(shí)域邊界積分方程轉(zhuǎn)換為初始位移加上基于位移增量的邊界積分方程。3. 根據(jù)權(quán)利要求1所述彈性動(dòng)力學(xué)邊界面計(jì)算方法，其特征在于，所述時(shí)域邊界積分方 程為：其中，當(dāng)源點(diǎn)P位于光滑邊界時(shí)Ci j (p) = 5i」/2，位于域內(nèi)時(shí)Ci j (p) = 5i」，當(dāng)i = j時(shí)，& j = 1，否則 Sij = 〇;<(P，叫-0是位移基本解，表示謝時(shí)作用在p點(diǎn)i方向的單位集中力脈沖所引起的q點(diǎn)t 瞬時(shí)j方向的位移分量為< 叫-r)是面力基本解〇)/改，是 速度基本解，為位移基本解對(duì)時(shí)間的一次導(dǎo)數(shù)。4. 根據(jù)權(quán)利要求1所述彈性動(dòng)力學(xué)邊界面計(jì)算方法，其特征在于，所述權(quán)系數(shù)的計(jì)算公 式如下：Y (幻=1.5-22+0.522為線性多步法中的二次向后差分公式，貨為誤差控制參數(shù)，計(jì)算公式為： 為拉普拉斯變換后的位移基本解和面力基本解。5. 根據(jù)權(quán)利要求2所述彈性動(dòng)力學(xué)邊界面計(jì)算方法，其特征在于，所述基于位移增量的 邊界積分方程為：其中，％ (化]],A/)、％ "(A;A/)及叫(％.p,A/)為權(quán)系數(shù)；A pj(q，tm-i+a S上任意一點(diǎn)的面力增量，A uj(q，tm-l+a A t)為邊界S上任意一點(diǎn)的位移增量，bj(q，tm-l+a A t)為單位質(zhì)量的體積力，flq)為單位質(zhì)量初始位移的等效虛擬力，//Iq)為單位質(zhì)量初 始速度的等效虛擬力，P是材料密度。6. 根據(jù)權(quán)利要求2所述彈性動(dòng)力學(xué)邊界面計(jì)算方法，其特征在于，所述單位質(zhì)量初始位 移的等效虛擬力的公式為，分別表示材料的壓縮性和剪切模 量，E為彈性模量，v為泊松比，H(t-O)為Heaviside函數(shù)，P是材料密度，Uj,ij(q，0)為初始時(shí) 亥Ijj方向位移的二階導(dǎo)數(shù)，m,ij(q，0)為初始時(shí)亥Iji方向位移的二階導(dǎo)數(shù)。7. 根據(jù)權(quán)利要求2所述彈性動(dòng)力學(xué)邊界面計(jì)算方法，其特征在于，所述單位質(zhì)量初始速 度的虛擬力的公式為利0，其中，S(t)是Dirac delta函數(shù)，為q點(diǎn)在i方 向的初始速度。
【文檔編號(hào)】G06F17/50GK106055743SQ201610340671
【公開(kāi)日】2016年10月26日
【申請(qǐng)日】2016年5月20日
【發(fā)明人】李源, 孫林, 張見(jiàn)明, 王世勛, 張仕光, 王艷軍
【申請(qǐng)人】河南師范大學(xué)