Orientation is the alignment of a particle in relation to Cartesian space. It is represented in EDEM as a 3-by-3 matrix that defines the rotation of the particle from its original orientation (as defined in the Creator and assigned an identity matrix).

由于浸入边界法(IBM)需要将固体颗粒的形状用函数来表示出来,这样 FLUENT 才能使用 UDF 来进行计算域与流体域的划分。而固体颗粒的堆积一般都是用离散单元法(DEM)来实现,譬如咱这边使用的 EDEM 的 颗粒工厂 功能。EDEM 中可以输出颗粒位置和方向的信息,但由于方向使用的是旋转矩阵(详见开头对方向的英语解释),在 UDF 中描述颗粒的堆积形状时不太方便(其实是不知道能不能用旋转矩阵来对此进行编写),因此这边找了方法来实现旋转矩阵与欧拉角的转换。

在之前,首先要说明的一点是,EDEM 中你可以设置颗粒初始的方向(旋转矩阵),默认为单位矩阵,而且有一个默认朝向,具体可以生成颗粒看看。然后经过颗粒的碰撞、堆积、稳定后,每个颗粒都会有一个方向(orientation),可以在 EDEM 的后处理中导出。

以下举例来说明。

创建颗粒工厂,选择 orientation 为 fixed,默认初始旋转矩阵为单位矩阵。

                                    | 1  0  0 |
                                  A=| 0  1  0 |
                                    | 0  0  0 |

堆积模拟结束后,导出颗粒的 orientation 信息,这边用一个颗粒举例。

                              | 0.765  -0.469  -0.441 |
                            M=| 0.469   0.875  -0.117 |
                              | 0.441  -0.117   0.89  |

打开 MATLAB ,由于咱使用的是 R2019a 版本,需要安装一个名为 Aerospace toolbox 的扩展,该扩展提供了旋转矩阵与欧拉角互转的函数,分别为 dcm2angleangle2dcm

M = [0.765, -0.469, -0.441; 0.469, 0.875, -0.117; 0.441, -0.117, 0.89];
[r1 r2 r3] = dcm2angle(M)

输出结果

r1 = -0.5500
r2 = 0.4567
r3 = -0.1307

使用 rad2deg 函数将弧度转为角度,结果如下:

r1 = -31.5113°
r2 = 26.1677°
r3 = -7.4892°

注意:r1,r2,r3 表示三次转动的转轴转角,默认转序为 ZYX,右手坐标系。
综上,咱就可以知道这个颗粒在堆积之后的方向为:先绕 Z 轴转动 -31°,再绕 Y 轴转动 26°,最后绕 X 轴转动 -7°。然后经过 亿点 处理就可以将该颗粒的函数表示出来啦。


Q.E.D

后记:后来发现不用将旋转矩阵转化成欧拉角,因为可以定义一个初始的颗粒的方向向量,然后用方向向量左乘旋转矩阵(因为原始方向是单位矩阵,所以旋转矩阵和终末颗粒的方向是一样的)得到新的颗粒的方向向量。然后用空间的距离公式什么的就可以完成颗粒形状函数的建立了。 :haha:

参考文献

[1] 欧拉角与旋转矩阵之间的转化公式及原理
[2] 旋转矩阵、四元数和欧拉角之间的转换——Matlab
[3] 三维坐标旋转矩阵算法


喜欢努力的雪集。