為了輕松生成隨機幾何表面，COMSOL Multiphysics? 軟件內置了一組功能強大的函數和算子，例如用于均勻和高斯隨機分布的函數以及非常有用的求和算子。在這篇文章中，我們將展示如何使用相當于“一個線性”表達式生成一個隨機表面，并對決定表面粗糙度性質的空間頻率分量進行詳細控制。

表征表面粗糙度

有許多方法可以表征粗糙表面。一種方法是使用它的近似分形維數，表面的分形維數值在 2 ~ 3 之間。分形維數為 2 的表面是一個普通的、近乎絕對光滑的表面; 分形維數為 2.5 表示相當崎嶇的表面; 分形維數接近 3 表示接近“3D 空間填充”的內容。相應地，分形維數為 1 的曲線幾乎在所有地方都是平滑的，1.5 表示一條非常崎嶇的線，接近 2表示接近“2D 空間填充”的內容。

曲線的分形維數值的范圍從 1(左)到大約 1.2(中)再到 1.6(右)。

使用分形維數度量可能是一個有用的近似值，但需要記住，真實的表面在超過幾個數量級的尺度上不是分形的。真實表面具有空間頻率“截止”點，這是因為它們的尺寸有限，并且當將它的細微部分“放大”后，其結構特征仍與原來的一樣。

表征表面粗糙度的另一種方法是它的空間頻率含量。這可以變成一種建設性方法，通過使用類似于傅里葉級數擴展的三角函數之和來合成表面數據。這種總和中的每個項都表示在空間中振蕩的某個頻率。這也是我們在本文中將使用的方法。在介紹三角級數之前，我們快速回顧一下空間頻率和基本波形的概念。

空間頻率

在物理學中，隨時間變化的振蕩頻率一般以數學表達式的形式出現，例如

其中頻率 f 的單位是 1/s，也稱為赫茲或 Hz。

通過空間的振蕩具有相應的空間頻率，如下面的表達式所示，我們只需將時間變量 t 替換為空間變量x，將時間頻率 f 替換為空間頻率 v。

其中空間頻率的 SI 單位為 1/m。

空間頻率通常由波數 k= 2πv 表示。

一個相關的量是波長，它與頻率和波長   有關，如下式所示：

空間可能有多個維度，因此可能存在多個空間頻率。在 2D 中，使用笛卡爾坐標表示：

其中，波矢量   并且 

波矢量   表示波的方向。

基本波

粗糙的表面 f（x，y）可以看作是由許多基本波組成的

其中，φ 是相位角。

由于，相位角   也可以表示正弦函數。

對于完全隨機的表面，應該認為相位角 φ 可以取任何值，例如，0 到 π 或 –π/2 到 π/2
之間。當為隨機表面合成基本波時，我們可以在這樣的長度為 π 的區間內均勻隨機分布中挑選 φ，因為我們允許表達式   取 -1 到 +1 之間的所有可能值。請注意，如果選擇大小大于 π 的間隔，則可能會出現終點或環繞效應。這是由于根據  ，余弦函數是它自己的鏡像，步長為 π。

為了獲得可用于模擬的有效表述，我們只允許一組離散的空間頻率：

νx = m， νy = n

其中，m 和 n 是整數。

考慮一個由以下形式的基本波組成的表面：

通過讓 m 和 n 以相等的概率取正值和負值，應該能夠得到一種方法來合成一個沒有首選振蕩方向的表面。

請注意，如果采用這種方式，每個波方向將表示兩次。例如，方向(-2,-3)與(2,3)相同;(2,-1)與(-2,1)相同；等等。

如果我們允許空間頻率 m 和 n 分別取最大整數值 M 和 N，那么對應于以下位置的高頻截止：

由于也允許負值，因此在以下位置存在負截止值：

在 x 方向上具有空間頻率截止   意味著可以表示的最短波長是  ，對于y方向也是如此，即  。

基本波的相關振幅

每個基本波都有一個相關的振幅，因此每個組成波分量具有以下形式：

最終表面將是這些波分量的總和：

最簡單的振幅選擇是選擇一個來自均勻分布或高斯分布的系數 Amn。但是，事實證明，這不會生成看起來特別自然的表面。在自然界中，不同的過程，如磨損和侵蝕，使得慢速振蕩比快速振蕩更有可能具有更大的振幅。在離散情況下，這對應于根據某種分布逐漸變細的振幅：

其中頻譜指數 β 表示較高頻率衰減的速度。根據 分形圖像科學（參考文獻1），譜指數與表面的分形維數相關，但僅適用于覆蓋任意高頻的無限系列波，并且僅適用于指數的某些范圍。在實踐中，合成表面的振幅 a(m，n)將使用有限數量的頻率乘以具有高斯分布的隨機函數 g(m，n)生成：

a(m,n) = g(m,n)h(m,n)

選擇高斯分布或正態分布獲得幅度平滑但隨機的變化，而幅度沒有限制。

相位角 φ 將從函數 u 采樣，函數 u 在 –π/2 和 π/2 之間具有均勻隨機分布：

φ(m,n) = u(m,n)

總結一下

為了表示粗糙表面，我們希望使用以下二重求和：

其中，x 和 y 是空間坐標;m 和 n 是空間頻率;a（m，n） 是振幅;φ（m，n）是相位角。該表達式類似于截斷的傅里葉級數。盡管該級數是用余弦函數表示的，但由于角度求和規則，相位角使得該求和可以表示一個非常一般的三角級數：

確定周期性

由于函數 f*(x，y) 的定義， 它將是周期性的。為了獲得自然的表面，應該通過讓 x 和 y 在某些有限值之間變化來“切出”適當的一小部分;否則，合成數據的周期性就會很明顯。這些值應該是什么？

整體周期性將由最慢振蕩決定，最慢振蕩分別對應于 x 方向和 y 方向的空間頻率 m= 1 或 n= 1。這使每個方向上的周期長度為 1。

我們可以在矩形 [a, a + 1] × [b, b + 1] 或更小的矩形上生成表面，來“避免”周期性。

在 COMSOL Multiphysics? 中定義參數和隨機函數

關于在 COMSOL Multiphysics中的實現，首先根據下圖定義幾個空間頻率分辨率和頻譜指數參數：

振幅的生成將需要在兩個變量中具有高斯分布的隨機函數。COMSOL 軟件的全局定義節點提供了此功能：

在這里，標簽 和函數名稱 已分別更改為高斯隨機 和 g1。此外，變元數 被設置為 2 而不是默認值 1，分布類型設置為“正態”分布，相當于于正態分布或高斯分布。

對于相位角，以類似的方式，我們需要在 –π/2 和 π/2 的區間內有一個均勻的隨機函數：

標簽 更改為均勻隨機，函數名稱 更改為 u1，變元數 更改為2，范圍 更改為 pi。

我們可以選擇使用隨機種子，以便在每次使用相同的輸入參數時獲得相同的表面。

定義參數化曲面

下一步是使用一個相當長的 z 坐標表達式，在幾何 下添加一個參數化曲面 節點，如下所示：

0.01*sum(sum(if((m!=0)||(n!=0),((m^2+n^2)^(-b/2))g1(m,n)*cos(2*pi(m*s1+n*s2)+u1(m,n)),0),m,-N,N),n,-N,N)

其中，x=s1 和 y=s2 在 0 和 1 之間變化。

系數 0.01 用于在縮放 z 方向上的數據。或者，該比例因子可以被應用到振幅系數中去。

參數化曲面幾何特征用于生成合成的隨機表面。

請注意，每當更新參數化曲面的任何參數或表達式時，都需要單擊設置窗口的高級設置 部分中的使用更新的功能重建 按鈕。

此表達式是整數參數 m 和 n 在 –N 到 N 之間的二重和。如果我們將其與前面的數學討論進行比較，可以看到已經設置了 M=N，從而產生了一個方形表面補丁。m 和 n 同時為零的項對應于不需要的“DC”項，并通過 if 語句從總和中消除。

sum()算子的語法如下：

sum(expr,index,lower,upper)

對于所有從低 到高 的索引，它計算通用表達式 expr 的總和。

if()算子的語法如下：

if(cond,expr1,expr2)

條件表達式 cond 的計算結果為 expr1 或 expr2，具體取決于條件的值。

在這個示例中，通過將最大節數 設置為 100（默認值為 20），提高了參數化表面的分辨率。此外，相對容差 被放寬到了 1e-3（默認值為 1e-6）。參數表面的底層表示基于非均勻有理 B 樣條曲線 （NURBS）。更多的節點對應于 NURBS 表示的更精細分辨率。因為我們并不過分關注本例中生成表面的近似精度，所以放寬了容差。

通過生成網格，我們可以獲得有用的表面可視化圖，如下所示。

網格隨機表面。

請注意，N = 20 意味著假設使用 SI 單位，最快的振蕩為 1/20 = 0.05 m。x 和 y 方向的周期性可以分別在 x= 0，x= 1 和 y= 0，y= 1 的平行于 y 軸和 x 軸的曲線處看到。

為了更清楚地看到周期性，我們可以在正方形 [0,2] × [0,2] 上繪制表面：

表面在正方形 [0，2] × [0，2] 上的周期性。表面高度由顏色表示。

在正方形 [0，1] × [0，1] 上生成的表面，方法是從左上角圖像順時針方向疊加 20 個頻率分量，幅度譜指數 β = 0.5、β = 1.0、β = 1.5 和 β = 1.8。表面高度由顏色表示。

在分析中使用表面數據

在 COMSOL Multiphysics 中，這種類型的隨機生成表面可用于任何類型的物理仿真環境，包括電磁學、結構力學、聲學、流體、熱或化學分析。二重和的表達式不僅限于幾何建模，還可用于材料數據、方程系數、邊界條件等。使用這個方法，可以在循環中使用大量表面來收集結果的統計信息。

通過將二重和推廣為三重和，可以合成 3D 非均勻材料數據。但是，在為 3D 模擬執行三重和時，必須做好耗時和占內存的計算準備。

基于合成生成的裂縫孔徑數據的裂縫流動模擬。巖石裂隙流模型是 COMSOL Multiphysics案例庫中 的一個案例模型。

基于文中描述的參數化表面，對兩個具有材料界面的 1 厘米大小的金屬塊進行通用熱膨脹分析。底部材料板是鋁，頂部材料板是鋼。可視化圖顯示了材料界面和鋁板表面的 von Mises 應力。

與離散余弦和傅里葉變換的關系

求和

類似于離散余弦變換或離散傅里葉變換的實部：

其中，下標 c 表示復數，x 和 y 現在采用離散值。這里，相位角信息以復數傅里葉系數編碼。

根據離散傅里葉變換的定義，我們可以執行索引的移位，用于生成以下我們更熟悉的形式：

或使用離散值：

更常見的是，離散傅里葉變換的索引如下：

其中，

請注意，為了生成實值數據，傅里葉系數需要滿足共軛對稱關系，以消除正弦函數的虛值貢獻。使用余弦函數的總和（即余弦變換）可以避免這個問題。

生成大量傅里葉系數的快速方法是使用快速余弦變換(FCT)或快速傅里葉變換(FFT)。這可以在另一個程序中完成，然后作為插值表導入到 COMSOL Desktop? 用戶界面中。上述三角插值方法計算較慢，但優點是可以直接在非結構化網格上使用，并且只需在用戶界面中細化網格就可以可自動細化。

一維和圓柱體示例

最后，我們來看看在 COMSOL Multiphysics 中由隨機表面生成的幾個有趣的特殊情況，包括曲線和圓柱體。

隨機曲線

在 2D 仿真中，可以使用以下表達式生成隨機曲線：

0.01*sum(if((m!=0),((m^2)^(-b/2))*g1(m)*cos(2*pi*m*s+u1(m)),0),m,-N,N)

其中， g1 和 u1 是一維隨機函數。

請注意，在生成曲線時，與“相同隨機水平”的表面相比，譜指數的值較低。

譜指數為 0.8 的隨機曲線。

隨機極性曲線

可以在極坐標中生成一條表示圓的隨機偏差的隨機曲線：

x=cos(2*pi*s)*(1+0.1*sum(if((m!=0),((m^2)^(-b/2))*g1(m)*cos(2*pi*m*s+u1(m)),0),m,-N,N))

y=sin(2*pi*s)*(1+0.1*sum(if((m!=0),((m^2)^(-b/2))*g1(m)*cos(2*pi*m*s+u1(m)),0),m,-N,N))

對應于二維極坐標中的參數曲線：

譜指數為 0.8 的隨機極性曲線。

隨機圓柱體

可以使用參數化曲面生成 3D 中的隨機圓柱體，參數如下：

x=cos(2*pi*s1)*(1+0.1*sum(sum(if((m!=0)||(n!=0),((m^2+n^2)^(-b/2))*g1(m,n)*cos(2*pi*(m*s1+n*s2)+u1(m,n)),0),m,-N,N),n,-N,N))

y=sin(2*pi*s1)*(1+0.1*sum(sum(if((m!=0)||(n!=0),((m^2+n^2)^(-b/2))*g1(m,n)*cos(2*pi*(m*s1+n*s2)+u1(m,n)),0),m,-N,N),n,-N,N))

z=s2*2*pi

其中，參數 s1 和 s2 在 0 和 1 之間變化。

這對應于圓柱坐標中的參數化曲面：

這種單一隨機圓柱體代表一種自相交表面，這在 COMSOL Multiphysics 中是不允許的。我們可以通過創建對應于參數 s1 的四個表面補丁來輕松解決此問題，這些表面補丁的范圍在 0 ~ 0.25、0.25 ~ 0.5、0.5~ 0.75 和 0.75 ~ 1.0 之間。一個這樣的補丁對應于大小為   的極角跨度。

使用極坐標的隨機管狀表面。

本文來自： COMSOL 博客