如何在 COMSOL Multiphysics® 中生成随机表面

2017年 6月 2日

为了轻松生成随机几何表面,COMSOL Multiphysics® 软件内置了一组功能强大的函数和算子,例如用于均匀和高斯随机分布的函数以及非常有用的求和 算子。在这篇博客中,我们将展示如何使用相当于“一个线性”表达式生成一个随机表面,并对决定表面粗糙度性质的空间频率分量进行详细控制。

表征表面粗糙度

有许多方法可以表征粗糙表面。一种方法是使用它的近似分形维数,表面的分形维数值在 2 ~ 3 之间。分形维数为 2 的表面是一个普通的、近乎绝对光滑的表面; 分形维数为 2.5 表示相当崎岖的表面; 分形维数接近 3 表示接近“3D 空间填充”的内容。相应地,分形维数为 1 的曲线几乎在所有地方都是平滑的,1.5 表示一条非常崎岖的线,接近 2表示接近“2D 空间填充”的内容。

分形维数值为1的曲线。
分形维数值为1.2的曲线。
分形维数值为1.6的曲线。

曲线的分形维数值的范围从 1(左)到大约 1.2(中)再到 1.6(右)。

使用分形维数度量可能是一个有用的近似值,但需要记住,真实的表面在超过几个数量级的尺度上不是分形的。真实表面具有空间频率“截止”点,这是因为它们的尺寸有限,并且当将它的细微部分“放大”后,其结构特征仍与原来的一样。

表征表面粗糙度的另一种方法是它的空间频率含量。这可以变成一种建设性方法,通过使用类似于傅里叶级数扩展的三角函数之和来合成表面数据。这种总和中的每个项都表示在空间中振荡的某个频率。这也是我们在本文中将使用的方法。在介绍三角级数之前,我们快速回顾一下空间频率和基本波形的概念。

空间频率

在物理学中,随时间变化的振荡频率一般以数学表达式的形式出现,例如

cos(2\pi f t)

其中频率 f 的单位是 1/s,也称为赫兹或 Hz。

通过空间的振荡具有相应的空间频率,如下面的表达式所示,我们只需将时间变量 t 替换为空间变量x,将时间频率 f 替换为空间频率 v

cos(2\pi \nu x)

其中空间频率的 SI 单位为 1/m

空间频率通常由波数 k= 2πv 表示。

一个相关的量是波长,它与频率和波长 \lambda=\frac{1}{\nu} 有关,如下式所示:

k=2\pi \nu=\frac{2\pi}{\lambda}

空间可能有多个维度,因此可能存在多个空间频率。在 2D 中,使用笛卡尔坐标表示:

cos(2\pi (\nu_x x + \nu_y y))=cos(\bf{k} \cdot \bf{x})

其中,波矢量 \bf{k}=(k_x,k_y)=(2\pi \nu_x,2\pi\nu_y) 并且 \bf{x}=(x,y)

波矢量 \bf{k} 表示波的方向。

基本波

粗糙的表面 fxy)可以看作是由许多基本波组成的

cos(\bf{k} \cdot \bf{x}+\phi)

其中,φ 是相位角。

由于,相位角 sin(\theta)=cos(\pi/2-\theta) 也可以表示正弦函数。

对于完全随机的表面,应该认为相位角 φ 可以取任何值,例如,0 到 π 或 –π/2 到 π/2
之间。当为随机表面合成基本波时,我们可以在这样的长度为 π 的区间内均匀随机分布中挑选 φ,因为我们允许表达式 cos(\phi) 取 -1 到 +1 之间的所有可能值。请注意,如果选择大小大于 π 的间隔,则可能会出现终点或环绕效应。这是由于根据 cos(\pi\theta)=-cos(\theta),余弦函数是它自己的镜像,步长为 π

为了获得可用于模拟的有效表述,我们只允许一组离散的空间频率:

νx = m νy = n

其中,mn 是整数。

考虑一个由以下形式的基本波组成的表面:

cos(\bf{k}_{mn} \cdot \bf{x}+\phi)= cos(2 \pi (mx+ny)+\phi) , \bf{k}_{mn}=2\pi(m,n)

通过让 mn 以相等的概率取正值和负值,应该能够得到一种方法来合成一个没有首选振荡方向的表面。

请注意,如果采用这种方式,每个波方向将表示两次。例如,方向(-2,-3)与(2,3)相同;(2,-1)与(-2,1)相同;等等。

如果我们允许空间频率 mn 分别取最大整数值 MN,那么对应于以下位置的高频截止:

\nu_{xmax}=M, \nu_{ymax}=N

由于也允许负值,因此在以下位置存在负截止值:

\nu_{xmin}=-M, \nu_{ymin}=-N

x 方向上具有空间频率截止 \nu_{xmax}=M 意味着可以表示的最短波长是 \lambda_{xmin}=\frac{1} {M},对于y方向也是如此,即 \lambda_{ymin}=\frac{1}{N}

基本波的相关振幅

每个基本波都有一个相关的振幅,因此每个组成波分量具有以下形式:

A_{mn}cos(\bf{k}_{mn} \cdot \bf{x}+\phi)

最终表面将是这些波分量的总和:

f(\bf{x})=\sum_{m,n}A_{mn}cos(\bf{k}_{mn} \cdot \bf{x}+\phi)

最简单的振幅选择是选择一个来自均匀分布或高斯分布的系数 Amn。但是,事实证明,这不会生成看起来特别自然的表面。在自然界中,不同的过程,如磨损和侵蚀,使得慢速振荡比快速振荡更有可能具有更大的振幅。在离散情况下,这对应于根据某种分布逐渐变细的振幅:

A_{mn} =a(m,n) \sim h(m,n)=\frac{1}{\vert m^2+n^2\vert^{\beta}}=\frac{1}{(m^2+n^2)^{\frac{\beta}{2}}}

其中频谱指数 β 表示较高频率衰减的速度。根据 分形图像科学参考文献1),谱指数与表面的分形维数相关,但仅适用于覆盖任意高频的无限系列波,并且仅适用于指数的某些范围。在实践中,合成表面的振幅 a(m,n)将使用有限数量的频率乘以具有高斯分布的随机函数 g(mn)生成:

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

选择高斯分布或正态分布获得幅度平滑但随机的变化,而幅度没有限制。

相位角 φ 将从函数 u 采样,函数 u 在 –π/2 和 π/2 之间具有均匀随机分布:

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

总结一下

为了表示粗糙表面,我们希望使用以下二重求和:

f(x,y)=\sum_{m=-M}^{M}\sum_{n=-N}^{N} a(m,n) cos(2 \pi(mx+ny)+\phi(m,n))

其中,xy 是空间坐标;mn 是空间频率;amn) 是振幅;φmn)是相位角。该表达式类似于截断的傅里叶级数。尽管该级数是用余弦函数表示的,但由于角度求和规则,相位角使得该求和可以表示一个非常一般的三角级数:

cos(\alpha+\beta)=cos(\alpha)cos(\beta)-sin(\alpha)sin(\beta)

确定周期性

由于函数 f*(xy) 的定义, 它将是周期性的。为了获得自然的表面,应该通过让 xy 在某些有限值之间变化来“切出”适当的一小部分;否则,合成数据的周期性就会很明显。这些值应该是什么?

整体周期性将由最慢振荡决定,最慢振荡分别对应于 x 方向和 y 方向的空间频率 m= 1 或 n= 1。这使每个方向上的周期长度为 1。

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

在 COMSOL Multiphysics® 中定义参数和随机函数

关于在 COMSOL Multiphysics中的实现,首先根据下图定义几个空间频率分辨率和频谱指数参数:

COMSOL Multiphysics 的参数设置的屏幕截图。

振幅的生成将需要在两个变量中具有高斯分布的随机函数。COMSOL 软件的全局定义 节点提供了此功能:

显示了高斯随机节点的COMSOL Multiphysics 模型开发器屏幕截图。

在这里,标签函数名称 已分别更改为高斯随机g1。此外,变元数 被设置为 2 而不是默认值 1,分布类型设置为“正态”分布,相当于于正态分布或高斯分布。

对于相位角,以类似的方式,我们需要在 –π/2 和 π/2 的区间内有一个均匀的随机函数:

显示了均匀随机设置的 COMSOL Multiphysics 屏幕截图。

标签 更改为均匀随机函数名称 更改为 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=s1y=s2 在 0 和 1 之间变化。

系数 0.01 用于在缩放 z 方向上的数据。或者,该比例因子可以被应用到振幅系数中去。

使用参数化曲面几何特征生成合成的随机表面。
参数化曲面几何特征用于生成合成的随机表面。

请注意,每当更新参数化曲面的任何参数或表达式时,都需要单击设置窗口的高级设置 部分中的使用更新的功能重建 按钮。

此表达式是整数参数 mn 在 –N 到 N 之间的二重和。如果我们将其与前面的数学讨论进行比较,可以看到已经设置了 M=N,从而产生了一个方形表面补丁。mn 同时为零的项对应于不需要的“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 表示的更精细分辨率。因为我们并不过分关注本例中生成表面的近似精度,所以放宽了容差。

通过生成网格,我们可以获得有用的表面可视化图,如下所示。

COMSOL Multiphysics 中随机曲面的网格。
网格随机表面。

请注意,N = 20 意味着假设使用 SI 单位,最快的振荡为 1/20 = 0.05 m。xy 方向的周期性可以分别在 x= 0,x= 1 和 y= 0,y= 1 的平行于 y 轴和 x 轴的曲线处看到。

为了更清楚地看到周期性,我们可以在正方形 [0,2] × [0,2] 上绘制表面:

A surface plot showing the periodicity of the square's surface.
表面在正方形 [02] × [02] 上的周期性。表面高度由颜色表示。

A surface plot for beta equals 0.5.
A surface plot for beta equals 1.
A surface plot for beta equals 1.8.
A surface plot for beta equals 1.5.

在正方形 [0,1] × [0,1] 上生成的表面,方法是从左上角图像顺时针方向叠加 20 个频率分量,幅度谱指数 β = 0.5、β = 1.0、β = 1.5 和 β = 1.8。表面高度由颜色表示。

在分析中使用表面数据

在 COMSOL Multiphysics 中,这种类型的随机生成表面可用于任何类型的物理仿真环境,包括电磁学、结构力学、声学、流体、热或化学分析。二重和的表达式不仅限于几何建模,还可用于材料数据、方程系数、边界条件等。使用这个方法,可以在循环中使用大量表面来收集结果的统计信息。

通过将二重和推广为三重和,可以合成 3D 非均匀材料数据。但是,在为 3D 模拟执行三重和时,必须做好耗时和占内存的计算准备。

基于合成裂缝孔径数据的示例模拟。
基于合成生成的裂缝孔径数据的裂缝流动模拟。岩石裂隙流模型 COMSOL Multiphysics案例库中 的一个案例模型。

基于此博客文章中描述的参数化曲面的模型。
基于本博客中描述的参数化表面,对两个具有材料界面的 1 厘米大小的金属块进行通用热膨胀分析。底部材料板是铝,顶部材料板是钢。可视化图显示了材料界面和铝板表面的 von Mises 应力。

与离散余弦和傅里叶变换的关系

求和

f(x,y)=\sum_{m=-M}^{M} \sum_{n=-N}^{N} a(m,n) cos(2 \pi(mx+ny)+\phi(m,n))

类似于离散余弦变换或离散傅里叶变换的实部:

f_c(x,y)=\sum_{m=-M}^{M} \sum_{n=-N}^{N} F_c(m,n)e^{i(2 \pi(mx+ny))}

其中,下标 c 表示复数,xy 现在采用离散值。这里,相位角信息以复数傅里叶系数编码。

根据离散傅里叶变换的定义,我们可以执行索引的移位,用于生成以下我们更熟悉的形式:

f_c(x,y)=\sum_{m=0}^{2M} \sum_{n=0}^{2N} F_c(m,n)e^{i(2 \pi(mx+ny))}

或使用离散值:

f_c(k,l)=\sum_{m=0}^{2M} \sum_{n=0}^{2N} F_c(m,n)e^{i(2 \pi(m \frac{k}{2M+1}+n \frac{l}{2N+1}))}

更常见的是,离散傅里叶变换的索引如下:

f_c(k,l)=\sum_{m=0}^{\mathfrak{M}-1} \sum_{n=0}^{\mathfrak{N}-1} F_c(m,n)e^{i(2 \pi(m \frac{k}{\mathfrak{M}}+n \frac{l}{\mathfrak{N}}))}

其中,

\mathfrak{M}=2M+1, \mathfrak{N}=2N+1.

请注意,为了生成实值数据,傅里叶系数需要满足共轭对称关系,以消除正弦函数的虚值贡献。使用余弦函数的总和(即余弦变换)可以避免这个问题。

生成大量傅里叶系数的快速方法是使用快速余弦变换(FCT)或快速傅里叶变换(FFT)。这可以在另一个程序中完成,然后作为插值表导入到 COMSOL Desktop® 用户界面中。上述三角插值方法计算较慢,但优点是可以直接在非结构化网格上使用,并且只需在用户界面中细化网格就可以可自动细化。

有关使用 FFT 合成表面的说明,请参见参考文献1

一维和圆柱体示例

最后,我们来看看在 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 是一维随机函数。

请注意,在生成曲线时,与“相同随机水平”的表面相比,谱指数的值较低。

COMSOL Multiphysics显示随机曲线的屏幕截图
谱指数为 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))

对应于二维极坐标中的参数曲线:

x=r(\phi) cos(\phi)
y=r(\phi) sin(\phi)

COMSOL Multiphysics 的屏幕截图显示了一个随机极性曲线。
谱指数为 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

其中,参数 s1s2 在 0 和 1 之间变化。

这对应于圆柱坐标中的参数化曲面:

x=r(\phi,z) cos(\phi)
y=r(\phi,z) sin(\phi)
z=z

这种单一随机圆柱体代表一种自相交表面,这在 COMSOL Multiphysics 中是不允许的。我们可以通过创建对应于参数 s1 的四个表面补丁来轻松解决此问题,这些表面补丁的范围在 0 ~ 0.25、0.25 ~ 0.5、0.5~ 0.75 和 0.75 ~ 1.0 之间。一个这样的补丁对应于大小为 \frac{\pi} {2} 的极角跨度。

COMSOL Multiphysics 的屏幕截图,显示随机管状表面。
使用极坐标的随机管状表面。

COMSOL Multiphysics® 5.5 版本中的新零件

从 5.5 版本开始,COMSOL Multiphysics 的零件库扩展了几个新零件,包括用于创建随机平面的零件。

零件库中随机平面零件的屏幕截图。
COMSOL Multiphysics 中零件库中 随机平面零件。

参考文献

  • The Science of Fractal Images, Editors: Peitgen, Heinz-Otto, Saupe, Dietmar. Eds.

评论 (20)

正在加载...
涛铭 卢
涛铭 卢
2023-03-27

工程师您好,该方法只能产生在[a, a + 1] × [b, b + 1]上的随机表面,请问应该如何避免周期性产生更大区域的随机表面呢?例如在[0,10]x[0,10]区域上产生非周期的随机表面呢?

Bjorn Sjodin
Bjorn Sjodin
2023-03-27 COMSOL 员工

The periodicity is intrinsic and cannot be avoided. However, this is not a problem. Instead, to get a larger patch you just use a unit patch and scale it up. In COMSOL using a Scale operation. There are no drawbacks with doing so. You may see that you now need additional randomness. The way to do this is to simply increase the number of frequencies: higher values of N and M. You may also need to adjust the scaling constant (0.01 above) In summary: 1) use higher values of M and N 2) scale the patch 3) adjust the scaling constant.

Bjorn

洁伯 范
洁伯 范
2023-05-15

工程师您好,这个方法的来源是什么?在论文的描述中应该如何体现或者引用?

Anran Wei
Anran Wei
2023-05-16 COMSOL 员工

有关分形图像科学与使用 FFT 合成表面的说明,请参考资料 The Science of Fractal Images, Editors: Peitgen, Heinz-Otto, Saupe, Dietmar. Eds.

洁伯 范
洁伯 范
2023-05-17

好的,谢谢

雪 李
雪 李
2023-06-22

工程师您好,建立一个长方体,长4cm,宽2cm,它的厚度随机分布在80-85mm之间,这种随机分布根据您上述所说可以实现吗?

Lei Cao
Lei Cao
2023-06-26 COMSOL 员工

李雪, 您好!

感谢您的评论。
您所说的需求可以实现。可以通过上文的方式获取随机表面,并用其分割厚度大于等于85mm的长方体,从而得到您想要的图形结果。上文的案例链接中,randomized_surface_structural.mph 模型与您的需求类似,可做参考。

如果有进一步问题,建议您联系COMSOL的技术支持团队:
在线支持中心:cn.comsol.com/support
Email: support@comsol.com
谢谢!

子康 赵
子康 赵
2023-07-02

工程师您好,当我想建立一个面积相对于其粗糙程度较大大的随机平面,此时显现出来图像的无限接近于一个平面,我该如何调整才可以使其粗糙程度很明显的显现出来呢?

Qihang Lin
Qihang Lin
2023-07-25 COMSOL 员工

可以尝试调整图形窗口中的光影效果进行图像展示。一般放大看应该能够看出存在粗糙的凹凸起伏,结果展示时也可截取部分图像结果进行展示。

Xun Zhang
Xun Zhang
2023-08-08

random_flat_surface_roughness.mph这个案例文件中存在Specify amplitude scale factor 1 和 Specify surface RMS height (Sq) 1这两个零件,如何使用他们创建指定幅值和RMS的表面呢

Anran Wei
Anran Wei
2023-08-11 COMSOL 员工

RMS高度指定了表面几何的统计特征,而幅值放大系数控制了表面形貌的直接值。两个参数之间没有直接的联系,所以很难写出直接形式的表达式以同时指定这两个参数。对于典型的作法,通常是选择指定RMS或者指定幅值放大系数,如random_flat_surface_roughness.mph案例中所示。

梦星 谢
梦星 谢
2023-11-23

你好,请问如何知道该随机表面的粗糙度呢

Anran Wei
Anran Wei
2023-11-28 COMSOL 员工

粗糙度可以由多种衡量准则来定义,比如轮廓算术平均偏差、轮廓最大高度等等。生成表面后,可以得到表面高度随空间的分布,然后结合您所使用的粗糙度衡量准则,在软件中编写相应的变量表达式,去计算对应的粗糙度值。

梦星 谢
梦星 谢
2023-11-23

Hi, how do I know the roughness of this random surface?

Anran Wei
Anran Wei
2023-11-29 COMSOL 员工

若您采用均方根RMS值作为粗糙度的衡量标准,在采用零件库中的“随机平面”创建几何时,会让您指定生成表面高度的RMS值

hh Sun
hh Sun
2023-12-03

您好,打开零件库中“随机平面”后,只有矩形宽度、矩形高度、比例因子、谱指数和空间频率分辨率这几个参数能改,请问RMS值在哪里改呢?

Anran Wei
Anran Wei
2023-12-06 COMSOL 员工

在 6.2 版本中,点击零件库中的“随机平面”,可以选择是“指定幅度比例因子”还是“指定表面RMS高度”的选项,后者可以在定义表面时输入RMS值

Zhu Hatcher
Zhu Hatcher
2023-12-12

工程师您好,请问这种方法可以用于生成具有随机粗糙度的曲面吗,比如球面或圆锥面等?

Qihang Lin
Qihang Lin
2023-12-14 COMSOL 员工

可以,但您需要获取表面的具体表达式代替使用的随机函数

Hy H
Hy H
2024-01-10

频谱指数和分形维数间有什么关系,如设定b,如何计算D呢?

浏览 COMSOL 博客