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
的扩展,该扩展提供了旋转矩阵与欧拉角互转的函数,分别为 dcm2angle
和 angle2dcm
。
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] 三维坐标旋转矩阵算法
Comments | NOTHING