diff --git a/65.Derivation_of_the_governing_equations_for_the_MRF_approach.md b/65.Derivation_of_the_governing_equations_for_the_MRF_approach.md
new file mode 100644
index 0000000000000000000000000000000000000000..a7c5be29a722ff25c6b1f9873604c98ac7195cad
--- /dev/null
+++ b/65.Derivation_of_the_governing_equations_for_the_MRF_approach.md
@@ -0,0 +1,410 @@
+### 65 MRF方法的控制方程推导
+
+#### 65.1 初步观察
+
+使用多参考系（MRF）方法，必须将网格划分为不同的区域。由于OpenFOAM中的MRF方法仅支持旋转参考系[^228]，因此区域上只能施加旋转。被设为非零旋转的区域必须关于旋转轴对称。此外，OpenFOAM的旧版本[^229]仅支持稳定旋转，即角速度$\omega$恒定。
+
+##### 65.1.1 陷阱
+
+**MRF旋转的轴对称区域**
+
+如上所述，参考系旋转的MRF区域必须是轴对称的。然而，OpenFOAM不会对此执行任何检查。因此，可对一个关于$z$轴对称的区域指定一个绕$x$轴旋转的参考系。结果显然是错误的，求解器可能会崩溃。但是，如果够倒霉的话，求解器将不会崩溃，且很不幸地计算出错误的结果，浪费时间和计算资源。
+
+![](images/53.PNG) 始终检查`MRFProperties`中的`axis`是否与MRF区域对称轴重合。否则，应确保`MRFProperties`的`origin`位于MRF区域对称轴上。
+
+#### 65.2 质量守恒方程
+
+不可压缩流的质量守恒方程
+
+$$
+\nabla \cdot \mathbf{u} = 0
+\tag{400}
+\label{eq400}
+$$
+
+在所有惯性参考系中均有效。惯性参考系在空间和时间上是固定的，或者以恒定速度平移。 
+
+鉴于上一节中列出的限制，我们在MRF方法中仅考虑旋转参考系。使用旋转矩阵$\mathbf{Q}$将矢量的表达从惯性参考系转换至旋转参考系。旋转矩阵具有一个特性，其逆同样也是其转置。
+
+<p>
+$$
+{_R\mathbf{u}} = \mathbf{Q} \mathbf{u}
+\tag{401}
+\label{eq401}
+$$
+</p>
+
+<p>
+$$
+\mathbf{u} = {\mathbf{Q}^T}  {_R\mathbf{u}}
+\tag{402}
+\label{eq402}
+$$
+</p>
+
+<p>
+$$
+\mathbf{Q}^{-1} = \mathbf{Q}^T
+\tag{403}
+\label{eq403}
+$$
+</p>
+
+符号$\mathbf{u}$前的下标$R$表示，矢量$_R\mathbf{u}$是在旋转参照系中给出的。如果符号前没有下标，则矢量在惯性坐标系中给出。下标$R$放在符号前，以防止矢量$_R\mathbf{u}$被误认为相对速度。
+
+从惯性坐标系中的质量守恒方程出发，推导旋转坐标系中的质量守恒方程。
+
+<p>
+$$
+\nabla \cdot \mathbf{u} = 0
+\tag{404}
+\label{eq404}
+$$
+</p>
+
+<p>
+$$
+\nabla \cdot (\underbrace{{\mathbf{Q}^T} \mathbf{Q}}_{=\mathbf{I}} \mathbf{u}) = 0
+\tag{405}
+\label{eq405}
+$$
+</p>
+
+<p>
+$$
+\nabla \cdot ({\mathbf{Q}^T} {_R\mathbf{u}}) = 0
+\tag{406}
+\label{eq406}
+$$
+</p>
+
+运用关系式$ \nabla \cdot (\mathbf{A} \cdot \mathbf{a}) = (\nabla \cdot \mathbf{A}) \cdot \mathbf{a} + \mathbf{A} : (\nabla \mathbf{a}) $，注意旋转矩阵在空间中是恒定不变的。
+
+<p>
+$$
+\nabla \cdot ({\mathbf{Q}^T} {_R\mathbf{u}}) = 
+{\underbrace{(\nabla \cdot {\mathbf{Q}^T})}_{=0}} \cdot {_R \mathbf{u}} + 
+{\mathbf{Q}^T} : (\nabla {_R \mathbf{u}})
+\tag{407}
+\label{eq407}
+$$
+</p>
+
+<p>
+$$
+\nabla \cdot ({\mathbf{Q}^T} {_R\mathbf{u}}) = {\mathbf{Q}^T} : (\nabla {_R \mathbf{u}})
+\tag{408}
+\label{eq408}
+$$
+</p>
+
+<p>
+$$
+{\mathbf{Q}^T} : (\nabla {_R \mathbf{u}}) = 0
+\tag{409}
+\label{eq409}
+$$
+</p>
+
+接下来将方程左侧乘以旋转矩阵。
+
+$$
+\mathbf{Q} {\mathbf{Q}^T} : (\nabla {_R \mathbf{u}}) = 0
+\tag{410}
+\label{eq410}
+$$
+
+$$
+\mathbf{I} : (\nabla {_R \mathbf{u}}) = 0
+\tag{411}
+\label{eq411}
+$$
+
+单位张量和速度梯度的缩并等于速度的散度。
+
+$$
+\mathbf{I} : (\nabla {_R \mathbf{u}}) = \nabla \cdot {_R \mathbf{u}}
+\tag{412}
+\label{eq412}
+$$
+
+$$
+\nabla \cdot {_R \mathbf{u}} = 0
+\tag{413}
+\label{eq413}
+$$
+
+因此，以旋转坐标系表示速度的质量守恒方程与惯性坐标系下的质量守恒方程具有相同的公式。
+
+#### 65.3 动量守恒方程
+
+使用旋转坐标系时，可将流速分解为两个分量。第一个分量缘于参考系的旋转，第二个分量则是粒子或考虑的流体块与旋转参考系之间的相对运动。
+
+$$
+\mathbf{u} = \mathbf{\omega} \times \mathbf{r} + {\mathbf{u}_R}
+\tag{414}
+\label{eq414}
+$$
+
+使用旋转张量$\mathbf{\Omega}$，等式$\eqref{eq414}$亦可写成如下形式。 
+
+$$
+\mathbf{u} = \mathbf{\Omega} \mathbf{r} + {\mathbf{u}_R}
+\tag{415}
+\label{eq415}
+$$
+
+该旋转张量是一个斜对称张量，含有角速度矢量$\mathbf{\omega}$的各个分量。
+
+<p> 
+$$
+\mathbf{\Omega} = 
+\begin{pmatrix}
+ 0          & -\omega_z     &  \omega_y    \\
+ \omega_z   &  0            & -\omega_x    \\
+-\omega_y   &  \omega_x     &  0           \\
+\end{pmatrix}
+\tag{416}
+\label{eq416}
+$$
+</p>
+
+接下来，我们从不可压缩牛顿流体的动量守恒方程出发，推导旋转区域速度表示的动量方程。
+
+<p> 
+$$
+\frac{\partial \mathbf{u}}{\partial t} +
+(\nabla \mathbf{u}) \cdot \mathbf{u} = 
+-\frac{\nabla p}{\rho} +
+\nabla \cdot (\nu {\nabla \mathbf{u}})
+\tag{417}
+\label{eq417}
+$$
+</p>
+
+左侧项为$\mathbf{u}$的时间全导数。因此，可写为
+
+<p> 
+$$
+\frac{\mathrm{d} \mathbf{u}}{\mathrm{d} t} = 
+-\frac{\nabla p}{\rho} +
+\nabla \cdot (\nu {\nabla \mathbf{u}})
+\tag{418}
+\label{eq418}
+$$
+</p>
+
+利用等式$\eqref{eq414}$，可得
+
+<p> 
+$$
+\frac{\mathrm{d}}{\mathrm{d} t} (\mathbf{\Omega} \mathbf{r} + {\mathbf{u}_R}) = 
+-\frac{\nabla p}{\rho} +
+\nabla \cdot (\nu {\nabla \mathbf{u}})
+\tag{419}
+\label{eq419}
+$$
+</p>
+
+由于仅考虑稳定旋转，因此旋转张量的时间导数为零
+
+<p> 
+$$
+\frac{\mathrm{d} {\mathbf{u}_R}}{\mathrm{d} t} +
+\underbrace{\frac{\mathrm{d} \mathbf{\Omega}}{\mathrm{d} t}}_{=0} \mathbf{r} +
+\mathbf{\Omega} \frac{\mathrm{d} \mathbf{r}}{\mathrm{d} t} = 
+-\frac{\nabla p}{\rho} +
+\nabla \cdot (\nu {\nabla \mathbf{u}})
+\tag{420}
+\label{eq420}
+$$
+</p>
+
+<p> 
+$$
+\frac{\mathrm{d} {\mathbf{u}_R}}{\mathrm{d} t} +
+\mathbf{\Omega} \mathbf{u} = 
+-\frac{\nabla p}{\rho} +
+\nabla \cdot (\nu {\nabla \mathbf{u}})
+\tag{421}
+\label{eq421}
+$$
+</p>
+
+下面展开${\mathbf{u}_R}$的时间全导数
+
+<p> 
+$$
+\frac{\partial {\mathbf{u}_R}}{\partial t} +
+\underbrace{\frac{\partial {\mathbf{u}_R}}{\partial \mathbf{x}}}_{= \nabla {\mathbf{u}_R}} 
+\underbrace{\frac{\mathrm{d} \mathbf{x}}{\mathrm{d} t}}_{= \mathbf{u}} +
+\mathbf{\Omega} \mathbf{u} = 
+-\frac{\nabla p}{\rho} +
+\nabla \cdot (\nu {\nabla \mathbf{u}})
+\tag{422}
+\label{eq422}
+$$
+</p>
+
+<p> 
+$$
+\frac{\partial {\mathbf{u}_R}}{\partial t} +
+(\nabla {\mathbf{u}_R}) \cdot \mathbf{u} +
+\mathbf{\Omega} \mathbf{u} = 
+-\frac{\nabla p}{\rho} +
+\nabla \cdot (\nu {\nabla \mathbf{u}})
+\tag{423}
+\label{eq423}
+$$
+</p>
+
+将等式$\eqref{eq414}$中的${\mathbf{u}_R}$插入局部导数
+
+<p> 
+$$
+\frac{\partial}{\partial t} ( \mathbf{u} - \mathbf{\Omega} \mathbf{r}) +
+(\nabla {\mathbf{u}_R}) \cdot \mathbf{u} +
+\mathbf{\Omega} \mathbf{u} = 
+-\frac{\nabla p}{\rho} +
+\nabla \cdot (\nu {\nabla \mathbf{u}})
+\tag{424}
+\label{eq424}
+$$
+</p>
+
+因为由参考系稳定旋转导致的速度分量为常数，所以$\frac{\partial}{\partial t} (\mathbf{\Omega} \mathbf{r})$项为零
+
+<p> 
+$$
+\frac{\partial \mathbf{u}}{\partial t} +
+(\nabla {\mathbf{u}_R}) \cdot \mathbf{u} +
+\mathbf{\Omega} \mathbf{u} = 
+-\frac{\nabla p}{\rho} +
+\nabla \cdot (\nu {\nabla \mathbf{u}})
+\tag{425}
+\label{eq425}
+$$
+</p>
+
+使用如下等式，可重组左侧第二项
+
+<p> 
+$$
+\nabla \cdot (\mathbf{a} \mathbf{b}) = (\nabla \cdot \mathbf{b}) \mathbf{a} + (\nabla \mathbf{a}) \cdot \mathbf{b}
+\tag{426}
+\label{eq426}
+$$
+</p> 
+
+<p> 
+$$
+\frac{\partial \mathbf{u}}{\partial t} +
+\nabla \cdot ({\mathbf{u}_R} \mathbf{u}) +
+\mathbf{\Omega} \mathbf{u} = 
+-\frac{\nabla p}{\rho} +
+\nabla \cdot (\nu {\nabla \mathbf{u}})
+\tag{427}
+\label{eq427}
+$$
+</p>
+
+因此，我们使用相对于局部参考系的通量推导出了绝对速度的控制方程。
+
+在该控制方程中，区域旋转仅对两项产生影响。方程$\eqref{eq428}$左侧第二项包含相对速度。方程$\eqref{eq428}$右侧最后一项含有科里奥利力。
+
+<p> 
+$$
+\frac{\partial \mathbf{u}}{\partial t} +
+\nabla \cdot ({\mathbf{u}_R} \mathbf{u}) =
+-\frac{\nabla p}{\rho} +
+\nabla \cdot (\nu {\nabla \mathbf{u}}) -
+\mathbf{\Omega} \mathbf{u}
+\tag{428}
+\label{eq428}
+$$
+</p>
+
+方程$\eqref{eq428}$与Fluent理论手册[7]中以绝对速度表示的动量方程相符。[46]为OpenFOAM中MRF方法的另一参考资源。
+
+
+#### 65.4 MRF方法的实现说明
+
+MRF方法中，对运动区域的单元施加科里奥利力并不是唯一的必要操作。转子的边界条件必须要进行调整。转子运动时，转子壁面处的流体速度不为零。壁面处流速须等于转子的刚体旋转速度。
+
+##### 65.4.1 OpenFOAM-2.*
+
+OpenFOAM-2.2.x中，MRF方法是*fvOptions*机制的一部分[^230]。*fvOptions*是一个通用机制，考虑运行时可选的一些物理模型。*fvOptions*框架是动量方程中源项的广义化。通过该框架，可以选择MRF方法，以及动量（如风轮机风轮）、孔隙率（用于多孔区域模型）和能量源（例如用于传热区域）。
+
+为提供这种灵活性，在*fvOptions*框架的实现中，使用抽象类来定义广义源项的行为。派生类则实现实际的物理模型，例如MRF方法。因此，`MRFSource`类派生自`option`。
+
+`MRFSource`类的构造函数调用了`initialise()`方法。该方法在`MRFSource`类中定义，并调用了`MRFZone`类的`correctBoundaryVelocity`方法。在`correctBoundaryVelocity`方法中，MRF区域内的边界流速被设为刚体旋转速度。否则，无滑移边界条件会强制将绝对速度设为零，这显然是错误的。
+
+清单506中，可见对于MRF区域内所有面上刚体旋转速度的指定。在这些面上，指定了刚体旋转速度。
+
+<p> 
+$$
+\mathbf{u}_{rot} = \mathbf{\omega} \times (\mathbf{r}_{face} - \mathbf{r}_{origin})
+\tag{429}
+\label{eq429}
+$$
+</p>
+
+```cpp
+void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
+{
+    const vector Omega = this->Omega();
+    // Included patches
+    forAll(includedFaces_, patchi)
+    {
+        const vectorField& patchC = mesh_.Cf().boundaryField()[patchi];
+        vectorField pfld(U.boundaryField()[patchi]);
+        forAll(includedFaces_[patchi], i)
+        {
+            label patchFacei = includedFaces_[patchi][i];
+            pfld[patchFacei] = (Omega ^ (patchC[patchFacei] - origin_));
+        }
+        U.boundaryField()[patchi] == pfld;
+    }
+}
+```
+清单506：`MRFZone.C`文件中的`correctBoundaryVelocity`方法
+
+**MRF方法的功能和限制**
+
+OpenFOAM中的MRF方法仅涉及旋转，不同于FLUENT，还可考虑平移运动[7]。在这两个CFD软件中，运动参考系的速度均必须为常数。这意味着OpenFOAM能够处理以恒定角速度运动的旋转参考系。
+
+参考系旋转的区域，其边界朝向必须符合一定的条件，使得垂直于边界的参考系速度分量为零。这意味着，对于旋转参考系，其起作用的区域必须为圆柱，且圆柱的轴平行于参考系的旋转轴。
+
+FLUENT理论手册表示，MRF方法严格来说仅适用于稳态情况[7]。但是，FLUENT也为非稳态模拟提供了这种方法[7]。
+
+相比于MRF方法，滑移网格技术可提供更准确的结果，尤其对于瞬态模拟。但是，与动网格技术相比，MRF方法的主要优点是其对于计算成本的低需求。
+
+##### 65.4.2 OpenFOAM-3.*
+
+OpenFOAM-3.0.0[^232]之后，MRF方法从*fvOptions*框架中剥离。OpenFOAM开发者对此给出如下理由：
+
+>fvOptions不具有合适的结构来支持MRF，因为它通过用户指定的场量进行选项选择，然而MRF必须应用于特定求解器中的所有速度场。fvOptions中特定设计选择的一个后果，将会是难以为多相流支持MRF。并且，分别支持与参考系相关和与场量相关的选项，是更容易实现的。\
+当前提供的MRF功能仅支持旋转，但该结构将会被通用化，以支持其它参考系运动，包括线性加速度、SRF旋转和6DoF，这些运动将会是运行时可选的。
+
+如上引用所指出，MRF框架也进行了通用化。因此，旋转参考系将不再局限于稳定旋转。清单507显示了此次更改中启动参考系旋转的方式。这或许能改善模拟初始阶段的数值行为。 
+
+```
+MRF1
+{
+    cellZone    rotor;
+
+    active      yes;
+
+    nonRotatingPatches ();
+
+    origin    (0 0 0);
+    axis      (0 0 1);
+
+    omega   table
+            2(
+                (0  0.01)
+                (0.5 104.72)
+            );
+}
+```
+清单507：将表格作为角速度传递到MRF框架。
\ No newline at end of file