時間域邊界元素法於二維多層次暫態波傳分析研究(I)

行政院國家科學委員會專題研究計畫成果報告

時間域邊界元素法於二維多層次暫態波傳分析研究（Ⅰ）

Time domain BEM for two-dimensional multi-region problems of

transient wave propagation（Ⅰ）

計畫編號：NSC 88-2211-E-009-017

執行期限：88 年 8 月 1 日至 89 年 7 月 31 日

主持人：劉俊秀 交通大學土木研究所

一、中文摘要 本計畫主要在利用二維時間域邊界元 素法來求解多層次（Multi-region）的問題。 而在推導式子的過程中，我們藉由次結構 物（Substructure）的觀念求解。在求解的 過程中，在不同區域的交界面上，有二項 條件必須滿足，一為力平衡關係，另一為 位移諧和關係。而為了驗証其精確性，我 們去比較結果的數值解和其理論解，由誤 差可知，用時間域邊界元素法來求解二維 多層次的波傳問題，可得到很好的數值 解。 關鍵詞：邊界元素法、多層次、次結構物 Abstr act

The project is devoted to solve the problems of wave propagation in 2-D multi-region using time-domain BEM. In the process of formulating 2-D time-domain BEM for multi-region. The concept of substructure is employed. In the substructure

formulations, two conditions, one is

equilibrium condition and the other is compatibility condition, must be applied at the interface of different regions. In order to demonstrate the accuracy of the presented method, the results were compared with some simple theoretical solutions. These comparisons show the presented method can give good numerical solutions for wave propagating in 2-D multi-region. Keywor ds: multi-region、BEM、substructure 二、目的 建立一以時間域邊界元素法來分析求 解二維多層次暫態波傳問題即為本計畫一 系列研究之目的。若將一區域內部切割成 若干個區域，則在交界面上有下列兩項條 件成立，一為力平衡關係，另一為位移諧 和關係。我們就是應用此兩項條件，且配 合時間域邊界元素法來求解邊界值（曳引 力、位移）。並且將我們所得到的數值解 和其理論解來比較對照，可知用時間域邊 界元素法來求解二維多層次的波傳問題， 可得到很好的數值解。 最後再利用得到的結果所畫出的位移 波形，來解釋波傳由一介質至另一不同介 質時的物理現象。 三、結果與討論 　　 首先我們先定義一參數β： l t C ∆ = 1 β 上式中，l 為一元素之長度，C 為壓力波傳₁ 速度， t∆ 為時階（Time Step）。在文獻[14] 中建議我們取用的β 值最好在 0.5~0.75，其 所得到的解會較為精確，若β 值超過此範 圍，過大或過小，所得到的解將會不精確， 甚至發散。 接著我們再建立出分析的問題模型： 我 們 討 論 一 矩 形 桿 件 ， 長 與 高 的 比 值 L/h=10，左端為固定端，右端為自由端， 在我們分析的例子中，可將左端的固定端 視成中間為簡支承（Hinge），上下兩端為 滾支承（Roller），且假設矩形桿上下表面 為自由表面，如圖 1 所示，且我們將矩形

桿件切割成 13 個二次邊界元素，其中兩端 各用一個元素，上下各 5 個元素，交界面 用一個元素，如圖 2 所示。 最後，若二塊為同一材質，且β 值在 0.5~0.75 時，受一單位靜力時，其 x 和 y 方向位移其數值解和理論解的誤差皆＜ 1%，若受動力時，剛開始的誤差為 2~3%， 但歷時越久，數值累積的誤差會越大。若 二塊為不同材質時，我們所採用β 值在 0.5~0.75 的那塊其誤差約為 3~4%，而另一 塊因β 值不落在 0.5~0.75 的範圍內，故其 值會較不準。由其數值解和理論解的比 較，在受靜力時，精度相當的好，在受動 力時，精度的誤差也在可接受的範圍內。 接下來，我們將再利用我們得到的結 果來畫出位移波形，並藉由正確的波傳物 理現象來反証我們所得到的結果是合理 的。而為了要畫出 x 方向較多點的位移波 形，故我們選用如圖 3 的切割模型，上、 下 各 20 個 元 素 ， 長 與 高 的 比 值 亦 是 L/h=10，左端為固定端，右端為自由端， 第一塊的剪力模數 G 為 6.24Pa（C1 = 40 m/sec），第二塊 G 為 2.4375Pa（C1 = 25 m/sec ） ， 二 塊 質 量 密 度ρ皆 為 0.0078 kg/m2 ，包松比ν皆為 0，，二塊交界面是在 點號 17 的位置，且受一動力 P（t）作用， P（t）如圖 4 所示。 入射波的傳遞方向是由第二塊傳入第 一塊，當入射波剛好到達二塊交界面時， 即如圖 5 所示，在時階 N=20 時，剛好到達 交界面，會產生一反射波，其產生位置剛 好在交界面位置 0.4 之處，且此時反射波振 幅為正，故會使位置 0.4 以後的各點 x 方向 位移上升，然而在事實上，上升之後的每 點其 x 方向的位移應都相同，會為一直線， 但在圖 5 上，在 N=20 尖峰處後會略為下降 後再保持一直線，會有此一差距，可能是 數值上的誤差所造成的。 當入射波剛到達交界面時，會產生一 反射波往返回去，即如同上面所說的，會 使位置 0.4 以後的各點 x 方向位移上升，另 外也產生一穿透波繼續向前傳，但因為穿 透波是在較剛性的介質，故穿透波所經過 的點 x 方向的位移會下降，且保持一直線， 如圖 6 所示。且在圖 6 中，N=21 之後，即 入射波剛過交界面後，可明顯的看出來波 行的方向，即反射波會往回走，和穿透波 會繼續向前傳的現象。 因為穿透波在 G 較大的介質中，相較 於反射波在 G 較小的介質，其波速會較 快，比較圖 6 和圖 7，可明顯的看出，當 N=27 時，穿透波已快到達固定端，即穿透 波從 N=21→N=27，走了將近 13 個點，但 反射波從 N=21→N=27，只走了約 8 個點。 又因為我們此一模型是受一三角衝擊波， 故穿透波在遇固定端反射往回後，反射的 穿透波過了之後的點，因後面沒有波再經 過，故其 x 方向的位移即會慢慢減少趨近 於 0，由圖 7 中，可明顯的看出，反射後的 穿透波遇固定端後，往回走的現象，且波 走了之後的那些點，其 x 方向的位移會慢 慢減少趨近於 0。而此時原先的入射波在交 界面產生的反射波約走到位置 0.75 處，距 自由端尚有一段距離。 由上面敘述可知，當穿透波遇固定端 後，會反射往回傳，如圖 8 所示，當 N=36 時，反射的穿透波又剛好遇到交界面，因 為是由較剛性的介質到較軟的介質，反射 的穿透波的反射係數 Cr為負值，可知反射 的穿透波所產生的反射波其振幅為負，造 成了圖 8 中位置 0.4 附近的 x 方向的位移為 負值。亦可從圖中可明顯的看出其反射波 產生的位置亦在位置 0.4 的交介面上，此時 原先的反射波已快要到達自由端。 由圖 9 可明顯的看出，當 N=37 時，原 先的反射波已到達自由端，其反射係數 Cr 為-1，穿透係數 Ct為 0，故原先的反射波在 自由端為一全反射現象，但振幅由正變為 負，所以原先的反射波在自由端所產生的 反射波，會使自由端附近點的 x 方向的位 移下降，從圖中的 N=37，N=38，N=39， N=40 可明顯的視出此一象現。此時穿透波

在交界面所產生的反射波由圖可知，又快 回到了固定端。 四、計畫成果自評 每一工作項目皆按預定進度百分之百 完成，並獲得重要結果，值得繼續進一步 深入相關研究與探討。 五、參考文獻

1. Jawson,M.A. ( 1963 ) , ’Integral equation methods in potential theory. ’ I. Proc. Roy. Soc. London A275,23-32.

2. Symm,G.T. ( 1963 ) , ’Integral equation methods in potential theory. ’ II. Proc. Roy. Soc. London A275,33-46.

3. C.A.Brebbia（1988）, ’Boundary elements X,vol.4：Geomechanics , Wave Propagation and Vibrations,’ Springer-Verlag , Berlin,Heidelberg.

4. G.D.Manolis and D.E.Beskos （1988）, ’Boundary Element Methods in Elastodynamics,’ Unwin-Hyman,London. 5. Y.Niwa,T.Fukui,S.Kito and K.Fujiki,’An

application of the integral equation method to two-dimensional elastodynamics,’ Theor.Appl.Mech.Tokyo Univ.,28 ： 281-290,1980.

6. W.J.Mansur, ’A time-stepping technique to solve wave propagation problems using the boundary element method,’ Ph.D.Thesis,Southampton University,1983. 7. H.Antes,’A boundary element procedure for

transient wave propagation in two-dimensional istropic elastic media,’Finite Elem.Anal.Des.,1：313-322,1985.

8. A.S.M.Israil and P.K.Banerjee, ‘Advanced time-domain formulation of BEM for two-dimensional transient elastodynamics,’Int.J.Numer.Meth.Eng.,29: 1421-1440,1990.

9. A.S.M.Israil and P.K.Banerjee,’Two-dimensional transient wave propagation problems by time-domain BEM,’Int.J.Solids Struct.’,26：851-864,1990.

10. A.S.M.Israil and P.K.Banerjee, ‘Interior stress calculations in 2-D time-domain transient BEM analysis,’ Int.J.Solids Struct.’,27：915-927,1991.

11. P.K.Banerjee and R.Butterfield (1981) , ‘Boundary Element Method in Engineering Science,’ McGraw-Hill,London and New York.

12. J.D.ACHENBACH (1980) , ‘Wave Propagation in Elastic Solids,’ The Technological Institute, Northwestern University, Evanston, Illinois.

13. 王忠成,“推導高階時間域邊界元素法求 解二維暫態應力波”,國立交通大學土木 工程研究所博士論文,1996. 14. 賴進益,“二維時間域邊界元素法：內應 力之求解與元素切割方式之探討”, 國立 交 通 大 學 土 木 工 程 研 究 所 碩 士 論 文,1998. 六、附圖： L=10 P(t) h=1 A B C 圖 1 矩形桿之固定端簡化後之模式 P(t) 圖 2 一矩形桿切割成 13 個元素之示意圖 1 2 17 41 P(t) １ ２ 圖 3 一矩形桿切割成 43 個元素之 示意圖與點號標示

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Position (x/L) 0.00 0.40 0.80 1.20 D is p la c e m e n t (u E /P L ) β=0.67 G 6.24-2.4375 ν 0-0 N=17 N=18 N=19 N=20 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Position (x/L) -0.40 0.00 0.40 0.80 1.20 D is p la c e m e n t (u E /P L ) β=0.67 G 6.24-2.4375 ν 0-0 N=33 N=34 N=35 N=36 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Position (x/L) -0.40 0.00 0.40 0.80 1.20 D is p la c e m e n t (u E /P L ) β=0.67 G 6.24-2.4375 ν 0-0 N=36 N=37 N=38 N=39 N=40 1.6 Time(Sec) 0.0 10 P(t) . 2∆t ∆t = 0.0133 圖 4 一三角形之衝擊載重，P(t)之形式 圖 5 受一三角形之衝擊載重，入射波 至交界面前之位移波形 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Position (x/L) 0.00 0.40 0.80 1.20 D is p la c e m e n t (u E /P L ) β=0.67 G 6.24-2.4375 ν 0-0 N=26 N=27 N=28 N=29 圖 8 受一三角形之衝擊載重，穿透波的 反射波回至交界面時之位移波形 圖 6 受一三角形之衝擊載重，入射波 至交界面時之位移波形 圖 9 受一三角形之衝擊載重，原先入射波 的反射波至自由端時之位移波形 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Position (x/L) 0.00 0.40 0.80 1.20 D is p la c e m e n t (u E /P L ) β=0.67 G 6.24-2.4375 ν 0-0 N=20 N=21 N=22 N=23 圖 7 受一三角形之衝擊載重，穿透波 至固定端後反射之位移波形