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

（中國電建集團中南勘測設計研究院有限公司，湖南 長沙 410116）

1 研究背景

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

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

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

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

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

在求解動水壓力時，多以動量方程求解，即：

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

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

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

2 虛擬面的概化

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

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

圖2 曲面概化方法示意圖

圖3 計算邏輯圖

3 計算思路

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

4 計算步驟

本計算方法的計算步驟如下：

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

步驟2：以網(wǎng)格坐標與網(wǎng)格要素（壓強、密度等）一一對應的形式導出數(shù)學模型計算結果。

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

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

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

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

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

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

步驟7：一次令i＝1，2，…，n，沿x方向對每一層網(wǎng)格做步驟4～步驟6，即可匯總到曲面在x方向投影的受力分布，即每一個法向為x方向虛擬面上的壓強差△px，i，j，k.

步驟8：對于所有法向為x方向的虛擬面，按式（6）可計算出單位y方向距離的曲面在x方向的總受力分布，即：

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

步驟9：對所有法向為x方向的虛擬面，按照式（7）可以計算出單位z方向距離的曲面在x方向總受力的分布，即：

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

步驟10：按照式（8）計算出每一個法向為x方向的虛擬面上的受力，即：

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

對所有法向為x方向的虛擬面上的受力求和，即可得到整個曲面在x方向的總受力Fx，即：

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

步驟11：沿y方向對每一層網(wǎng)格作與步驟4～步驟11相似的操作，可匯總到曲面y方向的壓強差分布，即每個法向為y方向虛擬面上的壓強差△px，i，j，k.同時還可計算每個法向為y方向的虛擬面的受力fy，i，j，k，法向為y方向的虛擬面在單位x方向的總受力分布fy，dx，法向為y方向的虛擬面在單位z方向的總受力分布fy，dz，曲面在y方向的總受力Fy.

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

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

同理，還需對步驟6的計算公式進行修正，即：

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

沿z方向對每一層網(wǎng)格采用修正后的公式進行與步驟4～步驟11相似的操作，可以匯總到曲面z方向的壓強差分布，即每一個法向為z方向的虛擬面上的壓強差△pz，i，j，k.同時，還可以計算得到每個法向為z方向的虛擬面的受力fz，i，j，k，法向為z方向的虛擬面在單位x方向的總受力分布fz，dx，法向為z方向的虛擬面在單位y方向的總受力分布fz，dy，曲面在z方向的總受力Fz.

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

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

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

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

由此，可以得到曲面的壓強分布。

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

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

5 結論

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

6 展望

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