如何在 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] 上绘制表面:

显示正方形表面周期性的表面图。
正方形 [0,2] × [0,2] 上的表面周期性。表面高度用彩色表示。

β=0.5 的表面图。
β=1 的表面图。
β=1.8 的表面图。
β=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.

评论 (29)

正在加载...
涛铭 卢
涛铭 卢
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 员工

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

子康 赵
子康 赵
2024-03-12

十分感谢!

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 员工

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

子康 赵
子康 赵
2024-03-12

工程师您好,我利用案例生成了一个随机表面,我该如何操作,才能得到这个随机表面的函数表达式呢?

梦星 谢
梦星 谢
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呢?

Anran Wei
Anran Wei
2024-03-05 COMSOL 员工

本文介绍的是利用频谱的方式来生成粗糙表面,不是利用分形维数,两者很难说有怎样直接的转换联系,有关频谱生成的方式文中已经有介绍。

琦 陈
琦 陈
2024-03-08

工程师您好,通过文章已知产生三维随机表面可以通过指定amplitude或者指定RMS(Sq)两种方法,那么产生一条指定粗糙度(Rq)为某个值的随机曲线该怎么操作呢?

Hao Li
Hao Li
2024-03-11 COMSOL 员工

陈琦,您好!
感谢您的评论。

需要根据粗糙度定义反向分解为随机函数的组合,再依据本文提供的随机函数方式定义。
随机函数的组合需要您自行推导或查阅文献获取。

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

熠杰 陈
熠杰 陈
2024-04-23

您好,二维轴对称垂直于电脑屏幕的平面是不是没法设置粗糙度呀

 Min Yuan
Min Yuan
2024-04-28 COMSOL 员工

根据您的描述是不能设置,二维轴对称模型的几何具有连续旋转对称性,建议确认您的模型是否可以简化为二维轴对称。

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

琦 陈
琦 陈
2024-05-16

工程师您好!如果使用以上方法生成的随机表面,组成这个随机表面的各个点高度分布是否符合高斯分布或者正态分布?还是说表面各个点高度分布完全是随机的

Anran Wei
Anran Wei
2024-05-17 COMSOL 员工

换言之,每个点的高度都符合高斯分布,因为在每个点,它是许多高斯随机变量的和。此外,由于使用一系列三角函数进行构造,因此某一点的随机值与相邻点的随机值会有一定的相关性。

浏览 COMSOL 博客