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] 三维坐标旋转矩阵算法