基于數(shù)值模擬的流體中曲面受力計(jì)算分析方法研究

李 翔，黃膺翰，顏劍波，李 璜

（中國(guó)電建集團(tuán)中南勘測(cè)設(shè)計(jì)研究院有限公司，湖南 長(zhǎng)沙 410116）

1 研究背景

在工程領(lǐng)域，攔污柵、隔水帷幕、弧形閘等外觀近曲面的工程設(shè)施具有形態(tài)復(fù)雜、受力計(jì)算困難、受力分布計(jì)算量大的特點(diǎn)，現(xiàn)有的計(jì)算方法多將這些曲面的受力分為靜水壓力和動(dòng)水壓力2方面求解。在求解靜水壓力時(shí)，認(rèn)為壓力符合靜壓分布，如果曲面一側(cè)為流體，一側(cè)為大氣，則靜水壓力計(jì)算公式為：

式（1）中：ps為靜水壓強(qiáng)；ρ為流體密度；g為重力加速度；h為水深。

如果將曲面置于流體中，兩側(cè)都有流體，則認(rèn)為曲面受到的靜水壓力為兩側(cè)靜水壓力之差，即：

式（2）中：pfrontside為曲面正面的靜水壓強(qiáng)；pbackside為曲面背面的靜水壓強(qiáng)。

該方法假設(shè)靜水壓力沿垂向變化，在水平方向相等，因此，僅適用于靜水壓力水平方向差異比較小的情況。如果流體在水平方向由于鹽度、含沙量、溫度等因素分布不均，則會(huì)呈現(xiàn)出密度差異，引起靜水壓力在水平方向的變化，那時(shí)，采用該方法難以正確求解。

在求解動(dòng)水壓力時(shí)，多以動(dòng)量方程求解，即：

該方法可求解簡(jiǎn)單情況下曲面承受的總動(dòng)水壓力受力，但是，不能精細(xì)求解動(dòng)水壓力的分布情況。另外，該方法也難以求解繞流等復(fù)雜流態(tài)下曲面的動(dòng)水壓力。

隨著數(shù)值模擬計(jì)算應(yīng)用的普及，當(dāng)前部分商業(yè)數(shù)值模擬軟件也具有曲面受力計(jì)算的功能，但是，這些軟件存在以下2個(gè)缺陷：①部分商業(yè)軟件僅可輸出曲面的整體受力，不能得到曲面的受力分布；②當(dāng)曲面至于流體中，兩側(cè)均有流體時(shí)，如果部分曲面的網(wǎng)格正面為流體，背面為固體，則該網(wǎng)格正面壓力會(huì)被計(jì)算為靜水壓力與動(dòng)水壓力之和，背面壓力會(huì)被識(shí)別為無(wú)壓力。而在實(shí)際情況下，如果曲面貼向固體，會(huì)受到固體的頂托力，以平衡流體帶來(lái)的壓力。因此，這種算法與實(shí)際受力不符，當(dāng)水深較深時(shí)，使用這個(gè)方法會(huì)導(dǎo)致非常大的計(jì)算誤差。

隨著工程技術(shù)的發(fā)展，對(duì)水工建筑物的受力分析提出了更高的要求，但現(xiàn)有的曲面受力計(jì)算方法已難以滿足工程應(yīng)用的需要。針對(duì)上述技術(shù)問題，本文提出了一種能在數(shù)值模擬計(jì)算結(jié)果的基礎(chǔ)上精確求解和分析流體中復(fù)雜曲面受力分布的計(jì)算方法。該計(jì)算方法計(jì)算邏輯與編程語(yǔ)言相符，能有效解決曲面受力分布計(jì)算量大的問題。

2 虛擬面的概化

圖1 網(wǎng)格劃分示意圖

如圖1所示，在數(shù)值模擬時(shí)，模擬空間會(huì)被劃分為若干排列有序的網(wǎng)格。如果模擬網(wǎng)格為結(jié)構(gòu)化網(wǎng)格，網(wǎng)格排列為矩陣形式。空間中的每一個(gè)網(wǎng)格都會(huì)被賦予坐標(biāo)根據(jù)其坐標(biāo)（x，y，z）的排序，賦予序號(hào)編號(hào)（i，j，k）。其中，i為網(wǎng)格沿y方向的序號(hào)，j為網(wǎng)格沿x方向的序號(hào)，k為網(wǎng)格沿z方向的序號(hào)。網(wǎng)格與網(wǎng)格之間存在網(wǎng)格交界面。如圖2所示，在網(wǎng)格解析實(shí)體曲面時(shí)，模型會(huì)根據(jù)網(wǎng)格尺寸和類型將實(shí)體的曲面（紅色曲線）概化為若干平面組成的近似于曲面的復(fù)雜形態(tài)，而這些平面即為特殊的網(wǎng)格交界面（以下稱“虛擬面”）。當(dāng)整個(gè)模型在迭代求解時(shí)，沒有被概化為曲面的網(wǎng)格交界面允許水流通過，而被概化為曲面的虛擬面不允許水流通過。也就是說(shuō)，在模型計(jì)算時(shí)，這些虛擬面就是曲面。因此，求出這些虛擬面的受力便能求解出曲面的受力。

圖2 曲面概化方法示意圖

圖3 計(jì)算邏輯圖

3 計(jì)算思路

對(duì)于結(jié)構(gòu)化網(wǎng)格，曲面周圍流場(chǎng)的數(shù)值模擬計(jì)算結(jié)果主要為位置與數(shù)據(jù)的關(guān)系。對(duì)于每一個(gè)網(wǎng)格，都有表示位置的坐標(biāo)和對(duì)應(yīng)的比如流速、壓強(qiáng)等數(shù)據(jù)。由于虛擬面上的壓強(qiáng)差、流速梯度等參數(shù)具有較大的區(qū)別，可作為虛擬面的特征之一。因此，按照一定的順序，根據(jù)曲面周圍的網(wǎng)格計(jì)算出所有網(wǎng)格的網(wǎng)格交界面后，可根據(jù)計(jì)算結(jié)果識(shí)別出虛擬面，對(duì)其進(jìn)行更精細(xì)的受力計(jì)算分析。此即為本方法的計(jì)算思路，詳細(xì)計(jì)算邏輯如圖3所示。

4 計(jì)算步驟

本計(jì)算方法的計(jì)算步驟如下：

步驟1：采用數(shù)學(xué)模型計(jì)算曲面周圍流體的流場(chǎng)、壓力場(chǎng)和密度場(chǎng)。

步驟2：以網(wǎng)格坐標(biāo)與網(wǎng)格要素（壓強(qiáng)、密度等）一一對(duì)應(yīng)的形式導(dǎo)出數(shù)學(xué)模型計(jì)算結(jié)果。

步驟3：將x方向、y方向網(wǎng)格點(diǎn)對(duì)應(yīng)壓力、密度、網(wǎng)格類型按照坐標(biāo)升序排列，垂向網(wǎng)格點(diǎn)對(duì)應(yīng)壓力、密度、網(wǎng)格類型按坐標(biāo)降序排列。

步驟4：對(duì)于任意第i層網(wǎng)格，將該層的每一個(gè)網(wǎng)格與相鄰的第i＋1層對(duì)應(yīng)網(wǎng)格的壓強(qiáng)值相減，即可得到第i層網(wǎng)格第i＋1層網(wǎng)格所有網(wǎng)格交界面（不論網(wǎng)格交界面是否為虛擬面，均要進(jìn)行計(jì)算）在y方向的壓強(qiáng)差△px，i，j，k，即：

式（4）中：△px，i，j，k為編號(hào)為（i，j，k）的網(wǎng)格在 x 方向網(wǎng)格交界面上的壓強(qiáng)差；pi，j，k為編號(hào)為（i，j，k）的網(wǎng)格在網(wǎng)格中心處的壓強(qiáng)值；pi＋1，j，k為編號(hào)為（i＋1，j，k）的網(wǎng)格在網(wǎng)格中心處的壓強(qiáng)值。

步驟5：識(shí)別第i層與第i＋1層每一個(gè)網(wǎng)格的網(wǎng)格類型，以及對(duì)應(yīng)網(wǎng)格交界面的壓強(qiáng)差△px，i，j，k，判定該交界面是否為幕墻虛擬面。如果交界面是曲面概化形成的虛擬面，則記錄下該交界面的位置與壓強(qiáng)差。

步驟6：檢查是否曲面概化時(shí)，出現(xiàn)圖2所示的灰色折線所示情況，即沿計(jì)算方向的連續(xù)2個(gè)面均為虛擬面。因此，將第i層網(wǎng)格與第i＋2層網(wǎng)格對(duì)應(yīng)的壓強(qiáng)值相減，并記錄下該交界面的位置與壓強(qiáng)差，該步驟計(jì)算公式為；

式（5）中：pi＋2，j，k為編號(hào)為（i＋2，j，k）的網(wǎng)格在網(wǎng)格中心處的壓強(qiáng)值。

步驟7：一次令i＝1，2，…，n，沿x方向?qū)γ恳粚泳W(wǎng)格做步驟4～步驟6，即可匯總到曲面在x方向投影的受力分布，即每一個(gè)法向?yàn)閤方向虛擬面上的壓強(qiáng)差△px，i，j，k.

步驟8：對(duì)于所有法向?yàn)閤方向的虛擬面，按式（6）可計(jì)算出單位y方向距離的曲面在x方向的總受力分布，即：

式（6）中：fx，dy為法向?yàn)閤方向的虛擬面在單位y方向距離的曲面在x方向的總受力分布；nz為計(jì)算區(qū)域在z方向的網(wǎng)格總數(shù)。

步驟9：對(duì)所有法向?yàn)閤方向的虛擬面，按照式（7）可以計(jì)算出單位z方向距離的曲面在x方向總受力的分布，即：

式（7）中：fx，dz為單位z方向距離的曲面在x方向的總受力分布；ny為計(jì)算區(qū)域在y方向的網(wǎng)格總數(shù)。

步驟10：按照式（8）計(jì)算出每一個(gè)法向?yàn)閤方向的虛擬面上的受力，即：

式（8）中：fx，i，j，k為編號(hào)為（i，j，k）的虛擬面在 x 方向上的受力；yi，j＋1，k為編號(hào)為（i，j＋1，k）的網(wǎng)格在網(wǎng)格中心處的 y 坐標(biāo)值；yi，j－1，k為編號(hào)為（i，j－1，k）的網(wǎng)格在網(wǎng)格中心處的 y 坐標(biāo)值；zi，j，k－1為編號(hào)為（i，j，k－1）的網(wǎng)格在網(wǎng)格中心處的 y 坐標(biāo)值；zi，j，k＋1為編號(hào)為（i，j，k＋1）的網(wǎng)格在網(wǎng)格中心處的y坐標(biāo)值。

對(duì)所有法向?yàn)閤方向的虛擬面上的受力求和，即可得到整個(gè)曲面在x方向的總受力Fx，即：

式（9）中：Fx為曲面在x方向的總受力；nz為計(jì)算區(qū)域在z方向的網(wǎng)格總數(shù)。

步驟11：沿y方向?qū)γ恳粚泳W(wǎng)格作與步驟4～步驟11相似的操作，可匯總到曲面y方向的壓強(qiáng)差分布，即每個(gè)法向?yàn)閥方向虛擬面上的壓強(qiáng)差△px，i，j，k.同時(shí)還可計(jì)算每個(gè)法向?yàn)閥方向的虛擬面的受力fy，i，j，k，法向?yàn)閥方向的虛擬面在單位x方向的總受力分布fy，dx，法向?yàn)閥方向的虛擬面在單位z方向的總受力分布fy，dz，曲面在y方向的總受力Fy.

步驟12：考慮到每一個(gè)網(wǎng)格的壓強(qiáng)值為網(wǎng)格中心處位置的壓強(qiáng)值，而非網(wǎng)格交界面上的壓強(qiáng)值，對(duì)于x方向和y方向，在一個(gè)網(wǎng)格尺度內(nèi)，網(wǎng)格中心處的壓強(qiáng)值與網(wǎng)格邊沿的壓強(qiáng)值無(wú)明顯變化，其誤差可以忽略。但是，在z方向，靜水壓力沿z方向有明顯變化，會(huì)給計(jì)算帶來(lái)明顯誤差。因此，在沿z方向?qū)γ恳粚泳W(wǎng)格作與步驟4～步驟7相似的操作時(shí)，需對(duì)步驟4的計(jì)算公式進(jìn)行如下修正：

式（10）中：△pz，i，j，k為編號(hào)為（i，j，k）的網(wǎng)格在 z 方向網(wǎng)格交界面上的壓強(qiáng)差；pi，j，k＋1為編號(hào)為（i，j，k＋1）的網(wǎng)格在網(wǎng)格中心處的壓強(qiáng)值；ρi，j，k為編號(hào)為（i，j，k）的網(wǎng)格在網(wǎng)格中心處的密度值；ρi，j，k＋1為編號(hào)為（i，j，k＋1）的網(wǎng)格在網(wǎng)格中心處的密度值；zi，j，k為編號(hào)為（i，j，k）的網(wǎng)格在網(wǎng)格中心處的 z 坐標(biāo)值；zi，j，k＋1為編號(hào)為（i，j，k＋1）的網(wǎng)格在網(wǎng)格中心處的z坐標(biāo)值。

同理，還需對(duì)步驟6的計(jì)算公式進(jìn)行修正，即：

式（11）中：pi，j，k＋2為編號(hào)為（i，j，k＋1）的網(wǎng)格在網(wǎng)格中心處的壓強(qiáng)值；ρi，j，k＋2為編號(hào)為（i，j，k＋2）的網(wǎng)格在網(wǎng)格中心處的密度值；zi，j，k＋2為編號(hào)為（i，j，k＋2）的網(wǎng)格在網(wǎng)格中心處的z坐標(biāo)值。

沿z方向?qū)γ恳粚泳W(wǎng)格采用修正后的公式進(jìn)行與步驟4～步驟11相似的操作，可以匯總到曲面z方向的壓強(qiáng)差分布，即每一個(gè)法向?yàn)閦方向的虛擬面上的壓強(qiáng)差△pz，i，j，k.同時(shí)，還可以計(jì)算得到每個(gè)法向?yàn)閦方向的虛擬面的受力fz，i，j，k，法向?yàn)閦方向的虛擬面在單位x方向的總受力分布fz，dx，法向?yàn)閦方向的虛擬面在單位y方向的總受力分布fz，dy，曲面在z方向的總受力Fz.

步驟13：對(duì)于曲面所在網(wǎng)格，如果網(wǎng)格僅有一個(gè)方向的虛擬面（例如x方向，其他方向與之類似），則該曲面在該網(wǎng)格處的壓強(qiáng)為此虛擬面的壓強(qiáng)，即：

式（12）中：pi，j，k為曲面在編號(hào)（i，j，k）網(wǎng)格處的壓強(qiáng)。

如果網(wǎng)格含有2個(gè)虛擬面，則曲面在該網(wǎng)格處的壓強(qiáng)為2個(gè)虛擬面（例如x方向和y方向）的壓強(qiáng)之和，即：

如果網(wǎng)格含有x方向、y方向、z方向的虛擬面，則曲面在該網(wǎng)格處的壓強(qiáng)為3個(gè)虛擬面（例如x方向和y方向）的壓強(qiáng)之和，即：

由此，可以得到曲面的壓強(qiáng)分布。

步驟14：對(duì)曲面在x方向、y方向、z方向的總受力進(jìn)行求和，可得到曲面的總受力為：

式（13）中：Fsurface為曲面總受力。

5 結(jié)論

針對(duì)當(dāng)前工程領(lǐng)域內(nèi)流體中曲面受力計(jì)算分析方法存在的不足，本文提出了一種能在數(shù)值模擬計(jì)算結(jié)果的基礎(chǔ)上精確求解和分析流體中復(fù)雜曲面受力分布的計(jì)算方法。該計(jì)算方法計(jì)算邏輯與編程語(yǔ)言相符，能有效解決曲面受力分布計(jì)算量大的問題。

6 展望

受文章篇幅限制，在此僅提出流體中曲面受力計(jì)算分析方法，該方法在工程上的應(yīng)用將另行他文介紹。