TWI566106B

TWI566106B - 抑制管湧現象之模擬方法

Info

Publication number: TWI566106B
Application number: TW104141136A
Authority: TW
Inventors: 黃良雄; 林孟郁; 張芯瑜; 賴柏元
Original assignee: 黃良雄
Priority date: 2015-12-08
Filing date: 2015-12-08
Publication date: 2017-01-11
Also published as: TW201721461A

Description

抑制管湧現象之模擬方法

本發明是有關於一種抑制管湧現象之模擬方法，特別是指一種用於模擬一設置於流體中及透水底床上且具有一包覆一透水層的不透水管材之模型的抑制管湧現象之模擬方法。

近年來，水底管線廣泛被運用在資源的運輸，在相關單位定期的檢測維修下仍有不可預期的事故發生，造成管線失效的原因可分做兩種，自然因素：波浪與海流、未支撐之跨距、沉積物移動及地震。人為因素：船隻錨泊、漁撈活動及海上廢棄物等。底床產生沖刷現象主要原因為圓管下方發生管湧(piping)現象，也就是土粒被滲流水帶走，導致土中形成集中水流之通道如管狀。若能抑制管湧現象發生，可望能大量減少事後修復工作及成本。

已知近期文獻包括Talebbeydokhti & Afzali(2008年)使用二維有限元素法(FEM)及邊界積分法計算管線所受之抬升力及管面壓力分佈。結構物貫穿地下含水層將會嚴重影響地下水流流況，張正緯等人(2009年“三維地下水模式之發展與應用”，中國土木水利工程學刊，第二十一卷，第二期，頁169-181)乃提出建議，設計結構物可採用透水性較佳之建材。另外，水流流經多孔底床上不透水圓管，在圓管上下游底床產生極大壓力梯度，水流急遽增大，張芯瑜(2012年“底床上圓管管湧現象之減緩方法”，國立臺灣大學碩士論文)推測為管湧現象發生之原因，並證實於底床上圓管外包覆透水材質能有效降低圓管下緣之壓力梯度，進而減少管湧現象發生。因而，將圓管外包覆一透水層進行模擬證實可以使水流較不受圓管阻礙，透水底床降低之速度大小可達數十倍至數百倍，可有效減緩管湧現象之發生。

本發明之其中一目的，即在提供一種解決先前技術問題的抑制管湧現象之模擬方法。

本發明抑制管湧現象之模擬方法在一些實施態樣中是應用於模擬一模型的管湧現象，該模型包括在一流體中的一透水底床、一不透水圓管，及一包覆於該不透水圓管且接觸該透水底床之透水層，該方法由一具有一運算單元及一記憶單元之電子裝置執行以下步驟：(A)該運算單元載入對應該模型的邊界值問題之方程式及邊界條件；(B)該運算單元讀取該記憶單元預先儲存之多個物理參數並依據正規擾動法將該等物理參數換算為二個無因次化參數，該等無因次化參數包括該透水層的一第一滲透係數及該透水底床的一第二滲透係數；及(C)該運算單元判斷該第一滲透係數遠大於該第二滲透係數或該第一滲透係數遠小於該第二滲透係數時，依據該透水底床及該透水層交界面的邊界條件修正成具有通量之連續邊界條件，並配合正規擾動法的相關方程式以產生該模型對應之流速、流速勢及/或流體壓力之分佈。

在一些實施態樣中，該等方程式具有一第一方程式組、一第二方程式組及一第三方程式組，若為該第一滲透係數接近該第二滲透係數之狀態則執行該第一方程式組、若為該第一滲透係數遠大於該第二滲透係數之狀態則執行該第二方程式組，以及若為該第一滲透係數遠小於該第二滲透係數之狀態則執行該第三方程式組，以分別產生該模型處於各狀態對應之流速、流速勢及/或流體壓力之分佈。

於是，本發明在一些實施態樣中，當電腦載入該電腦程式並執行後，可完成如所述抑制管湧現象之模擬方法。

本發明至少具有以下功效：本發明之抑制管湧現象之模擬方法提出相關方程式的處理流程及實驗結果，除了可依據各種物理參數進行模擬以進行抑制管湧現象之模擬，並提出評估方程式以最佳化配置降低管湧現象。

1‧‧‧不透水圓管

100‧‧‧模型

2‧‧‧透水層

3‧‧‧流體

4‧‧‧透水底床

5‧‧‧電子裝置

50‧‧‧運算單元

51‧‧‧輸入單元

52‧‧‧記憶單元

53‧‧‧輸出單元

S101~S106‧‧‧步驟

本發明之其他的特徵及功效，將於參照圖式的實施方式中清楚地呈現，其中：圖1是本發明抑制管湧現象之模擬方法的實施例的流程圖；圖2是該實施例應用於模擬管湧現象的模型的示意圖；圖3是該模型之各交界的示意圖；圖4是該模型的透水底床的邊界上假想薄層的示意圖；圖5是該模型的透水底床的薄層厚度與角度之關係的曲線圖；圖6是該模型於的透水層邊界上假想薄層的示意圖；圖7是該模型的透水層邊界上假想薄層流況的示意圖；圖8是透水層假想薄層厚度與角度之關係的曲線圖；圖9及圖10是邊界積分計算範圍的示意圖；圖11是線性邊界元素的示意圖；圖12A是不透水圓管之模型計算的流程圖；圖12B是接續圖12A之底床上不透水圓管包覆透水層之模型計算的流程圖；圖13是無包覆透水層的不透水圓管代表切線速度計算位置的示意圖；圖14是該模型包覆透水層的不透水圓管代表切線速度計算位置的示意圖；圖15是測試案例的|u_imp/U|與t/T之關係的折線圖；圖16是測試案例的|u_p/U|與流體區入射流流速之關係的折線圖；圖17是測試案例的|u_imp/u_p|與流體區入射流流速之關係的折線圖；圖18是測試案例的|u_p/U|與R₁/R₂之關係的折線圖；圖19是測試案例的|u_imp/u_p|與R₁/R₂之關係的折線圖；圖20是測試案例的|u_p/U|與α之關係的折線圖；圖21是測試案例的|u_imp/u_p|與α之關係的折線圖；圖22是測試案例的|u_p/U|與ε ₁/ε ₂之關係的折線圖；圖23是測試案例的|u_imp/u_p|與ε ₁/ε ₂之關係的折線圖；圖24是測試案例修正後的|u_p/U|與ε ₁/ε ₂之關係的折線圖；圖25是測試案例修正後的|u_imp/u_p|與ε ₁/ε ₂之關係的折線圖；圖26是測試案例的Π與ε ₂之關係的折線圖；圖27是測試案例的χ與ε ₂之關係的折線圖；圖28是測試案例的Π、χ與α之關係的長條及折線圖；圖29是測試案例的Π、χ與α之關係的長條及折線圖；圖30是測試案例的Π、χ與ε ₁/ε ₂之關係的折線圖；圖31是測試案例的Π、χ與ε ₁/ε ₂之關係的折線圖；圖32是測試案例的Π、χ與R₁/R₂之關係的長條及折線圖；圖33是測試案例的Π、χ與R₁/R₂之關係的長條及折線圖；圖34是測試案例的χ與α之關係的折線圖；圖35是執行本發明抑制管湧現象之模擬方法的電子裝置的電路方塊圖。

一、本發明抑制管湧現象之模擬方法的處理流程

參閱圖1及圖2，本發明抑制管湧現象之模擬方法在一些實施態樣中應用於模擬一模型100的管湧現象，該模型100包括在一流體3中的一透水底床4、一不透水圓管1，及一包覆於該不透水圓管1且接觸該透水底床4之透水層2。其中，第I區是位在流體3中，第II區是位在透水層2中，第III區是位在透水底床4中。茲將本發明的處理流程及實驗結果說明如下。

S101一建立模型100的邊界值問題之方程式及邊界條件。本實施例中，該流體3為具有非黏性及不可壓縮性，流體3底部深處的水流平穩而不易產生渦流，故假設水流運動為非旋性。透水的多孔介質為剛性結構且內部流場滿足達西定律。另外，流動為非旋性，故第I~III區的速度可用一流速位勢(velocity potential)表示且界定邊界值問題的方程式如下。

第I區方程式包含：連續方程式▽． V ₁=0，及動量方程式；代入V _j=▽Φ_j，j=1,2,3成為連續方程式▽²Φ₁=0，及動量方程式

第II區方程式包含：連續方程式▽． V ₂=0，及動量方程式；代入V _j=▽Φ_j,j=1,2,3成為連續方程式▽²Φ₂=0，及動量方程式

第III區方程式包含：連續方程式▽． V ₃=0，及動量方程式；代入V _j=▽Φ_j,j=1,2,3成為連續方程式▽²Φ₃=0，及動量方程式

參閱圖3，為模型100具有的各種交界面的示意圖。

交界面Γ₁&Γ₃的方程式包含：設定輻射邊界條件(r→∞),→outgoing or 0,j=1,3；通量連續式 V ₁． n =n ₃ V ₃． n 改寫成，壓力連續式P ₁=P ₃改寫成

交界面Γ₂的方程式包含：通量連續式 V ₁． n =n ₂ V ₂． n 可改寫成，及壓力連續式P ₁=P ₂可改寫成

交界面Γ₄的方程式包含：通量連續： V ₂． n =0可改寫成

交界面Γ₅的方程式如公式7包含：通量連續式n₂ V ₂．n=n₃ V ₃．n可改寫成，壓力連續式P ₂=P ₃可改寫成

S102-建立多個物理參數並依據正規擾動法將該等物理參數換算為二個無因次化參數，該等無因次化參數包括該透水層2的一第一滲透係數ε ₁及該透水底床4的一第二滲透係數ε ₂，相關公式推導如下。

依據多孔介質流體3壓力表示式：

因為純流體3區與多孔介質區壓力連續邊界條件之等號兩邊量階相差過大，在數值計算求解係數矩陣中可能會遭遇不良條件，無法正確求解邊界值問題，需要借助正規擾動法解決此量階相異之問題。正規擾動法可將問題利用微小參數做擾動展開，將原本不適合求解之問題依量階分作若干個邊界值問題，最後將所有量階之解疊加即當作原問題之結果。

本實施例定義無因次化參數如下，其中P、Φ、μ分別為壓力、流速勢及動力黏滯係數之無因次參數。

公式9中，ρ是該流體之密度，g為重力加速度，R₁是該透水層表面至該不透水圓管中心之距離，將公式9帶入公式8可得無因次化後之公式10。

藉由公式10可選取兩無因次微小擾動參數如公式11。

第一滲透係數，第二滲透係數

由O(n ₂)&O(n ₃)=10^-1~10⁰；O(R ₁)=10^-1~10⁰(m) 公式12

可估計第一滲透係數及第二滲透係數之量階範圍約為O(ε₁),O(ε₂)<O(10^-8)≪1 公式13

S103-判斷該第一滲透係數遠大於該第二滲透係數，或該第一滲透係數遠小於該第二滲透係數？若否，則進入步驟S104；若是，則進入步驟S105。

S104-配合邊界元素法及正規擾動法的相關方程式以產生該模型100對應之流速、流速勢及/或流體壓力之分佈。

本實施例中係提出正規擾動之邊界值問題，區分為CASE A-兩介質(人工透水材質及天然多孔底床)的滲透係數相近、CASE B-人工透水材質滲透係數遠大於天然多孔底床滲透係數，以及CASE C-人工透水材質滲透係數遠小於天然多孔底床滲透係數的三種狀態。

CASE A：

本實施例於兩介質(人工透水材質及天然多孔底床)的滲透係數相近之狀態係採用一第一方程式組進行運算，第一方程式組具有的公式14~23說明如下。

令ε ₁=ε ₂=ε，利用正規擾動法將第I、Ⅱ、Ⅲ區的速度勢展開如下，Φ₁為該流體之速度勢，Φ₂為該透水層之速度勢，Φ₃為該透水底床之速度勢。

公式14

正規擾動法展開在O(1)階，，

正規擾動法展開在O(ε)階，，

&；

第Ⅱ區積分矩陣式為：

第Ⅲ區積分矩陣式：

聯立第Ⅱ、Ⅲ區積分矩陣式：

交換矩陣位置：

補充方程式成為：

S105-依據該透水底床4及該透水層2交界面的邊界條件修正成具有通量之連續邊界條件。

CASE B：

第一滲透係數遠大於第二滲透係數的狀態下，因為考量流體3與多孔介質區之壓力連續條件，故將壓力首階限制在第一階，若以速度勢展開將會得到首階相同之速度勢。本實施例於人工透水材質滲透係數遠大於天然多孔底床滲透係數之狀態係採用一第二方程式組進行運算，第二方程式組具有的公式24~32說明如下。

利用雙擾動微小參數將第I區速度勢展開

以上述原則將第Ⅱ、Ⅲ區壓力展開如下。

將第Ⅱ、Ⅲ區速度勢以壓力表示

通量連續邊界條件之修正，圓管外透水層2與底床交界面通量連續條件如公式27，N2、N3為正向單位向量(unit normal vector)。

將通量連續式轉換成壓力表示為：

以量階分析將得到無通量條件：&與現實不相符。

令代入修正，

建立底床上圓管管湧現象之邊界值問題時，原本圓管外透水層與底床交界面之通量連續邊界條件，現因ε ₁>>ε ₂之緣故，若依量階分析其應為無通量條件，但現實中流體還是會穿透交界面，兩多孔介質之間仍有通量交換，故將通量連續邊界條件作修正，使其不為無通量。

參閱圖4及圖5，由於底床多孔介質滲透係數遠小於人工透水層滲透係數，流體雖然可以通過交界面，但在邊界附近徑向速度分量相對於切線方向速度分量極小，本實施例假設在天然底床邊界上有一薄層，如圖4所示，此薄層厚度為人工透水層流體所能到達之最遠處，即該薄層上之徑向流速為零。

以，R₁=1.0...,α=90°及為例之計算結果，如圖5所示。結果顯示假想薄層厚度於θ=0度，即圓管中心正下方幾乎等於零，應為u_r(R₁,θ=0)之緣故。薄層裡之流場因距離邊界過近，若使用邊界元素法計算將有很大之奇異性。

正規擾動展開法在O(1)階，，

正規擾動展開法在O(ε ₁)階，&；

，

CASE C：

第一滲透係數遠小於第二滲透係數的狀態下，因為考量流體3與多孔介質區之壓力連續條件，故將壓力首階限制在第一階，若以速度勢展開將會得到首階相同之速度勢。本實施例於人工透水材質滲透係數遠小於天然多孔底床滲透係數之狀態係採用一第三方程式組進行運算，第三方程式組具有的公式33~41說明如下。

利用雙擾動微小參數將第I區速度勢展開

以上述原則將第Ⅱ、Ⅲ區壓力展開

將第Ⅱ、Ⅲ區速度勢以壓力表示

圓管外透水層2與底床交界面通量連續條件：

將通量連續式轉換成壓力表示為

以量階分析將得到無通量條件：&與現實不相符。

令代入做修正：

參閱圖6及圖7，由於人工透水層滲透係數遠小於底床多孔介質滲透係數，流體雖然可以通過交界面，但在邊界附近徑向速度分量相對於切線方向速度分量極小，本實施例假設在人工多孔介質邊界上有一薄層，如圖6所示，此薄層厚度為底床多孔介質流體所能到達之最遠處，即該薄層上之徑向流速為零。

參閱圖8，以，R₁=1.0...,α=90°及為例計算算薄層厚度，顯示假想薄層厚度於θ=0度，即圓管中心正下方幾乎等於零，應為u_r(R₁,θ=0)之緣故。薄層裡之流場因距離邊界過近，若使用邊界元素法計算將有很大之奇異性。

正規擾動展開法在O(1)階，，

正規擾動展開法在O(ε ₂)階，&，

，

二、數值模擬方法

邊界元素法為模擬方法中較後期發展出來的方法，其優點是只需對計算領域的表面做元素分割和數值計算，即所求解問題之維度可以降低一階，大大減少數值計算的複雜性。邊界元素有常數元素、線性元素、二次元素或更高次元素，本實施例採用線性元素，將邊界上的物理量視為線性變化，並可以較為合理地處理不連續邊界條件且提供相當足夠的精度。然而，利用邊界元素法求解分階分區之邊界值問題數值解，評估各物理參數對管湧作用減緩之影響，因為純流體區與多孔介質區壓力連續邊界條件之等號兩邊量階相差過大，在數值計算求解係數矩陣中可能會遭遇不良條件，面臨無法正確求解邊界值問題。

邊界元素法為模擬方法中較後期發展出來的方法，其優點是只需對計算領域的表面做元素分割和數值計算，即所求解問題之維度可以降低一階，大大減少數值計算的複雜性。參閱圖9及圖10，Ω為計算域，Γ為計算域之表面，為場點(field point)座標，可以在計算域的任一位置，為源點(source point)座標，可以在計算域表面的任一位置，為場點與原點之距離，為表面法向量，計算內域時法向量指向外，如圖9。計算外域時法向量指向內，如圖10。以下利用格林定理轉換 Laplace方程式為邊界積分式，將Laplace方程式式(3.1)之未知物理量函數代換成格林函數(Green’s function)如公式42。

Laplace方程式：

公式42中為任一可二次微分函數，δ為Dirac delta函數，求解公式42可得Laplace方程式之基本解(fundamental solution)如公式43及公式44。

格林第二恆等式：

前述公式改寫後可得邊界積分方程式如公式45，為一隨場點之內角變化函數，其值為場點之內角徑度2π與之比值，例：若場點位於平滑邊界上，。

參閱圖11，以下將以線性元素切割邊界，線性元素將邊界上的物理量視為線性變化，並使用一維無因次化座標取代原本二維卡式座標。邊界元素有常數元素、線性元素、二次元素或更高次元素，本實施例採用線性元素，將邊界上的物理量視為線性變化，並可以較為合理地處理不連續邊界條件且提供相當足夠的精度。J表示原始元素長度與轉換後元素長度比值。

各元素座標範圍定為-1 η 1

定義線性內插函數：

元素上座標及物理量可用線性內插函數表示為

每一邊界元素座標轉換須在積分式內乘上一Jacobian行列式值

本實施例並採用二維Laplace線性邊界元素法，將邊界切割為N個線性元素，則域內點積分方程式改寫為公式49

將其中積分式表示為

最後將域內點積分方程式寫為

三、測試案例說明

參閱圖12A及12B，分別為不透水圓管之模型計算及底床上不透水圓管包覆透水層之模型計算的流程圖。

參閱圖13，顯示無透水層2之切線速度計算位置，圖14則為顯示透水層2之切線速度計算位置。以下說明各項物理參數對管湧作用減緩之影響及最佳配置。

參閱表1，為本實施例之模型100採用的物理參數。

以下將根據表1所列之測試案例模擬結果討論各物理參數影響管湧現象或透水材質減緩效用。

(一)時間/週期

考慮週期性入射流與圓管之交互關係，故將速度勢函數乘上一餘弦函數cos(ω t)，讓速度勢函數可反映週期性運動流場特性，流場隨時間而改變。若每次模擬皆須計算一個完整週期的流況，將是很不經濟的作法，故將以不透水圓管1的切線速度u_imp最大之時刻作為模擬時間代表，圖15為測試案例t/T=0~1之結果，可發現不透水圓管1切線速度u_imp亦是隨時間作週期性震盪，max |u_imp|發生在t/T=0及t/T=0.5，以此決定選取t/T=0作為本測試案例之模擬時間。

(二)震盪流角速度

若選取t/T=0作為模擬時間之代表值，控制方程式與邊界條件即簡化成穩態均勻入射流流況之情形，則震盪角速度ω將會在化簡過程中隨三角函數被消去，故在t/T=0之前提下，ω的改變將完全不影響模擬結果。

(三)流體區入射流流速

以測試案例改變入射流流速U=0.5、1.0、2.0、3.0(m/s)，其結果如圖16及圖17所示，圖16可觀察到透水層表面的切線速度u_p正比於入射流流速，圖17則顯示入射流流速的改變不影響減緩效益，故往後模擬選取U=1(m/s)作為入射流流速。

(四)圓管半徑

參閱圖2，圓管尺寸將依用途而有所不同，透水材質厚度是實際應用設計的要點之一。以測試案例改變包覆透水層2的圓管半徑R1=0.525、0.55、0.6、0.65(m)，即R1/R2=1.05、1.1、1.2、1.3，其結果如圖18與圖19所示，可知測試案例中包覆之透水層2愈厚，則切線速度u_p愈小，減緩效果愈大。

(五)圓管埋設角度

管線埋設角度(深度)將直接影響鋪設工程之成本及難易度，若將管線完全埋入底床中，管線將可以免於底床沖刷之風險，但其成本與難度也是最高，故該狀況不在本實施例討論範圍。以測試案例改變埋設角度30~150度，其結果如圖20與21所示，可知圓管埋設角度將不規則地影響切線速度u_p及透水層2減緩效果變化，測試案例在60度時有最小切線速度u_p及最大減緩效果。

(六)無因次擾動展開函數

在管線鋪設工作前，必定會先進行底床環境調查，並由底床土壤特性決定適當的施工方法，故代表底床透水性之無因次參數在未來工程設計中將具有重要的參考價值。以測試案例改變ε₁/ε₂=10^-3~10³，其結果如圖22與圖23所示，為合適地界定比值，使用兩種不同方式計算，第一種方式為將ε₁/ε₂=10^-3~10⁰視為ε₁≪ε₂，而將ε₁/ε₂=10⁰~10³視為ε₁≫ε₂，其計算結果如圖22及圖23中的實線；第二種方式為將ε₁/ε₂=10^-3~10³皆視為ε₁ ε₂，其計算結果如圖22及圖23中的虛線，其中的實線與虛線約莫在ε₁/ε₂>10後重合，故選取ε₁/ε₂=0.07作為ε₁≪ε₂之界定數值，即得到如圖示結果。

四、最佳配置模擬結果

透水層表面的切線速度u_p為公式52。

將公式52化簡成為公式53。

為了評估何種物理參數為最佳配置，定義兩無因次參數如公式54作為評估基礎，Π最小值之參數組合為減緩管湧作用之最佳配置，其中，u_p為透水層表面的切線速度，U為該流體入射流流速，ε ₁為該第一滲透係數，ε ₂為該第二滲透係數，α為該不透水圓管之埋設角度、R₁是該透水層表面至該不透水圓管中心之距離，R₂為該不透水圓管半徑。

當R₁=R₂=0時，即為不透水圓管之建模，故可利用公式54將不透水圓管切線速度u_imp無因次化。

除了以公式54尋找最佳配置外，公式55提供另一參數來表示透水材質減緩管湧作用之效率，讓工程師能綜合評估多種因素，例如施工難易度與成本再行決定是否採用最佳配置或者其他次佳配置，χ>1表示圓管包覆透水層後可降低切線速度，其值愈大表示減緩效率愈高，其中計算速度所對應的位置如圖14所示，u_imp為不透水圓管的切線速度，u_p為透水層表面的切線速度。

影響Π之參數有四種，選取適當的數值作為接續之模擬案例將可大量減少計算量。本實施例採用R₁/R₂>1.3，其於實際工程之施行可能性不大，故選取R₁/R₂=1.05、1.1、1.2、1.3作為模擬案例。本案例是在Π與χ變化較劇烈之ε₁/ε₂區間0.01~5選取若干個數值進行案例模擬。而ε₂的數值與底床多孔介質物理參數有關，參考宋長虹(1993)整理之天然多孔底床材料特性表，選取自然界常見材料(礫石、粗砂、細砂及粉土)作為模擬案例，並利用Carman-Kozeny方程式(Kaviany，1991)計算底床之滲透係數。

其中，n為多孔介質孔隙率，d為粒徑大小(mm)，模擬案例之多孔底床特性參數列於表2。

由於埋設角度α將不規則地影響Π與χ變化，故在α=0~180度均勻選取數值進行案例模擬，整理之模擬案例各參數列於表3。分析各物理參數對管湧作用之影響後，將模擬案例中固定不變之物理參數列於表4。

五、模擬結果與討論

若將表2與表3之模擬案例參數全數組合搭配進行模擬，其計算量相當大量亦有許多不必要結果產出，為使計算達最大效率，從埋設角度α、ε ₂、ε₁/ε₂及R₁/R₂再選取特定數值作為代表值計算，雖無法確切得到最佳配置，但預期觀察到不同參數搭配對於管湧作用之影響趨勢。

參閱圖26及27，為埋設角度α=45度，Π、χ和ε ₂變化之關係圖，三條線分別代表ε₁/ε₂為0.05、0.1及0.5。由圖26可觀察到ε ₂與Π成正比關係。從圖27可知在其他參數固定時，改變ε ₂並不會影響χ。故在選取最佳配置時不需考慮ε ₂。

參閱圖28及圖29，為ε₁/ε₂=0.05，Π、χ和ε₁/ε₂變化之關係圖，折線為Π、χ以柱狀圖呈現結果。由圖中可觀察在ε₁/ε₂=0.05，即透水層2相對於多孔介質的透水底床4不透水時，其χ大部分落在1左右，表示若不透水圓管1外所包覆之透水層2的材質相對不透水，則其對減緩管湧現象之效益幾乎是零作用。

參閱圖30及圖31，為埋設角度α=90度，Π、χ和ε₁/ε₂變化之關係，實線為Π，虛線為χ。Π、χ隨ε₁/ε₂變化影響劇烈且最佳配置之位置亦是不規則，但在選取區間內Π、χ皆只有一個極值出現，可在選取最佳配置時利用該特性作為判斷依據。

參閱圖32及圖33，為α=45度，Π、χ與R₁/R₂變化之關係，折線為Π，χ以柱狀圖呈現結果，因R₁/R₂改變會影響R₁大小，故ε ₂不為固定值。在某些特定案例中，Π，χ隨R₁/R₂變化影響劇烈，且並不是增加包覆透水層2厚度即可得更佳的減緩效果，有些情況越薄的透水層2更能減緩管湧作用，故在選取最佳配置時要考慮R₁/R₂之變化。因ε ₂之改變不影響最佳配置的選取，故在以下之最佳配置模擬中將假設透水底床4的材質為粗砂進行模擬，Π，χ隨α與ε₁/ε₂變化較劇烈且不規則，故本實施例先尋找到最佳之搭配，再改變R₁/R₂得最佳配置。

本實施例中採用埋設角度α=30、45、60、75、90、105、120、135及150度，ε₁/ε₂=0.1~1.2，R₁/R₂=1.2之最佳配置模擬結果，如表5。

參閱圖34，並配合表5，本實施例依埋設角度挑選最佳ε₁/ε₂比值後，接著尋找最佳內外徑比，改變R₁/R₂=1.05、1.1、1.15、1.2、1.25及1.3，其結果當α≧60度之各角度配置皆有數百倍之減緩效率，在α=75度、R₁/R₂=1.3時的減緩效率可達兩千多倍，而減緩效果與R₁/R₂並無一定規則。

參閱圖35，一電子裝置5用以執行本實施例的抑制管湧現象之模擬方法，具有一運算單元50、一輸入單元51、一記憶單元52及一輸出單元53。其中，輸入單元51可以是(但不限於)一鍵盤，用以輸入各種模擬參數以提供運算單元50運算之用。運算單元 50載入對應如前述各種邊界值問題之方程式及邊界條件，並讀取記憶單元52預先儲存之多個物理參數以依據正規擾動法將該等物理參數換算為無因次化參數，輸出單元53可以是(但不限於)一顯示器，用以顯示運算單元50運算後的輸出結果。

綜上所述，本發明之抑制管湧現象之模擬方法針對設置於流體3中及透水底床4上且包覆透水層2的不透水圓管1之模型100，提出相關方程式的處理流程及實驗結果，除了可依據各種物理參數進行模擬以進行抑制管湧現象之模擬，並提出評估方程式以最佳化配置降低管湧現象，故確實能達成本發明之目的。

惟以上所述者，僅為本發明之實施例而已，當不能以此限定本發明實施之範圍，凡是依本發明申請專利範圍及專利說明書內容所作之簡單的等效變化與修飾，皆仍屬本發明專利涵蓋之範圍內。

S101~S106‧‧‧步驟

Claims

一種抑制管湧現象之模擬方法，應用於模擬一模型的管湧現象，該模型包括在一流體中的一透水底床、一不透水圓管，及一包覆於該不透水圓管且接觸該透水底床之透水層，該方法由一具有一運算單元及一記憶單元之電子裝置執行以下步驟：(A)該運算單元載入對應該模型的邊界值問題之方程式及邊界條件；(B)該運算單元讀取該記憶單元預先儲存之多個物理參數並依據一正規擾動法將該等物理參數換算為二個無因次化參數，該等無因次化參數包括該透水層的一第一滲透係數及該透水底床的一第二滲透係數；及(C)該運算單元判斷該第一滲透係數遠大於該第二滲透係數或該第一滲透係數遠小於該第二滲透係數時，將該透水底床及該透水層交界面的邊界條件修正成具有通量之連續邊界條件，並配合該正規擾動法以產生該模型對應之流速、流速勢及/或流體壓力之分佈。
如請求項1所述抑制管湧現象之模擬方法，其中，該流體之流體壓力P、流速勢Φ及動力黏滯係數μ之公式分別為： ,,，其中，ρ是該流體之密度，g為重力加速度，R₁是該透水層表面至該不透水圓管中心之距離。
如請求項2所述抑制管湧現象之模擬方法，其中，該第一滲透係數之公式為及該第二滲透係數之公式為，其中，k₂為該透水層之孔隙率，n₂ 為該透水層之孔隙率，k₃為該透水底床之孔隙率，n₃為該透水底床之孔隙率。
如請求項3所述抑制管湧現象之模擬方法，其中，該等方程式具有一第一方程式組，若為該第一滲透係數接近該第二滲透係數之狀態則執行該第一方程式組包括：令ε ₁=ε ₂=ε，Φ₁為該流體之速度勢，Φ₂為該透水層之速度勢，Φ₃為該透水底床之速度勢，,j=1,2,3；其中，O代表量階(order of magnitude)。
如請求項3所述抑制管湧現象之模擬方法，其中，該等方程式具有一第二方程式組，若為該第一滲透係數遠大於該第二滲透係數之狀態則執行該第二方程式組包括：Φ₁為該流體之速度勢，Φ₂為該透水層之速度勢，Φ₃為該透水底床之速度勢，,j=1,2,3；令代入，N₂、N₃為正向單位向量，得到
如請求項3所述抑制管湧現象之模擬方法，其中，該等方程式具有一第三方程式組，若為該第一滲透係數遠小於該第二滲透係數之狀態則執行該第三方程式組包括：Φ₁為該流體之速度勢，Φ₂為該透水層之速度勢，Φ₃為該透水底床之速度勢，,j=1,2,3； ,j=1,2,3；以及令代入，N₂、N₃為正向單位向量，得到
如請求項1所述抑制管湧現象之模擬方法，還包含下述步驟：(D)定義無因次參數的一第一評估公式作為最佳配置之評估基礎，該第一評估公式為，其中，u_p為透水層表面的切線速度，U為該流體入射流流速，ε ₁為該第一滲透係數，ε ₂為該第二滲透係數，α為該不透水圓管之埋設角度、R₁是該透水層表面至該不透水圓管中心之距離，R₂為該不透水圓管半徑。
如請求項1所述抑制管湧現象之模擬方法，還包含下述步驟：(E)定義無因次參數的一第二評估公式作為最佳配置之評估基礎，該第二評估公式為，其中，u_imp為不透水圓管的切線速度，u_p為透水層表面的切線速度。
一種電腦程式產品，當電腦載入該電腦程式並執行後，可完成如請求項1至8中任一請求項所述抑制管湧現象之模擬方法。