消除多次透射公式高頻振蕩失穩(wěn)的一種措施1

董 青 周正華 蘇 杰 劉厚毅 宋加密

?

消除多次透射公式高頻振蕩失穩(wěn)的一種措施1

董 青 周正華 蘇 杰 劉厚毅 宋加密

（南京工業(yè)大學(xué)，交通運輸工程學(xué)院，南京 210009）

多次透射公式（MTF）物理概念簡單，便于在計算機上實現(xiàn)時空解藕的高精度波動數(shù)值模擬。然而，MTF與其它局部人工邊界條件類似，存在數(shù)值模擬失穩(wěn)問題，如高頻振蕩便是可能出現(xiàn)的失穩(wěn)現(xiàn)象。本文在分析MTF高頻振蕩失穩(wěn)機理的基礎(chǔ)上，提出了在波動有限元數(shù)值模擬中消除MTF高頻振蕩失穩(wěn)的一種措施，即在整個有限元數(shù)值模擬區(qū)內(nèi)施加與應(yīng)變速率成正比的較小粘性阻尼；同時，討論了這一穩(wěn)定措施的有效性及其對數(shù)值計算精度的影響，并通過數(shù)值試驗檢驗了這一穩(wěn)定措施的可行性。結(jié)果表明，消除高頻振蕩失穩(wěn)的措施行之有效，且只對波動有限元數(shù)值模擬中無意義的高頻分量具有抑制作用，而對有意義的較低頻段內(nèi)的波動有限元數(shù)值模擬精度影響較小。

單側(cè)波 粘性阻尼 數(shù)值模擬 散射問題 波源問題

引言

在無限介質(zhì)或半無限介質(zhì)的有限元波動數(shù)值模擬中，需人為引入一種人工邊界，以獲得有限計算模型。這種人為引入的邊界會在邊界上產(chǎn)生外行波的假反射。因此，需在人工邊界上建立人工邊界條件，以消除人工邊界產(chǎn)生的反射。自20世紀(jì)60年代末以來，國內(nèi)外學(xué)者針對全局與局部人工邊界已進行了大量研究，并獲得了一系列研究成果（廖振鵬，1997）。在已建立的人工邊界條件中，基于對單側(cè)波一般運動學(xué)特征一維描述的、直接模擬的多次透射公式（MTF），物理概念簡單，便于在計算機上實現(xiàn)時空解藕的高精度波動有限元或有限差分數(shù)值模擬（Liao等，1984；廖振鵬等，1984；Liao，1996；廖振鵬，1996）。

將集中質(zhì)量動力時域有限元方法與MTF結(jié)合，可實現(xiàn)時空解藕的波動數(shù)值模擬，且精度可控。與其它局部人工邊界條件一樣，穩(wěn)定實現(xiàn)MTF是波動有限元或有限差分數(shù)值模擬的研究重點。MTF結(jié)合集中質(zhì)量動力時域有限元方法模擬失穩(wěn)現(xiàn)象，包括高頻振蕩失穩(wěn)和低頻飄移失穩(wěn)，周正華等（2001）給出了消除多次透射公式飄移失穩(wěn)的措施。本文將對MTF在波動有限元數(shù)值模擬中引起的高頻振蕩失穩(wěn)問題進行討論，并提出相應(yīng)的消除高頻振蕩失穩(wěn)的措施，以穩(wěn)定實現(xiàn)MTF；在此基礎(chǔ)上，通過三維波動數(shù)值模擬，檢驗了提出的消除高頻振蕩失穩(wěn)措施的有效性。

1 多次透射公式高頻失穩(wěn)機理及穩(wěn)定措施

有關(guān)MTF的基本理論及基本公式的導(dǎo)出已有詳細介紹（Liao，1996），文中不再贅述。Liao等（1992）和廖振鵬等（1992）基于有限元離散網(wǎng)格對波動數(shù)值模擬影響的分析，導(dǎo)出MTF反射系數(shù)的一般表達公式，并由此解釋了MTF產(chǎn)生高頻振蕩失穩(wěn)的主要特征及其機理。此外，謝志南等（2008）通過分析一維波動模型的離散模型，進一步討論了MTF高頻振蕩失穩(wěn)機理。研究結(jié)果表明：在對波動有限元數(shù)值模擬或有限差分數(shù)值模擬有意義的頻段內(nèi)，MTF一般不會產(chǎn)生高頻振蕩失穩(wěn)，而對數(shù)值模擬無意義的高頻段則會產(chǎn)生高頻振蕩失穩(wěn)；MTF高頻振蕩失穩(wěn)的實質(zhì)是數(shù)值模擬中無意義的高頻波動在人為引入的人工邊界上發(fā)生反射放大，并在數(shù)值模擬有限區(qū)域內(nèi)，由人工邊界引起多次反射導(dǎo)致在人工邊界上發(fā)生不斷的反射放大，數(shù)值模擬值越來越大，最終產(chǎn)生高頻振蕩失穩(wěn)現(xiàn)象。在動力有限元或有限差分數(shù)值模擬中，MTF隨著計算時間的增長有可能會引起數(shù)值模擬結(jié)果的高頻振蕩失穩(wěn)，且振蕩的頻率遠遠超過了波動數(shù)值模擬有意義的頻段，這種現(xiàn)象一般由動力有限元或有限差分數(shù)值模擬引起。引起MTF高頻振蕩失穩(wěn)的高頻波動分量，可能源于入射波中包含此高頻成分，或因時域逐步數(shù)值積分舍入誤差而隨機產(chǎn)生的高頻成份。若在數(shù)值模擬的波動中不含數(shù)值模擬無意義的高頻波動分量，或這樣的波動高頻分量甚小，加之與計算隨機舍入誤差產(chǎn)生的波動高頻分量在形成之初其值甚小，則在波動有限元或有限差分數(shù)值模擬中，所有這些極其微小的高頻波動分量起初并不顯現(xiàn)，只有在計算區(qū)內(nèi)經(jīng)過多次反射放大后，才能在有限元或有限差分數(shù)值模擬的結(jié)果中顯示出來（圖1）。因此，區(qū)別于動力有限元穩(wěn)定條件不滿足引起的失穩(wěn)，由MTF引起的高頻振蕩失穩(wěn)并非在數(shù)值模擬之初就開始，而是經(jīng)歷一段時間逐步放大直至失穩(wěn)。基于對MTF高頻振蕩失穩(wěn)的上述認識，認為可通過引入高頻濾波措施消除MTF產(chǎn)生的高頻振蕩失穩(wěn)。廖振鵬（1996）提出了1種消除MTF高頻振蕩失穩(wěn)的平滑措施，即通過平滑方法來消除波動的高頻分量。由于平滑濾波措施是在MTF波動高頻成份形成之后再作濾波，其濾波措施效果不總有效，并不能完全消除MTF在數(shù)值模擬中引起的高頻振蕩失穩(wěn)。鑒于平滑濾波措施不夠理想，關(guān)慧敏等（1997）提出了加邊界阻尼層的高頻濾波措施，但這一措施對于抑制MTF波動高頻成份仍不理想。

圖1 MTF數(shù)值模擬高頻振蕩失穩(wěn)

基于平滑濾波與加邊界阻尼層措施的啟發(fā)，嘗試通過引入措施對MTF波動高頻成份在形成之初就進行壓制，避免因人工邊界對數(shù)值模擬中無意義的高頻波動成份的放大而影響數(shù)值模擬結(jié)果，并通過在整個計算區(qū)內(nèi)引入阻尼效應(yīng)以消除高頻振蕩失穩(wěn)。大量的數(shù)值模擬結(jié)果顯示，在動力有限元數(shù)值模擬中，若考慮分析模型介質(zhì)的阻尼性質(zhì)，則能有效地抑制MTF引起的高頻振蕩失穩(wěn)（楊宇等，2014），同時亦可通過具有高頻能耗特性的積分格式消除MTF引起的高頻振蕩失穩(wěn)（李小軍等，2007；唐暉等，2010；楊宇等，2014）。因此，在分析模型中引入阻尼效應(yīng)或積分格式的能耗特性，可將無意義的波動高頻成份在形成之初消除。

已有的阻尼模型研究結(jié)果表明，粘性阻尼與應(yīng)變速率成正比，且隨著頻率變高而變大，尤其對高頻波動具有較好的抑制效能，但對波動的低頻成份影響卻很小。因此，本研究將在整個模型計算區(qū)內(nèi)施加粘性阻尼，濾掉數(shù)值模擬中不需要考慮的波動高頻成份，來實現(xiàn)MTF在動力有限元或有限差分數(shù)值模擬中的應(yīng)用。若在數(shù)值模擬分析模型中已具有此類阻尼，且介質(zhì)阻尼夠大，MTF高頻振蕩失穩(wěn)現(xiàn)象會自動得到抑制，無需再施加粘性阻尼；若分析模型中無此類阻尼，則可適當(dāng)施加，其實施辦法是在內(nèi)節(jié)點運動方程中，加上用單元阻尼陣[C]所表示的阻尼力，[C]與有限單元剛度矩陣[k]成正比，即：

式中，*為參考圓頻率（rad/s），取值通常為大于波動模擬有意義的截止頻率；為無量綱阻尼系數(shù)，取值為小正數(shù)，決定了引入阻尼的大小；[C]為單元阻尼就矩陣；[k]為單元剛度矩陣。

針對引入阻尼效應(yīng)消除MTF高頻振蕩失穩(wěn)的措施，可通過模態(tài)分析對集中質(zhì)量時域動力有限元方法內(nèi)節(jié)點計算精度的影響進行簡單說明。由模態(tài)分析可知，動力有限元內(nèi)節(jié)點的運動可通過疊加各振型運動而求得。對于給定的和*值，第振型的阻尼比d可表示為：

式中，為第振型的固有圓頻率，小于*。當(dāng)取值很小，且*取值比高時，阻尼比d很小，則對動力有限元數(shù)值模擬精度的影響很小，可忽略不計。因此，引入與應(yīng)變速率成正比的粘性阻尼這一措施，只對波動數(shù)值模擬中無意義的高頻分量具有較強的抑制作用，而不影響動力有限元數(shù)值模擬中有意義的較低頻段內(nèi)的計算精度。

2 數(shù)值試驗

本文采用時域波動有限元數(shù)值模擬的時空解耦方法（廖振鵬等，1984；Liao，1998），即MTF與集中質(zhì)量有限元相結(jié)合的顯式時域逐步積分方法，通過典型算例的穩(wěn)定實現(xiàn)，驗證所提措施對MTF引起的高頻振蕩失穩(wěn)的抑制效果，同時證明該措施對數(shù)值模擬精度影響很小，工程意義上可忽略不計。

分析研究基于三維波源問題的典型算例，考慮在均勻、各向同性線彈性半無限介質(zhì)自由表面，作用一豎向近似脈沖力產(chǎn)生的波動。采用直角坐標(biāo)系，其中平面位于彈性半空間自由表面，軸垂直向下，且坐標(biāo)系的原點與豎向近似脈沖力的作用點相同。設(shè)定豎向近似脈沖力的作用方向與軸一致，幅值為104kN，歸一化豎向近似脈沖力的時間函數(shù)為（），脈沖寬度為0.4s，如圖2所示。通過傅里葉分析得到豎向近似脈沖力的截止頻率c約為10Hz。

圖2 歸一化豎向近似脈沖力時程

為檢驗消除MTF高頻振蕩失穩(wěn)措施對數(shù)值模擬計算精度的影響，本文通過擴大有限計算區(qū)域，獲得數(shù)值模擬精確解。擴大有限計算區(qū)域旨在減小甚至消除人工邊界的引入對數(shù)值模擬計算精度的影響，幾何尺寸為240m×240m×120m，其離散單元大小、模型介質(zhì)力學(xué)參數(shù)、數(shù)值積分方法、計算時間步距以及MTF的實施均與上述分析相同，且未引入與應(yīng)變速率成正比的粘性阻尼，作用力仍然在模型自由表面的中點，數(shù)值模擬結(jié)果如圖5所示。為方便對比分析，圖中還給出了采用消除MTF高頻振蕩失穩(wěn)措施的計算結(jié)果。若以最大體波波速p=350m/s計算，則在0.6s內(nèi)，作用點處的數(shù)值模擬結(jié)果未受到人工邊界的影響。圖5結(jié)果表明，2種分析模型的計算結(jié)果基本相同；數(shù)值精確解的峰值為0.043286m，而采用消除MTF高頻振蕩失穩(wěn)措施的峰值為0.043276m，相對誤差約為2.3‰，說明本文提出的措施對計算精度的影響很小，可以忽略不計。

圖3 豎向近似脈沖力作用下作用點處的位移反應(yīng)時程

圖4 引入措施后豎向近似脈沖力作用下作用點處的位移反應(yīng)時程

圖5 2種模型數(shù)值模擬結(jié)果對比

綜上所述，本文提出的消除MTF高頻振蕩失穩(wěn)措施行之有效，且對數(shù)值模擬分析精度影響甚微。

3 結(jié)語

本文通過對多次透射公式（MTF）高頻振蕩失穩(wěn)機理的討論，提出了在有限元或有限差分波動數(shù)值模擬中針對MTF高頻振蕩失穩(wěn)的一種消除措施，即在整個計算區(qū)域內(nèi)施加與應(yīng)變速率成正比的粘性阻尼；并以此為基礎(chǔ)，從理論上論證了這一穩(wěn)定措施的有效性及其對數(shù)值模擬精度的影響；最后，通過數(shù)值模擬分析進行了檢驗。結(jié)果表明：

（1）引入與應(yīng)變速率成正比的粘性阻尼，能有效地消除MTF引起的高頻振蕩失穩(wěn)，且這一消除措施只對波動數(shù)值模擬無意義的高頻分量具有抑制作用。

（2）若的取值很小，且參考圓頻率*大于波動數(shù)值模擬截止圓頻率c時，這一消除高頻振蕩失穩(wěn)措施對有意義的較低頻段內(nèi)的波動數(shù)值模擬精度影響很小，可忽略不計。

關(guān)慧敏，廖振鵬，1997．一種改善多次透射邊界穩(wěn)定性的措施．地震工程與工程振動，17（4）：1—8．

李小軍，唐暉，2007．結(jié)構(gòu)體系動力方程求解的顯式積分格式的能耗特征．工程力學(xué)，24（2）：28—33．

廖振鵬，黃孔亮，楊柏坡等，1984．暫態(tài)波透射邊界．中國科學(xué)：A輯，14（6）：556—564．

廖振鵬，劉晶波，1992．波動有限元模擬的基本問題．中國科學(xué)：B輯，22（8）：874—882．

廖振鵬，1996．工程波動理論導(dǎo)引．北京：科學(xué)出版社．

廖振鵬，1997．近場波動的數(shù)值模擬．力學(xué)進展，27（2）：193—216．

唐暉，李小軍，李真，2010．顯式積分格式對局部透射邊界高頻失穩(wěn)的抑制和消除作用．世界地震工程，26（4）：50—54．

謝志南，廖振鵬，2008．人工邊界高頻振蕩失穩(wěn)機理的一點注記．地震學(xué)報，30（3）：302—306．

楊宇，李小軍，賀秋梅等，2014．散射問題中消除多次透射邊界高頻振蕩失穩(wěn)措施比較分析．地震工程學(xué)報，36（3）：476—481．

周正華，廖振鵬，2001．消除多次透射公式飄移失穩(wěn)的措施力學(xué)學(xué)報，33（4）：550—554．

Liao Z. P., Wong H. L., 1984. A transmitting boundary for the numerical simulation of elastic wave propagation. International Journal of Soil Dynamics and Earthquake Engineering, 3(4): 174—183.

Liao Z. P., Liu J. B., 1992. Numerical instabilities of a local transmitting boundary. Earthquake Engineering & Structural Dynamics, 21(1): 65—77.

Liao Z. P., 1996. Extrapolation non-reflecting boundary conditions. Wave Motion, 24(2): 117—138.

Liao Z. P., 1998. A decoupling numerical simulation of wave motion. Developments in Geotechnical Engineering, 83: 125—140.

The Measure Against High Frequency Oscillating Instability of Multi-transmitting Formula

Dong Qing, Zhou Zhenghua, Su Jie, Liu Houyi and Song Jiami

(College of Transportation Science and Engineering, Nanjing Tech University, Nanjing 210009, China)

Multi-Transmitting Formula (MTF) based on one-dimensional description on general kinematic characteristics of one-way wave and direct numerical simulation has simple physical concepts, and it is easy to realize the high precision and decoupling wave motion numerical simulation on the computer. However, Multi-Transmitting Formula, similar as other local artificial boundary conditions, is also a local artificial boundary condition which existed numerical instability problems, in which high-frequency oscillation instability is one of the instability phenomenon. Through the discussion of the mechanism of high frequency oscillation instability of Multi-Transmitting Formula, we present a measure to eliminate high frequency oscillation instability in wave numerical simulation. Through numerical simulation, the feasibility of the stabilization measure is verified. The results show that this stabilization measure is effective and of significantly cutting effectiveness on high-frequency component that is meaningless for wave motion numerical simulation. However, its effect on the wave motion numerical simulation was too small to ignore in the lower frequency band.

One-way wave; Viscous damping; Numerical simulation; The scattering problem; Source problem

董青，周正華，蘇杰，劉厚毅，宋加密，2018．消除多次透射公式高頻振蕩失穩(wěn)的一種措施．震災(zāi)防御技術(shù)，13（3）：571—577．

10.11899/zzfy20180308

國家自然科學(xué)基金項目（41374049）

2018-02-11

董青，女，生于1992年。碩士研究生。主要從事巖土力學(xué)研究。E-mail：2458810997@qq.com

周正華，男，生于1962年。研究員。研究領(lǐng)域：防災(zāi)減災(zāi)與防護工程。E-mail：bjsmoc@163.com