Merge remote-tracking branch 'origin/master'

copilot-conversations/20251103_101647__[[Tower]]分析主要内容.md

@@ -0,0 +1,984 @@
+---
+epoch: 1762136207636
+modelKey: gemini-2.5-flash|google
+tags:
+  - copilot-conversation
+---
+
+**user**: [[Tower]]分析主要内容
+[Timestamp: 2025/11/03 10:16:47]
+
+**ai**: 根据笔记 [[Tower]] 的内容，其主要分析了在动力学模型中，风力机塔筒的分布式特性所产生的各种广义力。以下是其核心内容的总结：
+
+### 1. 核心思想与模型简化
+笔记首先阐述了塔筒的分布式特性会带来与**塔筒弹性、阻尼、气动特性和重量**相关的广义惯性力和广义主动力。
+
+为了简化模型并使其更直观，笔记作者做出了以下关键决定：
+*   **移除质量调谐器 (Mass Tuners)**：在调整塔筒频率时，同时使用质量和刚度调谐器是多余的。保留刚度调谐器更为有利，因为质量调谐器会直接影响塔基载荷的计算，使相关方程变得非常复杂和不直观。
+*   **移除刚度调谐器对重力去刚化 (Gravitational Destiffening) 的影响**：同样是为了简化塔基载荷方程。这一做法是合理的，因为重力去刚化的影响相对于塔筒的整体刚度来说很小。
+
+### 2. 广义力的推导
+笔记详细推导了由不同物理效应产生的广义力，以及它们如何贡献于系统的运动方程 $[C(\mathbf{q},t)]\ddot{\mathbf{q}} = -f(\dot{\mathbf{q}},\mathbf{q},t)$。
+
+*   **广义惯性力 ($F_r^* \Big|_T$)**:
+    *   由塔筒自身的质量分布和塔顶的偏航轴承质量 (YawBrMass) 产生。
+    *   这部分贡献了系统的**质量矩阵 $[C(\mathbf{q},t)]_T$** 和与速度相关的**科里奥利/向心力项 $\{-f(\dot{q},q,t)\}|_T$**。
+
+*   **广义弹性力 ($F_{r \Big| ElasticT}$)**:
+    *   源于塔筒的弹性变形，通过应变能对广义坐标求偏导得出。
+    *   定义了前后 (Fore-Aft, TFA) 和侧向 (Side-to-Side, TSS) 的广义刚度矩阵 $k_{ij}^{'TFA}$ 和 $k_{ij}^{'TSS}$。这些刚度考虑了模态刚度调谐器 (`FAStTunr`, `SSStTunr`)，以便更好地匹配塔筒的固有频率。
+    *   这部分贡献了系统的**力向量 $\{-f(\dot{q},q,t)\}|_{\text{Elastic}T}$**。
+
+*   **广义阻尼力 ($F_r |_{\text{Damp}T}$)**:
+    *   采用瑞利阻尼 (Rayleigh damping) 模型，假设阻尼与刚度成正比。
+    *   阻尼力的大小取决于各阶模态的结构阻尼比 ($\zeta$)、广义刚度 ($k'$) 和不考虑塔顶质量影响的固有频率 ($f'$)。
+    *   这部分也贡献了系统的**力向量 $\{-f(\dot{q},q,t)\}|_{\text{Damp}T}$**。
+
+*   **广义重力 ($F_r |_{\text{Grav}T}$)**:
+    *   由塔筒沿高度分布的质量和塔顶质量在重力场中所受的力产生。
+    *   这部分贡献了系统的**力向量 $\{-f(\dot{q},q,t)\}|_{\text{Grav}T}$**。
+
+*   **广义气动力 ($F_r |_{\text{Aero}T}$)**:
+    *   由作用在塔筒上的分布式气动力和气动力矩产生。
+    *   这部分贡献了系统的**力向量 $\{-f(\dot{q},q,t)\}|_{\text{Aero}T}$**。
+
+总而言之，该笔记系统地分解了塔筒在动态响应中的受力情况，并为每一种力提供了详细的数学表达式，最终将它们整合到统一的动力学方程框架中。
+[Timestamp: 2025/11/03 10:17:07]
+
+**user**: [[Tower]]展开推导过程
+[Timestamp: 2025/11/03 10:22:19]
+
+**ai**: 好的，我们来详细展开笔记 [[Tower]] 中的推导过程。
+
+该笔记的推导遵循**凯恩方法 (Kane's Method)**，其核心思想是将所有力分为**广义主动力 (Generalized Active Forces)** $F_r$ 和**广义惯性力 (Generalized Inertia Forces)** $F_r^*$。系统的动力学方程为：
+$F_r + F_r^* = 0$
+
+在FAST v8的符号体系中，方程被写为 $[C(\mathbf{q},t)]\ddot{\mathbf{q}} = -f(\dot{\mathbf{q}},\mathbf{q},t)$。其中：
+*   $[C]\ddot{\mathbf{q}}$ 部分来自于广义惯性力中与广义加速度 $\ddot{q}$ 相关的项。
+*   $-f$ 部分则包含了广义惯性力中与广义速度 $\dot{q}$ 相关的项（如科里奥利力、离心力）以及所有的广义主动力（弹性力、阻尼力、重力、气动力等）。
+
+下面我们逐一分析每个部分的推导过程。
+
+---
+
+### 1. 广义惯性力 ($F_r^* \Big|_T$)
+
+这部分力是由于塔筒自身的质量和加速度产生的。
+
+*   **推导起点 (Eq. 1)**:
+    $$ F_r^* \Big|_T = - \int_0^{TwFlexL} \mu^T (h) \cdot {}^E v_r^T (h) \cdot {}^E a^T (h)dh - YawBrMass {}^E v_r^O \cdot {}^E a^O $$
+    这是凯恩方法中广义惯性力的标准定义。它由两部分组成：
+    1.  **分布式质量的贡献**：对整个柔性塔筒长度 (`TwrFlexL`) 进行积分。被积函数是塔筒上任意一点 `h` 处的达朗贝尔力 ($-\mu^T(h) \cdot {}^E a^T(h)$) 在第 `r` 个广义速度方向上的投影。
+        *   $\mu^T(h)$: 塔筒在高度 `h` 处的线密度（单位长度质量）。
+        *   ${}^E a^T(h)$: 塔筒在高度 `h` 处的绝对加速度。
+        *   ${}^E v_r^T(h)$: 塔筒在高度 `h` 处的点相对于第 `r` 个广义速度的**偏速度 (partial velocity)**。
+    2.  **集中质量的贡献**：位于塔顶的偏航轴承质量 (`YawBrMass`) 产生的惯性力。其形式与积分项类似，但作用于塔顶 `O` 点。
+
+*   **加速度项展开 (Eq. 2)**:
+    为了分离出与 $\ddot{q}$ 和 $\dot{q}$ 相关的项，需要将绝对加速度 ${}^E a$ 进行展开。任意一点的绝对速度 ${}^E \mathbf{v}$ 可以表示为广义速度的线性组合：${}^E \mathbf{v} = \sum_{i} {}^E \mathbf{v}_i \dot{q}_i + {}^E \mathbf{v}_t$。其对时间的导数即为加速度：
+    $$ {}^E \mathbf{a} = \frac{d}{dt}({}^E \mathbf{v}) = \sum_{i=1}^{10} {}^E \mathbf{v}_i \ddot{q}_i + \sum_{i=4}^{10} \frac{d}{dt}({}^E \mathbf{v}_i) \dot{q}_i $$
+    将这个展开式代入 (Eq. 1) 就得到了 (Eq. 2)。这个式子明确地将加速度分成了两部分：
+    1.  与 $\ddot{q}_i$ 成线性关系的部分。
+    2.  与 $\dot{q}_i$ 成线性关系的部分（包含偏速度的时间导数）。
+
+*   **分离质量矩阵和力向量 (Eq. 3 & 4)**:
+    将 (Eq. 2) 代入动力学方程 $F_r + F_r^* = 0$，并将所有项移到一边，整理成 $[C]\ddot{\mathbf{q}} + f = 0$ 的形式。
+    *   **质量矩阵 $[C]_T$ (Eq. 3)**：所有乘以 $\ddot{q}_{Col}$ 的系数项构成了质量矩阵的第 `(Row, Col)` 个元素。从 (Eq. 2) 中提取这些系数，就得到了 (Eq. 3)。
+    *   **力向量 $\{-f\}_T$ (Eq. 4)**：(Eq. 2) 中剩余的、与 $\dot{q}_i$ 相关的项，就构成了广义惯性力对力向量 $f$ 的贡献。注意 $f$ 前的负号，所以 (Eq. 4) 的表达式与 (Eq. 2) 中的第二部分符号相同。
+
+---
+
+### 2. 广义弹性力 ($F_{r \Big| ElasticT}$)
+
+这部分力来自于塔筒的弹性变形。
+
+*   **推导起点 (Eq. 5)**:
+    $$ F_{r \Big| ElasticT} = - \frac{\partial V'^T}{\partial q_r} $$
+    对于保守力（如弹性力），广义力可以由势能 $V'^T$ (这里指塔筒的弹性应变能) 对相应的广义坐标 $q_r$ 求偏导并取负号得到。
+
+*   **具体形式 (Eq. 6)**:
+    对于梁的弯曲变形，其应变能与位移的二阶导数（曲率）的平方成正比。经过推导，可以得到广义力是广义坐标的线性组合，即 $F = -[K]q$。 (Eq. 6) 就是这种线性关系的分量形式，其中 $k'_{ij}$ 是广义刚度矩阵的元素。它只对塔筒的自由度 (TFA1, TSS1, TFA2, TSS2) 有非零值。
+
+*   **广义刚度矩阵 (Eq. 7-10)**:
+    广义刚度 $k'_{ij}$ 的表达式来自于结构动力学中的标准推导：
+    $$ k'_{ij} = \int_0^{L} EI(h) \frac{d^2 \phi_i(h)}{dh^2} \frac{d^2 \phi_j(h)}{dh^2} dh $$
+    *   $EI(h)$: 沿塔高的抗弯刚度。
+    *   $\phi_i(h)$: 第 `i` 阶振动模态的振型函数。
+    *   $\frac{d^2 \phi(h)}{dh^2}$: 振型的曲率。
+    笔记中的 (Eq. 7) 和 (Eq. 9) 在此基础上增加了**刚度调谐器** `FAStTunr(i)` 和 `SSStTunr(i)`，允许用户微调每个模态的刚度，以更好地匹配实验或仿真得到的塔筒频率。
+
+*   **对系统方程的贡献 (Eq. 15 & 16)**:
+    *   弹性力只与位移 $q$ 有关，与 $\ddot{q}$ 无关，因此它对质量矩阵 $[C]$ 的贡献为零 (Eq. 15)。
+    *   弹性力是主动力的一部分，因此它贡献给力向量 $f$。将 (Eq. 6) 写入向量形式，就得到了 (Eq. 16)。
+
+---
+
+### 3. 广义阻尼力 ($F_r |_{\text{Damp}T}$)
+
+这里采用与刚度成正比的**瑞利阻尼**模型。
+
+*   **推导起点 (Eq. 11)**:
+    模态阻尼的标准形式是 $F_{damp} = -[C_{damp}]\dot{q}$，其中阻尼矩阵 $[C_{damp}]$ 的对角项为 $c_{ii} = 2\zeta_i \omega_i m_{ii}$。其中 $\omega_i = 2\pi f_i$ 是固有圆频率。
+    $c_{ii} = 2\zeta_i (2\pi f_i) m_{ii}$。
+    又因为 $k_{ii} = \omega_i^2 m_{ii} = (2\pi f_i)^2 m_{ii}$，所以 $m_{ii} = \frac{k_{ii}}{(2\pi f_i)^2}$。
+    代入可得：
+    $c_{ii} = 2\zeta_i (2\pi f_i) \frac{k_{ii}}{(2\pi f_i)^2} = \frac{\zeta_i k_{ii}}{\pi f_i}$。
+    因此，广义阻尼力为 $F_{r, damp} = -\sum_j c_{rj} \dot{q}_j = -\sum_j \frac{\zeta_j k'_{rj}}{\pi f'_j} \dot{q}_j$。这就是 (Eq. 11) 的来源。
+
+*   **频率和广义质量的定义 (Eq. 12-14)**:
+    为了计算阻尼系数，需要定义未耦合的固有频率 $f'$ (Eq. 12) 和广义质量 $m'$ (Eq. 13, 14)。这些都是模态分析中的标准定义，其中广义质量是通过对质量密度与振型函数乘积的积分得到的。
+
+*   **对系统方程的贡献 (Eq. 17 & 18)**:
+    *   阻尼力只与速度 $\dot{q}$ 有关，与 $\ddot{q}$ 无关，因此它对质量矩阵 $[C]$ 的贡献为零 (Eq. 17)。
+    *   阻尼力是主动力的一部分，贡献给力向量 $f$。将 (Eq. 11) 写入向量形式，得到 (Eq. 18)。
+
+---
+
+### 4. 广义重力 ($F_r |_{\text{Grav}T}$)
+
+*   **推导起点 (Eq. 19)**:
+    这是广义主动力的标准定义：$F_r = \int \mathbf{F}_{applied} \cdot \mathbf{v}_r dV$。
+    *   作用在塔筒微元上的重力为 $d\mathbf{F}_{grav} = -\mu^T(h) g \mathbf{z}_2 dh$。
+    *   将其与偏速度 ${}^E \mathbf{v}_r^T(h)$ 点乘，并沿塔高积分，得到分布式质量的贡献。
+    *   再加上塔顶集中质量的重力贡献 $(-\text{YawBrMass} \cdot g \mathbf{z}_2) \cdot {}^E \mathbf{v}_r^O$，就得到了 (Eq. 19)。
+
+*   **对系统方程的贡献 (Eq. 20 & 21)**:
+    *   重力与 $\ddot{q}$ 和 $\dot{q}$ 都无关，因此对 $[C]$ 的贡献为零 (Eq. 20)。
+    *   重力是主动力，贡献给力向量 $f$ (Eq. 21)。
+
+---
+
+### 5. 广义气动力 ($F_r |_{\text{Aero}T}$)
+
+*   **推导起点 (Eq. 22)**:
+    这同样是广义主动力的标准定义，但同时考虑了力和力矩：
+    $F_r = \int (\mathbf{F}_{applied} \cdot \mathbf{v}_r + \mathbf{M}_{applied} \cdot \boldsymbol{\omega}_r) dh$
+    *   $\mathbf{F}_{\text{Aero}T}^T(h)$: 作用在塔筒单位高度上的气动力。
+    *   $\mathbf{M}_{\text{Aero}T}^F(h)$: 作用在塔筒单位高度上的气动力矩。
+    *   ${}^E \boldsymbol{\omega}_r^F(h)$: 塔筒截面相对于第 `r` 个广义速度的**偏角速度 (partial angular velocity)**。
+    将它们点乘并沿塔高积分，得到总的广义气动力。
+
+*   **对系统方程的贡献 (Eq. 23 & 24)**:
+    *   气动力是外部施加的力，与 $\ddot{q}$ 无关，因此对 $[C]$ 的贡献为零 (Eq. 23)。
+    *   气动力是主动力，贡献给力向量 $f$ (Eq. 24)。
+
+通过以上步骤，笔记将塔筒的惯性、弹性、阻尼、重力和气动效应全部分解，并推导了它们各自对系统动力学方程（质量矩阵 $[C]$ 和力向量 $f$）的贡献。
+[Timestamp: 2025/11/03 10:23:05]
+
+**user**: 现在在此基础上做简化，仅考虑稳态过程，加速度为0
+[Timestamp: 2025/11/03 10:55:11]
+
+**ai**: 好的，我们基于笔记 [[Tower]] 的完整动力学方程，进行“仅考虑稳态过程，加速度为0”的简化。
+
+这个简化的核心物理意义是：我们不再分析塔筒的振动和动态响应，而是求解其在**恒定**外部载荷（如恒定的风）作用下的**静态平衡**位置。
+
+### 1. 简化假设
+
+1.  **稳态过程 (Steady-State)**: 系统的状态不随时间变化。
+2.  **加速度为零 ($\ddot{q} = 0$)**: 这是稳态的直接结果，意味着广义加速度为零。
+3.  **速度为零 ($\dot{q} = 0$)**: 在结构静力学中，稳态通常指**静态平衡 (static equilibrium)**，即系统不仅没有加速度，也没有速度。因此，所有广义速度也为零。
+
+### 2. 对动力学方程的影响
+
+原始的动力学方程是：
+$F_r^* + F_r = 0$
+或者写成矩阵形式：
+$[C(\mathbf{q},t)]\ddot{\mathbf{q}} = -f(\dot{\mathbf{q}},\mathbf{q},t)$
+
+根据我们的简化假设，我们来分析方程的各个组成部分会发生什么变化。
+
+#### a. 广义惯性力 ($F_r^* \Big|_T$)
+
+广义惯性力包含与加速度 ($\ddot{q}$) 和速度 ($\dot{q}$) 相关的项。
+
+*   **与加速度相关的项**:
+    由于我们假设 $\ddot{q} = 0$，整个质量矩阵项 $[C(\mathbf{q},t)]\ddot{\mathbf{q}}$ 变为零。这意味着在静态平衡中，质量和惯性不再起作用。
+
+*   **与速度相关的项 (Eq. 4)**:
+    这一项（代表科里奥利力和离心力）与 $\dot{q}_i$ 成正比。由于我们假设 $\dot{q} = 0$，这一项也变为零。
+    $$ \left\{ - f(\dot{q},q,t) \right\} \Big|_T = 0 $$
+
+**结论：在稳态（静态）假设下，整个广义惯性力 $F_r^* \Big|_T$ 为零。**
+
+#### b. 广义阻尼力 ($F_r |_{\text{Damp}T}$)
+
+阻尼力 (Eq. 11) 与广义速度 $\dot{q}$ 成正比。
+
+*   由于 $\dot{q} = 0$，所有阻尼力项都变为零。
+    $$ F_r |_{\text{Damp}T} = 0 $$
+
+**结论：在静态分析中，结构阻尼不产生任何力。**
+
+#### c. 保持不变的力
+
+以下这些力仅与系统的位置/位移 ($q$) 或外部条件有关，与运动状态（速度和加速度）无关，因此在简化后**依然存在**。
+
+1.  **广义弹性力 ($F_{r \Big| ElasticT}$)**:
+    该力由塔筒的变形产生，仅是广义坐标 $q$ 的函数 (Eq. 6)。它代表了结构抵抗变形的恢复力。
+    $$ F_{r \Big| ElasticT} = -[K']\mathbf{q} $$
+
+2.  **广义重力 ($F_r |_{\text{Grav}T}$)**:
+    该力由塔筒自身重量产生 (Eq. 19)。它的大小和方向取决于塔筒的几何形状和位置 $q$。
+
+3.  **广义气动力 ($F_r |_{\text{Aero}T}$)**:
+    在稳态分析中，我们假设风速是恒定的，因此作用在塔筒上的气动力和力矩也是恒定的 (Eq. 22)。这成为一个**恒定的外部载荷**。
+
+### 3. 简化后的最终方程
+
+原始的动力学方程 $F_r^* + F_r = 0$ 中，广义主动力 $F_r$ 是所有非惯性力的总和：
+$F_r = F_{r \Big| ElasticT} + F_r |_{\text{Damp}T} + F_r |_{\text{Grav}T} + F_r |_{\text{Aero}T}$
+
+将简化结果代入，我们得到：
+$0 + (F_{r \Big| ElasticT} + 0 + F_r |_{\text{Grav}T} + F_r |_{\text{Aero}T}) = 0$
+
+因此，描述塔筒稳态（静态）平衡的最终方程为：
+$$ F_{r \Big| ElasticT} + F_{r \Big| GravT} + F_{r \Big| AeroT} = 0 $$
+
+### 总结
+
+这个简化后的方程具有非常清晰的物理意义：
+
+**塔筒的内部弹性恢复力 ($F_{r \Big| ElasticT}$) 与作用在其上的外部载荷（重力 $F_{r \Big| GravT}$ 和稳态气动力 $F_{r \Big| AeroT}$）相平衡。**
+
+这本质上是一个**非线性静力学问题**。我们需要求解一组代数方程（而非微分方程）来找到唯一的广义坐标 $\mathbf{q}$，这个 $\mathbf{q}$ 就描述了塔筒在恒定风载和重力作用下的最终变形形状。
+[Timestamp: 2025/11/03 10:55:43]
+
+**user**: [[FillAugmat]][[Tower]] 解释其中与塔架动力学计算的内容
+[Timestamp: 2025/11/03 11:30:30]
+
+**ai**: 好的，我们来结合理论笔记 [[Tower]] 和代码笔记 [[FillAugmat]]，详细解释其中与塔架动力学计算相关的内容。
+
+[[Tower]] 笔记提供了风力机塔筒动力学模型的**理论基础和数学推导**。它详细描述了塔筒的广义惯性力、弹性力、阻尼力、重力和气动力是如何用数学公式表达的。
+
+[[FillAugmat]] 笔记中的 Fortran 子程序则是将这些理论公式**付诸实践的数值实现**。这个程序的核心任务是填充一个**增广矩阵 (Augmented Matrix)** `AugMat`。这个矩阵代表了整个风力机系统的动力学方程：
+
+$[C(\mathbf{q},t)]\ddot{\mathbf{q}} = -f(\dot{\mathbf{q}},\mathbf{q},t)$
+
+其中 `AugMat` 的结构可以理解为 `[ [C] | {-f} ]`。它的前 N 列是质量矩阵 `[C]`，最后一列 (`p%NAug`) 是力向量 `{-f}`。
+
+下面我们将 [[Tower]] 中的理论公式与 [[FillAugmat]] 中的代码实现一一对应起来进行解释。
+
+---
+
+### 1. 广义惯性力 (质量矩阵和惯性力)
+
+这部分对应于 $[C(\mathbf{q},t)]_T$ 和 $\{-f(\dot{q},q,t)\}|_T$ 的计算。
+
+#### **理论 ([[Tower]], Eq. 3 & 4):**
+
+*   **质量矩阵 $[C]_T$**: 由两部分组成：
+    1.  对塔筒分布式质量的积分：$\int_0^{TwrFlexL} \mu^T (h) {}^E \mathbf{v}_{Row}^T (h) \cdot {}^E \mathbf{v}_{Col}^T (h)dh$
+    2.  塔顶偏航轴承质量的贡献：$YawBrMass \cdot {}^E \mathbf{v}_{Row}^O \cdot {}^E \mathbf{v}_{Col}^O$
+*   **力向量 $\{-f\}_T$**: 同样由分布式质量和偏航轴承质量产生，与速度 $\dot{q}$ 相关（科里奥利力/离心力）。
+
+#### **代码实现 ([[FillAugmat]]):**
+
+代码通过将塔筒离散为多个节点 (`TwrNodes`)，用求和代替积分来实现。
+
+1.  **偏航轴承质量 (`YawBrMass`) 的贡献**:
+    ```fortran
+    ! Initialize the portions of the mass matrix ... using the yaw bearing mass effects.
+    DO L = 1,p%DOFs%NPTTE
+       DO I = L,p%DOFs%NPTTE
+          AugMat(p%DOFs%PTTE(I),p%DOFs%PTTE(L)) = p%YawBrMass*DOT_PRODUCT( RtHSdat%PLinVelEO(p%DOFs%PTTE(I),0,:), &     ! [C(q,t)]T of YawBrMass
+                                                       RtHSdat%PLinVelEO(p%DOFs%PTTE(L),0,:)    )
+       ENDDO
+    ENDDO
+    ```
+    *   **解释**: 这段代码完全对应理论公式的第二部分。它计算了偏航轴承质量 (`p%YawBrMass`) 对质量矩阵的贡献。`RtHSdat%PLinVelEO` 就是塔顶 O 点的偏速度 ${}^E \mathbf{v}^O$。`DOT_PRODUCT` 实现了向量点乘。
+
+2.  **塔筒分布式质量的贡献**:
+    ```fortran
+    ! Integrate to find the portions of the mass matrix ... associated with purely tower DOFs
+    DO J = 1,p%TwrNodes
+       ...
+       DO L = 1,p%DOFs%NPTTE
+          DO I = L,p%DOFs%NPTTE
+             AugMat(p%DOFs%PTTE(I),p%DOFs%PTTE(L)) = AugMat(...) + p%TElmntMass(J) * DOT_PRODUCT( RtHSdat%PLinVelET(J,...), RtHSdat%PLinVelET(J,...) )
+          ENDDO
+       ENDDO
+       ...
+    ENDDO ! J - Tower nodes / elements
+    ```
+    *   **解释**: 这段代码通过 `DO J = 1,p%TwrNodes` 循环，对所有塔筒单元进行遍历，实现了数值积分。`p%TElmntMass(J)` 相当于理论中的 $\mu^T(h)dh$，即单元质量。`RtHSdat%PLinVelET(J,...)` 则是塔筒 J 节点的偏速度 ${}^E \mathbf{v}^T(h)$。这完全对应了理论公式的第一部分。
+
+3.  **惯性相关的力向量贡献**:
+    ```fortran
+    ! Forcing vector from YawBrMass
+    TmpVec1 = -p%YawBrMass*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccEOt )
+    DO I = 1,p%DOFs%NPTTE
+          AugMat(p%DOFs%PTTE(I), p%NAug) = DOT_PRODUCT( RtHSdat%PLinVelEO(...), TmpVec1 )
+    ENDDO
+    
+    ! Forcing vector from distributed mass (integrated)
+    DO J = 1,p%TwrNodes
+        TmpVec1 = ... - p%TElmntMass(J)*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccETt(:,J) )
+        ...
+        DO I = 1,p%DOFs%NPTTE
+            AugMat(p%DOFs%PTTE(I), p%NAug) = AugMat(...) + DOT_PRODUCT( RtHSdat%PLinVelET(J,...), TmpVec1 )
+        ENDDO
+    ENDDO
+    ```
+    *   **解释**: 这部分计算了力向量 `{-f}`。代码中的 `LinAccEOt` 和 `LinAccETt` 是与 $\ddot{q}$ 无关的加速度项（即科里奥利/离心加速度）。`TmpVec1` 计算了达朗贝尔力（`-m*a`）和重力（`-m*g`）的总和。然后通过与偏速度 `PLinVelE...` 点乘，计算出它们在力向量 `{-f}` 中的贡献。这与 [[Tower]] 中的 Eq. 4 (惯性力) 和 Eq. 21 (重力) 的实现相对应。
+
+---
+
+### 2. 广义弹性和阻尼力
+
+这部分对应于 $\{-f(\dot{q},q,t)\}|_{\text{Elastic}T}$ 和 $\{-f(\dot{q},q,t)\}|_{\text{Damp}T}$ 的计算。
+
+#### **理论 ([[Tower]], Eq. 6, 11, 16, 18):**
+
+*   **弹性力**: 形式为 $-k_{ij}q_j$，是广义位移 $q$ 的线性函数。
+*   **阻尼力**: 形式为 $-c_{ij}\dot{q}_j$，是广义速度 $\dot{q}$ 的线性函数。
+*   它们都只贡献给力向量 `{-f}`，对质量矩阵 `[C]` 没有贡献。
+
+#### **代码实现 ([[FillAugmat]]):**
+
+```fortran
+! If the associated DOFs are enabled, add the tower elasticity and damping forces to the forcing vector
+IF ( p%DOF_Flag(DOF_TFA1) )  THEN
+   AugMat(DOF_TFA1,p%NAug) = AugMat(DOF_TFA1,p%NAug) &
+                           - p%KTFA(1,1)*x%QT( DOF_TFA1) - p%KTFA(1,2)*x%QT( DOF_TFA2) &  ! Elasticity
+                           - p%CTFA(1,1)*x%QDT(DOF_TFA1) - p%CTFA(1,2)*x%QDT(DOF_TFA2)    ! Damping
+ENDIF
+! ... (similar blocks for TSS1, TFA2, TSS2) ...```
+*   **解释**: 这段代码是理论公式的**直接翻译**。
+    *   `p%KTFA` 和 `p%KTSS` 是广义刚度矩阵 $k'_{ij}$。
+    *   `p%CTFA` 和 `p%CTSS` 是广义阻尼矩阵 $c_{ij}$ (即理论中的 $\frac{\zeta k'}{\pi f'}$ 项)。
+    *   `x%QT(...)` 是广义位移 $q$。
+    *   `x%QDT(...)` 是广义速度 $\dot{q}$。
+    *   计算结果被直接累加到增广矩阵的最后一列 (`p%NAug`)，即力向量 `{-f}` 中，这与理论完全吻合。
+
+---
+
+### 3. 耦合项：机舱和叶轮对塔架的影响
+
+[[Tower]] 笔记主要关注由塔筒自身属性产生的力。但在一个完整的系统中，塔筒的运动还受到其上部结构（机舱、发电机、轮毂、叶片等）的惯性影响。
+
+#### **理论**:
+
+当塔筒运动时（例如，一阶前后弯曲 `TFA1`），它会带动整个上部结构（Rotor-Nacelle Assembly, RNA）一起运动。根据达朗贝尔原理，RNA 的质量和加速度会产生一个惯性力，这个力反作用于塔筒，从而影响塔筒的动力学方程。这体现为质量矩阵 `[C]` 和力向量 `{-f}` 中的**耦合项**。
+
+#### **代码实现 ([[FillAugmat]]):**
+
+在程序的后半部分，通过所谓的**“局部载荷 (partial loads)”**来计算这些耦合项。
+
+```fortran
+! Now that all of the partial loads have been found, let's fill in the portions of the mass matrix ...
+IF ( p%DOF_Flag (DOF_TFA1) )  THEN
+   DO I = p%DOFs%Diag(DOF_TFA1),p%DOFs%NActvDOF
+      AugMat(p%DOFs%SrtPS(I),DOF_TFA1) = AugMat(...) &
+                                       - DOT_PRODUCT( RtHSdat%PLinVelEO(DOF_TFA1,0,:), RtHSdat%PFrcONcRt(:,p%DOFs%SrtPS(I)) ) &
+                                       - DOT_PRODUCT( RtHSdat%PAngVelEB(DOF_TFA1,0,:), RtHSdat%PMomBNcRt(:,p%DOFs%SrtPS(I)) )
+   ENDDO
+   AugMat(DOF_TFA1, p%NAug) = AugMat(...) &
+                           + DOT_PRODUCT( RtHSdat%PLinVelEO(DOF_TFA1,0,:), RtHSdat%FrcONcRtt  ) &
+                           + DOT_PRODUCT( RtHSdat%PAngVelEB(DOF_TFA1,0,:), RtHSdat%MomBNcRtt  )
+ENDIF
+```
+*   **解释**:
+    *   `RtHSdat%PFrcONcRt` 和 `RtHSdat%PMomBNcRt` 是关键变量，它们代表了由 RNA 运动产生的**局部力和力矩**。
+    *   `RtHSdat%PLinVelEO(DOF_TFA1,0,:)` 是塔筒 `TFA1` 自由度在塔顶 O 点产生的偏速度。
+    *   代码通过将塔顶的偏速度/偏角速度与 RNA 的局部载荷做点积，计算出 RNA 的惯性对塔筒自由度 (`DOF_TFA1`) 对应方程的贡献。这部分计算了质量矩阵中**塔筒自由度**与其他**非塔筒自由度**之间的耦合项，以及 RNA 惯性力对塔筒力向量的贡献。
+
+### 总结
+
+*   [[Tower]] 提供了“做什么”的理论蓝图，定义了所有需要计算的物理量。
+*   [[FillAugmat]] 提供了“怎么做”的数值方法，它通过以下步骤实现了理论：
+    1.  **离散化**: 将连续的塔筒模型变为一系列离散的单元和节点。
+    2.  **数值积分**: 用对所有单元的求和来代替理论公式中的积分。
+    3.  **直接计算**: 对于像弹簧和阻尼这样的集中力，直接根据代数公式计算。
+    4.  **模块化/局部载荷**: 使用预先计算好的上部结构载荷来高效地处理复杂的耦合动力学问题。
+
+通过这种方式，代码精确地将 [[Tower]] 笔记中描述的连续体力学模型转化为了计算机可以求解的离散线性代数问题。
+[Timestamp: 2025/11/03 11:31:18]
+
+**user**: 现在在此基础上做简化，仅考虑稳态过程，加速度为0，求解q，代码该如何整理
+[Timestamp: 2025/11/03 14:28:44]
+
+**ai**: 好的，我们来探讨如何基于 [[Tower]] 的理论和 [[FillAugmat]] 的代码实现，进行简化以求解稳态（静态）下的广义位移 `q`。
+
+### 1. 理论基础：从动态到静态
+
+如前所述，稳态分析的核心是移除所有与时间变化相关的项。
+原始动态方程:
+$[C(\mathbf{q},t)]\ddot{\mathbf{q}} + f_{damp}(\dot{\mathbf{q}}) + f_{inertial}(\dot{\mathbf{q}}, \mathbf{q}) + f_{elastic}(\mathbf{q}) = F_{ext}(\mathbf{q})$
+
+简化假设:
+*   加速度 $\ddot{\mathbf{q}} = 0$
+*   速度 $\dot{\mathbf{q}} = 0$
+
+简化后的静态平衡方程:
+$f_{elastic}(\mathbf{q}) = F_{ext}(\mathbf{q})$
+
+其中：
+*   $f_{elastic}(\mathbf{q})$ 是由结构弹性产生的恢复力，可以写成 $[K]\mathbf{q}$ 的形式，其中 $[K]$ 是**刚度矩阵**。
+*   $F_{ext}(\mathbf{q})$ 是所有外部载荷（重力、稳态气动力等）产生的广义力，组成**外力向量**。
+
+因此，我们的最终目标是求解一个线性代数方程组：
+$$ [K]\mathbf{q} = \mathbf{F}_{ext} $$
+解出这个方程，就能得到稳态位移向量 $\mathbf{q}$。
+
+### 2. 代码整理思路：从 `FillAugmat` 到 `SolveStatic`
+
+`FillAugmat` 的任务是构建增广矩阵 `[ [C] | {-f} ]`。我们的新任务是构建刚度矩阵 `[K]` 和外力向量 `F_ext`。
+
+这意味着我们需要修改 `FillAugmat` 的逻辑，将其转变为一个新的子程序，例如 `SolveStaticTower`。这个新程序将执行以下操作：
+
+1.  **忽略所有质量/惯性项**：所有填充 `AugMat` 前 N 列（质量矩阵 `[C]`）的代码都可以被移除。所有与 `LinAcc...t` (科里奥利/离心加速度) 相关的计算也可以移除。
+2.  **忽略所有阻尼项**：所有与 `x%QDT` (广义速度) 相乘的项（如 `p%CTFA * x%QDT`）都可以被移除。
+3.  **分离刚度项和外力项**：
+    *   **刚度矩阵 `[K]`**: 来自于弹性力。在代码中，形如 `- p%KTFA(1,1)*x%QT(DOF_TFA1)` 的项，其系数 `p%KTFA(1,1)` 就是刚度矩阵的一个元素。
+    *   **外力向量 `F_ext`**: 来自于重力、气动力等。在代码中，这些力被计算并累加到 `AugMat` 的最后一列。我们需要把这些项收集起来，形成 `F_ext`。
+
+### 3. 伪代码实现 (`SolveStaticTower`)
+
+下面是一个基于 `FillAugmat` 逻辑整理出的伪代码，用于求解塔筒的静态位移 `q`。
+
+```fortran
+SUBROUTINE SolveStaticTower( p, u, q_solution )
+! p: Parameters, u: Inputs (aero forces), q_solution: Output static displacements
+
+    ! 1. 初始化
+    INTEGER, PARAMETER :: N_DOF = ... ! 系统自由度总数
+    REAL(R8Ki) :: GlobalStiffnessMatrix(N_DOF, N_DOF)
+    REAL(R8Ki) :: ExternalForceVector(N_DOF)
+    REAL(R8Ki) :: q_solution(N_DOF)
+    
+    GlobalStiffnessMatrix = 0.0
+    ExternalForceVector   = 0.0
+
+    ! ==========================================================================
+    ! 2. 填充外力向量 F_ext (来自重力和气动力)
+    !    这部分逻辑从 FillAugmat 中提取，但只保留与外力相关的计算
+    ! ==========================================================================
+
+    ! a. 塔顶偏航轴承的重力贡献
+    !    (原FillAugmat中计算TmpVec1的部分)
+    TmpVec_Grav_YawBr = -p%YawBrMass * p%Gravity * CoordSys%z2
+    DO I = 1, p%DOFs%NPTTE
+        ExternalForceVector(p%DOFs%PTTE(I)) = ExternalForceVector(p%DOFs%PTTE(I)) + &
+            DOT_PRODUCT( RtHSdat%PLinVelEO(p%DOFs%PTTE(I),0,:), TmpVec_Grav_YawBr )
+    ENDDO
+
+    ! b. 塔筒分布式质量的重力贡献 (通过积分/求和)
+    DO J = 1, p%TwrNodes
+        ! 假设有稳态气动力 FSAero_static
+        TmpVec_Forces_Elmnt = ( FSAero_static(:,J) )*abs(p%DHNodes(J)) &  ! 稳态气动力
+                              - p%TElmntMass(J) * p%Gravity * CoordSys%z2      ! 重力
+        ! 假设有稳态气动力矩 MMAero_static
+        TmpVec_Moments_Elmnt = ( MMAero_static(:,J) )*abs(p%DHNodes(J))
+
+        DO I = 1, p%DOFs%NPTTE
+            ExternalForceVector(p%DOFs%PTTE(I)) = ExternalForceVector(p%DOFs%PTTE(I)) + &
+                DOT_PRODUCT( RtHSdat%PLinVelET(J,p%DOFs%PTTE(I),0,:), TmpVec_Forces_Elmnt ) + &
+                DOT_PRODUCT( RtHSdat%PAngVelEF(J,p%DOFs%PTTE(I),0,:), TmpVec_Moments_Elmnt )
+        ENDDO
+    ENDDO
+
+    ! c. 来自机舱和叶轮的重力和稳态气动载荷 (原FillAugmat中与partial loads相关的部分)
+    !    这部分会更复杂，需要计算RNA的重力和稳态气动力对塔筒自由度的贡献
+    !    ... (此处省略，但逻辑是相似的：力 * 偏速度) ...
+
+
+    ! ==========================================================================
+    ! 3. 填充刚度矩阵 [K]
+    !    这部分逻辑从 FillAugmat 中提取，但只保留与弹性力相关的系数
+    ! ==========================================================================
+
+    ! a. 塔筒自身的刚度
+    !    原代码是 AugMat(DOF, NAug) = AugMat(...) - K*q
+    !    我们现在把 K 填入 GlobalStiffnessMatrix
+    IF ( p%DOF_Flag(DOF_TFA1) ) THEN
+        GlobalStiffnessMatrix(DOF_TFA1, DOF_TFA1) = GlobalStiffnessMatrix(DOF_TFA1, DOF_TFA1) + p%KTFA(1,1)
+        GlobalStiffnessMatrix(DOF_TFA1, DOF_TFA2) = GlobalStiffnessMatrix(DOF_TFA1, DOF_TFA2) + p%KTFA(1,2)
+    ENDIF
+    IF ( p%DOF_Flag(DOF_TSS1) ) THEN
+        GlobalStiffnessMatrix(DOF_TSS1, DOF_TSS1) = GlobalStiffnessMatrix(DOF_TSS1, DOF_TSS1) + p%KTSS(1,1)
+        GlobalStiffnessMatrix(DOF_TSS1, DOF_TSS2) = GlobalStiffnessMatrix(DOF_TSS1, DOF_TSS2) + p%KTSS(1,2)
+    ENDIF
+    IF ( p%DOF_Flag(DOF_TFA2) ) THEN
+        GlobalStiffnessMatrix(DOF_TFA2, DOF_TFA1) = GlobalStiffnessMatrix(DOF_TFA2, DOF_TFA1) + p%KTFA(2,1)
+        GlobalStiffnessMatrix(DOF_TFA2, DOF_TFA2) = GlobalStiffnessMatrix(DOF_TFA2, DOF_TFA2) + p%KTFA(2,2)
+    ENDIF
+    IF ( p%DOF_Flag(DOF_TSS2) ) THEN
+        GlobalStiffnessMatrix(DOF_TSS2, DOF_TSS1) = GlobalStiffnessMatrix(DOF_TSS2, DOF_TSS1) + p%KTSS(2,1)
+        GlobalStiffnessMatrix(DOF_TSS2, DOF_TSS2) = GlobalStiffnessMatrix(DOF_TSS2, DOF_TSS2) + p%KTSS(2,2)
+    ENDIF
+    
+    ! 注意：实际的刚度矩阵可能还包含重力引起的几何刚度效应，这里为了简化暂未列出。
+
+    ! ==========================================================================
+    ! 4. 求解线性方程组 [K]q = F_ext
+    ! ==========================================================================
+    
+    ! 调用一个线性代数库 (如 LAPACK) 的求解器来求解 q
+    ! 例如，使用 DGESV 求解 A*X = B, 这里 A=GlobalStiffnessMatrix, B=ExternalForceVector, X=q_solution
+    CALL LAPACK_DGESV( N_DOF, 1, GlobalStiffnessMatrix, N_DOF, ipiv, ExternalForceVector, N_DOF, info )
+    
+    ! 求解后，ExternalForceVector 会被解 q_solution 覆盖
+    q_solution = ExternalForceVector
+
+END SUBROUTINE SolveStaticTower
+```
+
+### 总结与关键点
+
+1.  **核心转变**：从构建 `[ [C] | {-f} ]` 并求解 $\ddot{q}$，转变为构建 `[K]` 和 `F_ext` 并求解 `q`。
+2.  **代码复用**：`FillAugmat` 中的大部分计算逻辑（如偏速度、点积）是可以复用的，但需要重新组织。
+3.  **外力向量 `F_ext`**：主要来源于**重力**和**稳态气动力**。在原代码中，它们是 `{-f}` 的一部分。在整理时，需要注意符号，因为 $F_{ext} = -f_{ext}$。
+4.  **刚度矩阵 `[K]`**：主要来源于**结构弹性**。在原代码中，它们是 `{-f}` 中与 `x%QT` 相乘的系数。
+5.  **求解**：最后一步是调用标准的线性方程组求解器，这是与原动态时域积分完全不同的数学步骤。
+
+通过这样的整理，就可以将一个复杂的动态仿真代码简化为一个静态结构分析程序，用于计算风力机在特定工况下的静态变形。
+[Timestamp: 2025/11/03 14:29:25]
+
+**user**: [[FillAugmat]]中\left[ \sum_{i=4}^{10} \frac{d}{dt} \left( {}^E \mathbf{v}_i^O \right) \dot{q}_i \right] 在哪里体现
+[Timestamp: 2025/11/03 14:44:54]
+
+**ai**: Excellent question. This gets to the heart of how the theoretical dynamics equations are translated into a computational algorithm.
+
+The term you've highlighted, $\left[ \sum_{i=4}^{10} \frac{d}{dt} \left( {}^E \mathbf{v}_i^O \right) \dot{q}_i \right]$, represents the part of the absolute acceleration of the tower top (point O) that is **not** proportional to the generalized accelerations ($\ddot{q}$). This term physically corresponds to the **Coriolis and centrifugal accelerations** of point O.
+
+In the `FillAugmat` code, this complex mathematical term is not calculated directly within the subroutine itself. Instead, it is pre-calculated in a preceding kinematics routine (part of the `RtHndSid` module) and passed into `FillAugmat` through the `RtHSdat` data structure.
+
+The specific variable that embodies this term is:
+
+**`RtHSdat%LinAccEOt`**
+
+Let's break down how this is implemented in the code.
+
+### 1. Theoretical Context ([[Tower]], Eq. 4)
+
+The contribution of the yaw bearing mass (`YawBrMass`) to the force vector `{-f}` is given by:
+$$ \{-f...\} = ... - YawBrMass \cdot {}^E \mathbf{v}_{Row}^O \cdot \left[ \sum_{i=4}^{10} \frac{d}{dt} \left( {}^E \mathbf{v}_i^O \right) \dot{q}_i \right] $$
+This formula states that we need to compute the Coriolis/centrifugal acceleration of point O, multiply it by the negative of the yaw bearing mass to get an inertial force, and then project that force onto the partial velocity vector ${}^E \mathbf{v}_{Row}^O$ to get its contribution to the generalized force.
+
+### 2. Code Implementation ([[FillAugmat]])
+
+The code implements this in two steps:
+
+**Step 1: Calculate the total inertial force associated with the yaw bearing.**
+
+This happens in the following line:
+```fortran
+TmpVec1 = -p%YawBrMass*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccEOt ) ! The portion of FrcT0Trbt associated with the YawBrMass
+```
+Here:
+*   `p%YawBrMass` is the `YawBrMass` from the theory.
+*   `p%Gravity*CoordSys%z2` is the gravitational acceleration vector $\mathbf{g}$.
+*   **`RtHSdat%LinAccEOt`** is the crucial variable. Its name can be decoded as:
+    *   `LinAcc`: **Lin**ear **Acc**eleration
+    *   `E`: In the **E**arth (inertial) frame
+    *   `O`: Of point **O** (the yaw bearing/tower top)
+    *   `t`: The part of the acceleration that is **t**ime-dependent but **independent of $\ddot{q}$**. This is precisely the Coriolis and centrifugal acceleration term.
+
+So, `TmpVec1` calculates the total force ($-m\mathbf{a} - m\mathbf{g}$) acting on the yaw bearing, where $\mathbf{a}$ is the Coriolis/centrifugal part of the acceleration.
+
+**Step 2: Project this force to get the generalized force.**
+
+This is done in the subsequent loop:
+```fortran
+DO I = 1,p%DOFs%NPTTE    ! Loop through all active (enabled) tower DOFs...
+      AugMat(p%DOFs%PTTE(I),   p%NAug) = DOT_PRODUCT( RtHSdat%PLinVelEO(p%DOFs%PTTE(I),0,:), &     ! {-f(qd,q,t)}T + {-f(qd,q,t)}GravT of YawBrMass
+                                                   TmpVec1                   )
+ENDDO
+```
+Here:
+*   `AugMat(..., p%NAug)` is the force vector `{-f}`.
+*   `RtHSdat%PLinVelEO(...)` is the partial velocity of point O, which corresponds to ${}^E \mathbf{v}_{Row}^O$ in the theory.
+*   `DOT_PRODUCT(..., TmpVec1)` performs the vector dot product.
+
+When you expand the `DOT_PRODUCT` operation, you get:
+`DOT_PRODUCT(v_r^O, -YawBrMass*(g + LinAccEOt))`
+$= -YawBrMass \cdot \text{DOT_PRODUCT}(v_r^O, g) - YawBrMass \cdot \text{DOT_PRODUCT}(v_r^O, \text{LinAccEOt})$
+
+This perfectly matches the theoretical contributions to the force vector from:
+1.  The generalized gravity force (from [[Tower]], Eq. 21).
+2.  The generalized Coriolis/centrifugal force (from [[Tower]], Eq. 4).
+
+### Summary
+
+The term $\left[ \sum_{i=4}^{10} \frac{d}{dt} \left( {}^E \mathbf{v}_i^O \right) \dot{q}_i \right]$ is **represented by the variable `RtHSdat%LinAccEOt`** in the `FillAugmat` subroutine. The code uses a modular approach where complex kinematic quantities like this are calculated beforehand and then used within `FillAugmat` to construct the final terms of the system's equations of motion.
+[Timestamp: 2025/11/03 14:45:22]
+
+**user**: [[FillAugmat]]基于fillaugmat函数简化出一个fill_steady_augmat函数，变量名保持一致，求解塔架稳态q，后面用LU分解求解q
+[Timestamp: 2025/11/03 15:17:00]
+
+**ai**: 好的，我们来基于 `FillAugMat` 函数，简化出一个 `fill_steady_augmat` 函数，用于在稳态（加速度为0，速度为0）条件下构建刚度矩阵 `[K]` 和外部载荷向量 `F_ext`，以便后续使用 LU 分解求解广义位移 `q`。
+
+### 1. 简化原则
+
+1.  **稳态条件**: 广义加速度 $\ddot{\mathbf{q}} = 0$ 和广义速度 $\dot{\mathbf{q}} = 0$。
+2.  **方程形式**: 最终目标是构建并返回 `[K]` 和 `F_ext`，以求解线性方程组 `[K]q = F_ext`。
+3.  **变量名一致**: 尽可能保持原始 `FillAugMat` 中的变量名。
+4.  **忽略动态项**:
+    *   所有与质量矩阵 `[C]` 相关的计算（即 `AugMat` 的前 N 列）将被移除，因为在静态分析中，惯性力为零。
+    *   所有与广义加速度 `RtHSdat%LinAccESt`, `RtHSdat%LinAccEOt`, `RtHSdat%AngAccEGt` 相关的项将被视为零。
+    *   所有与广义速度 `x%QDT` 相关的项（阻尼力、科里奥利力、离心力）将被移除，因为在静态分析中，速度为零。
+5.  **分离刚度项和外力项**:
+    *   原始 `FillAugMat` 中，弹性力项形如 `-K*q`，并被加到 `AugMat` 的最后一列 (`p%NAug`)。在新函数中，`K` 的系数将被提取出来填充 `GlobalStiffnessMatrix`。
+    *   原始 `FillAugMat` 中，重力、气动力、水动力、刹车力等外部载荷被加到 `AugMat` 的最后一列 (`p%NAug`)。在新函数中，这些项将被累加到 `ExternalForceVector`。
+
+### 2. `fill_steady_augmat` 子程序伪代码
+
+```fortran
+SUBROUTINE fill_steady_augmat( p, x_static, CoordSys, u_static, HSSBrTrq_static, RtHSdat_static, &
+                               GlobalStiffnessMatrix, ExternalForceVector )
+!..................................................................................................................................
+! This subroutine constructs the global stiffness matrix [K] and the external force vector F_ext
+! for a steady-state (static) analysis, where accelerations and velocities are assumed to be zero.
+! The output can then be used to solve [K]q = F_ext for static displacements q.
+!..................................................................................................................................
+
+      ! Passed variables
+   TYPE(ED_ParameterType),       INTENT(IN   )  :: p           !< Parameters
+   TYPE(ED_ContinuousStateType), INTENT(IN   )  :: x_static    !< Continuous states (only x_static%QT for displacements is relevant)
+   TYPE(ED_CoordSys),            INTENT(IN   )  :: CoordSys    !< The coordinate systems
+   TYPE(ED_InputType),           INTENT(IN   )  :: u_static    !< The static aero blade forces/moments
+   TYPE(ED_RtHndSide),           INTENT(IN   )  :: RtHSdat_static ! Data from RtHndSid module (assumed to contain static forces/moments and zero accelerations)
+   REAL(ReKi),                   INTENT(IN )    :: HSSBrTrq_static ! Static HSS brake torque
+   REAL(R8Ki),                   INTENT(OUT)    :: GlobalStiffnessMatrix(:,:) !< The global stiffness matrix [K]
+   REAL(R8Ki),                   INTENT(OUT)    :: ExternalForceVector(:)     !< The external force vector F_ext
+   
+      ! Local variables (keeping original names for consistency, but some will be effectively zero)
+   REAL(ReKi)                   :: TmpVec    (3)                                   ! A temporary vector used in various computations.
+   REAL(ReKi)                   :: TmpVec1   (3)                                   ! A temporary vector used in various computations.
+   REAL(ReKi)                   :: TmpVec3   (3)                                   ! A temporary vector used in various computations.
+   REAL(ReKi)                   :: GBoxTrq                                         ! Gearbox torque on the LSS side in N-m (calculated from inputs and parameters).
+   REAL(ReKi)                   :: GBoxEffFac2                                     ! A second gearbox efficiency factor = ( 1 / GBoxEff^SgnPrvLSTQ - 1 )
+
+   INTEGER(IntKi)               :: I                                               ! Loops through some or all of the DOFs
+   INTEGER(IntKi)               :: J                                               ! Counter for elements
+   INTEGER(IntKi)               :: K                                               ! Counter for blades
+   INTEGER(IntKi)               :: L                                               ! Generic index
+
+   ! Initialize the output matrices:
+   GlobalStiffnessMatrix = 0.0
+   ExternalForceVector   = 0.0
+
+   ! Calculate static gearbox torque (HSSBrTrq_static is assumed to be a static value)
+   GBoxTrq = ( u_static%GenTrq + HSSBrTrq_static ) * ABS(p%GBRatio) 
+   
+   ! NOTE: For this static function, we assume that all acceleration-related fields
+   !       in RtHSdat_static (e.g., LinAccESt, LinAccEOt, AngAccEGt) are already zero.
+   !       Also, force/moment fields (e.g., FSAero, MMAero, FTHydrot, MFHydrot, FrcZAllt, MomXAllt, FrcONcRtt, MomBNcRtt, MomNGnRtt)
+   !       are assumed to represent only their static (non-inertial, non-damping) components.
+
+   ! Loop through all blades
+   DO K = 1,p%NumBl 
+      ! Blade tip forces (only static aero and gravity contributions)
+      ! Original: TmpVec1 = RtHSdat%FSTipDrag(:,K) - p%TipMass(K)*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccESt(:,K,p%TipNode) )
+      TmpVec1 = RtHSdat_static%FSTipDrag(:,K) - p%TipMass(K)*( p%Gravity*CoordSys%z2 ) ! LinAccESt is 0
+      DO I = 1,p%DOFs%NPSBE(K)
+            ExternalForceVector(p%DOFs%PSBE(K,I)) = ExternalForceVector(p%DOFs%PSBE(K,I)) + &
+                                                     DOT_PRODUCT( RtHSdat_static%PLinVelES(K,p%TipNode,p%DOFs%PSBE(K,I),0,:), TmpVec1 )
+      ENDDO
+
+      ! Loop through blade nodes / elements
+      DO J = 1,p%BldNodes 
+         ! Blade element forces (static aero, gravity, hydro)
+         ! Original: TmpVec1 = RtHSdat%FSAero(:,K,J)*p%DRNodes(J) - p%BElmntMass(J,K)*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccESt(:,K,J) )
+         TmpVec1 = ( RtHSdat_static%FSAero(:,K,J)*p%DRNodes(J) ) - p%BElmntMass(J,K)*( p%Gravity*CoordSys%z2 ) ! LinAccESt is 0
+         TmpVec3 = ( RtHSdat_static%MMAero(:,K,J)*p%DRNodes(J) )
+         DO I = 1,p%DOFs%NPSBE(K)
+               ExternalForceVector(p%DOFs%PSBE(K,I)) = ExternalForceVector(p%DOFs%PSBE(K,I)) + &
+                                                         DOT_PRODUCT( RtHSdat_static%PLinVelES(K,J,p%DOFs%PSBE(K,I),0,:), TmpVec1 ) + &
+                                                         DOT_PRODUCT( RtHSdat_static%PAngVelEM(K,J,p%DOFs%PSBE(K,I),0,:), TmpVec3 )
+         ENDDO
+      ENDDO ! J - Blade nodes / elements      
+      
+      ! Blade elasticity forces (damping forces are zero in static analysis)
+      IF ( p%DOF_Flag(DOF_BF(K,1)) )  THEN
+         GlobalStiffnessMatrix(DOF_BF(K,1),DOF_BF(K,1)) = GlobalStiffnessMatrix(DOF_BF(K,1),DOF_BF(K,1)) + p%KBF(K,1,1)
+         GlobalStiffnessMatrix(DOF_BF(K,1),DOF_BF(K,2)) = GlobalStiffnessMatrix(DOF_BF(K,1),DOF_BF(K,2)) + p%KBF(K,1,2)
+      ENDIF
+      IF ( p%DOF_Flag(DOF_BF(K,2)) )  THEN
+         GlobalStiffnessMatrix(DOF_BF(K,2),DOF_BF(K,1)) = GlobalStiffnessMatrix(DOF_BF(K,2),DOF_BF(K,1)) + p%KBF(K,2,1)
+         GlobalStiffnessMatrix(DOF_BF(K,2),DOF_BF(K,2)) = GlobalStiffnessMatrix(DOF_BF(K,2),DOF_BF(K,2)) + p%KBF(K,2,2)
+      ENDIF
+      IF ( p%DOF_Flag(DOF_BE(K,1)) )  THEN
+         GlobalStiffnessMatrix(DOF_BE(K,1),DOF_BE(K,1)) = GlobalStiffnessMatrix(DOF_BE(K,1),DOF_BE(K,1)) + p%KBE(K,1,1)
+      ENDIF
+   END DO ! K - Blades
+         
+   ! Yaw bearing gravity force
+   ! Original: TmpVec1 = -p%YawBrMass*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccEOt )
+   TmpVec1 = -p%YawBrMass*( p%Gravity*CoordSys%z2 ) ! LinAccEOt is 0
+   DO I = 1,p%DOFs%NPTTE
+         ExternalForceVector(p%DOFs%PTTE(I)) = ExternalForceVector(p%DOFs%PTTE(I)) + &
+                                                DOT_PRODUCT( RtHSdat_static%PLinVelEO(p%DOFs%PTTE(I),0,:), TmpVec1 )
+   ENDDO
+   
+   ! Loop through tower nodes / elements
+   DO J = 1,p%TwrNodes
+      ! Tower element forces (static hydro, gravity)
+      ! Original: TmpVec1 = ( RtHSdat%FTHydrot(:,J) )*abs(p%DHNodes(J)) - p%TElmntMass(J)*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccETt(:,J) )
+      TmpVec1 = ( RtHSdat_static%FTHydrot(:,J) )*abs(p%DHNodes(J)) - p%TElmntMass(J)*( p%Gravity*CoordSys%z2 ) ! LinAccETt is 0
+      TmpVec3 = ( RtHSdat_static%MFHydrot(:,J) )*abs(p%DHNodes(J))
+      DO I = 1,p%DOFs%NPTTE
+            ExternalForceVector(p%DOFs%PTTE(I)) = ExternalForceVector(p%DOFs%PTTE(I)) + &
+                                                  DOT_PRODUCT( RtHSdat_static%PLinVelET(J,p%DOFs%PTTE(I),0,:), TmpVec1 ) + &
+                                                  DOT_PRODUCT( RtHSdat_static%PAngVelEF(J,p%DOFs%PTTE(I),0,:), TmpVec3 )
+      ENDDO
+   ENDDO ! J - Tower nodes / elements   
+   
+   ! Tower elasticity forces (damping forces are zero in static analysis)
+   IF ( p%DOF_Flag(DOF_TFA1) )  THEN
+      GlobalStiffnessMatrix(DOF_TFA1,DOF_TFA1) = GlobalStiffnessMatrix(DOF_TFA1,DOF_TFA1) + p%KTFA(1,1)
+      GlobalStiffnessMatrix(DOF_TFA1,DOF_TFA2) = GlobalStiffnessMatrix(DOF_TFA1,DOF_TFA2) + p%KTFA(1,2)
+   ENDIF
+   IF ( p%DOF_Flag(DOF_TSS1) )  THEN
+      GlobalStiffnessMatrix(DOF_TSS1,DOF_TSS1) = GlobalStiffnessMatrix(DOF_TSS1,DOF_TSS1) + p%KTSS(1,1)
+      GlobalStiffnessMatrix(DOF_TSS1,DOF_TSS2) = GlobalStiffnessMatrix(DOF_TSS1,DOF_TSS2) + p%KTSS(1,2)
+   ENDIF
+   IF ( p%DOF_Flag(DOF_TFA2) )  THEN
+      GlobalStiffnessMatrix(DOF_TFA2,DOF_TFA1) = GlobalStiffnessMatrix(DOF_TFA2,DOF_TFA1) + p%KTFA(2,1)
+      GlobalStiffnessMatrix(DOF_TFA2,DOF_TFA2) = GlobalStiffnessMatrix(DOF_TFA2,DOF_TFA2) + p%KTFA(2,2)
+   ENDIF
+   IF ( p%DOF_Flag(DOF_TSS2) )  THEN
+      GlobalStiffnessMatrix(DOF_TSS2,DOF_TSS1) = GlobalStiffnessMatrix(DOF_TSS2,DOF_TSS1) + p%KTSS(2,1)
+      GlobalStiffnessMatrix(DOF_TSS2,DOF_TSS2) = GlobalStiffnessMatrix(DOF_TSS2,DOF_TSS2) + p%KTSS(2,2)
+   ENDIF
+   
+   ! Other DOF forces (only static components)
+   ! Mass matrix contributions from partial loads are removed.
+   ! Coriolis/centrifugal terms are removed.
+   IF ( p%DOF_Flag (DOF_Sg  ) )  THEN
+         ExternalForceVector(DOF_Sg  ) = ExternalForceVector(DOF_Sg  ) + DOT_PRODUCT( RtHSdat_static%PLinVelEZ(DOF_Sg ,0,:), RtHSdat_static%FrcZAllt )
+   ENDIF
+   IF ( p%DOF_Flag (DOF_Sw  ) )  THEN
+         ExternalForceVector(DOF_Sw  ) = ExternalForceVector(DOF_Sw  ) + DOT_PRODUCT( RtHSdat_static%PLinVelEZ(DOF_Sw ,0,:), RtHSdat_static%FrcZAllt )
+   ENDIF
+   IF ( p%DOF_Flag (DOF_Hv  ) )  THEN
+         ExternalForceVector(DOF_Hv  ) = ExternalForceVector(DOF_Hv  ) + DOT_PRODUCT( RtHSdat_static%PLinVelEZ(DOF_Hv ,0,:), RtHSdat_static%FrcZAllt )
+   ENDIF
+   IF ( p%DOF_Flag (DOF_R   ) )  THEN
+         ExternalForceVector(DOF_R   ) = ExternalForceVector(DOF_R   ) + DOT_PRODUCT( RtHSdat_static%PAngVelEX(DOF_R  ,0,:), RtHSdat_static%MomXAllt )
+   ENDIF
+   IF ( p%DOF_Flag (DOF_P   ) )  THEN
+         ExternalForceVector(DOF_P   ) = ExternalForceVector(DOF_P   ) + DOT_PRODUCT( RtHSdat_static%PAngVelEX(DOF_P  ,0,:), RtHSdat_static%MomXAllt )
+   ENDIF
+   IF ( p%DOF_Flag (DOF_Y   ) )  THEN
+         ExternalForceVector(DOF_Y   ) = ExternalForceVector(DOF_Y   ) + DOT_PRODUCT( RtHSdat_static%PAngVelEX(DOF_Y  ,0,:), RtHSdat_static%MomXAllt ) + u_static%YawMom
+   ENDIF
+
+   ! Coupling terms for TFA/TSS DOFs (only static force contributions)
+   IF ( p%DOF_Flag (DOF_TFA1) )  THEN
+         ExternalForceVector(DOF_TFA1) = ExternalForceVector(DOF_TFA1) + &
+                                          DOT_PRODUCT( RtHSdat_static%PLinVelEO(DOF_TFA1,0,:), RtHSdat_static%FrcONcRtt  ) + &
+                                          DOT_PRODUCT( RtHSdat_static%PAngVelEB(DOF_TFA1,0,:), RtHSdat_static%MomBNcRtt  )
+   ENDIF
+   IF ( p%DOF_Flag (DOF_TSS1) )  THEN
+         ExternalForceVector(DOF_TSS1) = ExternalForceVector(DOF_TSS1) + &
+                                          DOT_PRODUCT( RtHSdat_static%PLinVelEO(DOF_TSS1,0,:), RtHSdat_static%FrcONcRtt  ) + &
+                                          DOT_PRODUCT( RtHSdat_static%PAngVelEB(DOF_TSS1,0,:), RtHSdat_static%MomBNcRtt  )
+   ENDIF
+   IF ( p%DOF_Flag (DOF_TFA2) )  THEN
+         ExternalForceVector(DOF_TFA2) = ExternalForceVector(DOF_TFA2) + &
+                                          DOT_PRODUCT( RtHSdat_static%PLinVelEO(DOF_TFA2,0,:), RtHSdat_static%FrcONcRtt  ) + &
+                                          DOT_PRODUCT( RtHSdat_static%PAngVelEB(DOF_TFA2,0,:), RtHSdat_static%MomBNcRtt  )
+   ENDIF
+   IF ( p%DOF_Flag (DOF_TSS2) )  THEN
+         ExternalForceVector(DOF_TSS2) = ExternalForceVector(DOF_TSS2) + &
+                                          DOT_PRODUCT( RtHSdat_static%PLinVelEO(DOF_TSS2,0,:), RtHSdat_static%FrcONcRtt  ) + &
+                                          DOT_PRODUCT( RtHSdat_static%PAngVelEB(DOF_TSS2,0,:), RtHSdat_static%MomBNcRtt  )
+   ENDIF
+   
+   ! Yaw, RFrl, GeAz, DrTr, TFrl, Teet DOFs (only static force/moment contributions)
+   ! TmpVec related to generator inertia is effectively zero for static.
+   IF ( p%DOF_Flag (DOF_GeAz) )  THEN
+         ExternalForceVector(DOF_GeAz) = ExternalForceVector(DOF_GeAz) + &
+                                          DOT_PRODUCT( RtHSdat_static%PAngVelEL(DOF_GeAz,0,:), RtHSdat_static%MomLPRott ) - GBoxTrq
+   ENDIF
+   IF ( p%DOF_Flag (DOF_DrTr) )  THEN
+         ExternalForceVector(DOF_DrTr) = ExternalForceVector(DOF_DrTr) + &
+                                          DOT_PRODUCT( RtHSdat_static%PAngVelEL(DOF_DrTr,0,:), RtHSdat_static%MomLPRott ) - p%DTTorSpr*x_static%QT(DOF_DrTr)
+         GlobalStiffnessMatrix(DOF_DrTr, DOF_DrTr) = GlobalStiffnessMatrix(DOF_DrTr, DOF_DrTr) + p%DTTorSpr
+   ENDIF
+   IF ( p%DOF_Flag (DOF_TFrl) )  THEN
+         ExternalForceVector(DOF_TFrl) = ExternalForceVector(DOF_TFrl) + &
+                                          DOT_PRODUCT( RtHSdat_static%PAngVelEA(DOF_TFrl,0,:), RtHSdat_static%MomNTailt ) + RtHSdat_static%TFrlMom
+   ENDIF
+   IF ( ( p%NumBl == 2 ) ) THEN 
+      IF ( p%DOF_Flag(DOF_Teet) )  THEN
+         ExternalForceVector(DOF_Teet) = ExternalForceVector(DOF_Teet) + &
+                                          DOT_PRODUCT( RtHSdat_static%PAngVelEH(DOF_Teet,0,:), RtHSdat_static%MomLPRott ) + RtHSdat_static%TeetMom
+      ENDIF
+   ENDIF
+
+   ! Fill in the upper triangle of the stiffness matrix (it should be symmetric)
+   DO L = 2,p%DOFs%NActvDOF
+      DO I = 1,L-1
+         GlobalStiffnessMatrix(p%DOFs%SrtPS(I),p%DOFs%SrtPS(L)) = GlobalStiffnessMatrix(p%DOFs%SrtPS(L),p%DOFs%SrtPS(I))
+      ENDDO
+   ENDDO
+
+   ! Gearbox friction (non-conservative, contributes to F_ext)
+   GBoxEffFac2 = ( 1.0/RtHSdat_static%GBoxEffFac - 1.0 )
+   ExternalForceVector(DOF_GeAz) = ExternalForceVector(DOF_GeAz) - GBoxEffFac2 * GBoxTrq ! AngAccEGt is 0
+
+END SUBROUTINE fill_steady_augmat
+```
+
+### 3. 使用 LU 分解求解 `q` (在主程序中调用)
+
+在你的主程序中，调用 `fill_steady_augmat` 后，你需要一个线性代数库（如 LAPACK）来执行 LU 分解并求解 `[K]q = F_ext`。
+
+```fortran
+PROGRAM StaticAnalysis
+  USE ED_Types ! Assuming ED_Types contains definitions for ED_ParameterType, ED_ContinuousStateType, etc.
+  ! USE LAPACK_MODULE ! Assuming you have a LAPACK interface module
+
+  IMPLICIT NONE
+
+  ! Declare variables (matching the subroutine interface)
+  TYPE(ED_ParameterType)        :: p
+  TYPE(ED_ContinuousStateType)  :: x_static
+  TYPE(ED_CoordSys)             :: CoordSys
+  TYPE(ED_InputType)            :: u_static
+  TYPE(ED_RtHndSide)            :: RtHSdat_static
+  REAL(ReKi)                    :: HSSBrTrq_static
+  
+  INTEGER(IntKi)                :: N_DOF
+  REAL(R8Ki), ALLOCATABLE       :: GlobalStiffnessMatrix(:,:)
+  REAL(R8Ki), ALLOCATABLE       :: ExternalForceVector(:)
+  REAL(R8Ki), ALLOCATABLE       :: q_solution(:)
+  INTEGER(IntKi), ALLOCATABLE   :: ipiv(:)
+  INTEGER(IntKi)                :: info
+
+  ! --- 1. 初始化参数和静态状态 ---
+  ! 这一步需要根据你的具体风机模型和稳态工况来填充
+  ! 例如：
+  ! CALL Init_Parameters(p)
+  ! CALL Init_StaticStates(x_static) ! x_static%QT 应该包含初始猜测的位移，x_static%QDT 应该为0
+  ! CALL Init_CoordSys(CoordSys)
+  ! CALL Init_StaticInputs(u_static) ! u_static%GenTrq 等应为稳态值
+  ! HSSBrTrq_static = 0.0_ReKi ! 或其他稳态刹车力
+  ! CALL Init_RtHSdat_Static(RtHSdat_static, p, x_static, CoordSys, u_static) 
+  !    ! 这个函数需要你自己实现，确保 RtHSdat_static 中的加速度项为0，力/力矩项为稳态值
+
+  ! 获取自由度数量
+  N_DOF = p%DOFs%NActvDOF ! 假设 p%DOFs%NActvDOF 包含了所有激活的自由度数量
+
+  ! 分配矩阵和向量
+  ALLOCATE(GlobalStiffnessMatrix(N_DOF, N_DOF))
+  ALLOCATE(ExternalForceVector(N_DOF))
+  ALLOCATE(q_solution(N_DOF))
+  ALLOCATE(ipiv(N_DOF))
+
+  ! --- 2. 调用 fill_steady_augmat 构建 [K] 和 F_ext ---
+  CALL fill_steady_augmat( p, x_static, CoordSys, u_static, HSSBrTrq_static, RtHSdat_static, &
+                           GlobalStiffnessMatrix, ExternalForceVector )
+
+  ! --- 3. 求解线性方程组 [K]q = F_ext ---
+  ! LAPACK 的 DGESV 子程序会修改 GlobalStiffnessMatrix (存储 LU 分解结果) 和 ExternalForceVector (存储解)
+  q_solution = ExternalForceVector ! 将 F_ext 复制到 q_solution，因为 DGESV 会覆盖它
+
+  ! 调用 LAPACK 的 DGESV 子程序
+  ! N: 矩阵的阶数
+  ! NRHS: 右端向量的数量 (这里是1)
+  ! A: 输入矩阵 [K]，输出为 LU 分解结果
+  ! LDA: 矩阵 A 的第一维大小
+  ! IPIV: 存储行交换信息的向量
+  ! B: 输入向量 F_ext，输出为解 q
+  ! LDB: 矩阵 B 的第一维大小 (这里是 N_DOF)
+  ! INFO: 错误信息 (0表示成功)
+  CALL DGESV( N_DOF, 1, GlobalStiffnessMatrix, N_DOF, ipiv, q_solution, N_DOF, info )
+
+  IF ( info == 0 ) THEN
+    PRINT *, 'Static displacements (q_solution):'
+    PRINT *, q_solution
+    ! 可以将 q_solution 更新回 x_static%QT，以便后续使用
+    ! x_static%QT = q_solution
+  ELSE
+    PRINT *, 'Error solving linear system: INFO = ', info
+  END IF
+
+  ! 释放内存
+  DEALLOCATE(GlobalStiffnessMatrix)
+  DEALLOCATE(ExternalForceVector)
+  DEALLOCATE(q_solution)
+  DEALLOCATE(ipiv)
+
+END PROGRAM StaticAnalysis
+```
+
+**重要提示**:
+
+*   **`Init_RtHSdat_Static` 函数**: 这个函数是关键。在实际的 FAST 仿真中，`RtHSdat` 是由 `RtHndSid` 模块在每个时间步计算的。对于静态分析，你需要确保 `RtHSdat_static` 中的所有**加速度**相关项（如 `LinAccESt`, `LinAccEOt`, `AngAccEGt`）都被设置为零。同时，所有**力/力矩**相关项（如 `FSAero`, `MMAero`, `FTHydrot`, `MFHydrot`, `FrcZAllt`, `MomXAllt`, `FrcONcRtt`, `MomBNcRtt`, `MomNGnRtt`）必须准确反映**稳态**的外部载荷（例如，恒定风速下的气动力，或静水中的水动力）。
+*   **`x_static%QT`**: 在调用 `DGESV` 之前，`x_static%QT` 应该包含一个初始猜测的位移向量。如果 `[K]` 是非线性的（例如，包含几何刚度效应），你可能需要一个迭代过程（如牛顿-拉普森法）来求解，每次迭代都更新 `x_static%QT` 并重新构建 `[K]` 和 `F_ext`。这里提供的代码是针对线性 `[K]` 的直接求解。
+*   **LAPACK 接口**: `DGESV` 是 LAPACK 库中的一个标准双精度通用矩阵求解器。你需要确保你的 Fortran 编译器能够链接到 LAPACK 库，并且 `DGESV` 的接口声明正确。
+*   **DOF 索引**: `p%DOFs%PSBE`, `p%DOFs%PTTE`, `p%DOFs%SrtPS` 等数组用于将物理自由度映射到矩阵的行/列索引。这些在 `FillAugMat` 中已经定义，并应在 `fill_steady_augmat` 中保持一致。
+
+通过这些步骤，你就可以将原始的动态 `FillAugMat` 函数转化为一个高效的静态分析工具，并利用 LU 分解来求解塔架在稳态条件下的变形。
+[Timestamp: 2025/11/03 15:18:00]
+
+**user**: 求解稳态的情况下，我的augmat矩阵是什么形式，最后一列前的矩阵是只有对角线元素？
+[Timestamp: 2025/11/03 15:21:04]
+
+**ai**: 在稳态（静态）情况下，`AugMat` 矩阵的形式会发生根本性的变化。
+
+原始的 `AugMat` 实际上是一个增广矩阵，代表了凯恩方程的线性化形式：
+`[ [M] | [C_coriolis] ] * [ [ddot_q] ] = [ [-K] | [-D] | [F_ext] ] * [ [q] ]`
+或者更常见的形式是：
+`[M]ddot_q + [D]dot_q + [K]q = F_ext`
+
+其中 `AugMat` 的结构是 `[ [M_eff] | [F_eff] ]`，其中 `M_eff` 是有效质量矩阵，`F_eff` 是包含所有非惯性力和速度相关惯性力的有效力向量。
+
+### 稳态情况下的 `AugMat` 形式
+
+在稳态（静态）情况下，我们有以下假设：
+*   广义加速度 $\ddot{q} = 0$
+*   广义速度 $\dot{q} = 0$
+
+将这些假设代入动力学方程，所有与 $\ddot{q}$ 和 $\dot{q}$ 相关的项都将变为零。因此，方程简化为：
+`[K]q = F_ext`
+
+在这种情况下，我们不再需要构建一个包含质量矩阵 `[M]` 的增广矩阵。相反，我们只需要构建一个**刚度矩阵 `[K]`** 和一个**外部载荷向量 `F_ext`**。
+
+所以，在 `fill_steady_augmat` 函数中，`AugMat` 的概念会转化为两个独立的输出：
+1.  `GlobalStiffnessMatrix` (对应 `[K]`)
+2.  `ExternalForceVector` (对应 `F_ext`)
+
+**`GlobalStiffnessMatrix` 的形式：**
+
+这个矩阵**通常不是只有对角线元素**。它是一个**全矩阵**（或至少是一个稀疏矩阵，但包含非对角线元素）。
+
+原因如下：
+
+*   **模态耦合**: 塔筒的弯曲振动（例如，前后弯曲 TFA1 和 TFA2，侧向弯曲 TSS1 和 TSS2）通常是耦合的。这意味着一个自由度的位移会影响另一个自由度的弹性力。在 [[Tower]] 笔记的 (Eq. 6) 和 (Eq. 7, 9) 中，你可以看到广义刚度 $k_{ij}$ 包含非对角线项（例如 $k_{12}^{'TFA}$），这表示不同模态之间存在弹性耦合。
+    ```
+    F_{r \Big| ElasticT} = \left\{
+    \begin{array}{ll}
+    -k_{11}^{'TFA} q_{TFA1} - k_{12}^{'TFA} q_{TFA2} & \text{for } r=TFA1 \\
+    ...
+    -k_{21}^{'TFA} q_{TFA1} - k_{22}^{'TFA} q_{TFA2} & \text{for } r=TFA2 \\
+    ...
+    \end{array}
+    \right.
+    ```
+    这些非对角线项 $k_{12}$ 和 $k_{21}$ 会直接填充到 `GlobalStiffnessMatrix` 的非对角线位置。
+
+*   **结构连接**: 即使是不同的物理部件（如塔筒和机舱），它们的自由度也可能通过连接点产生耦合。例如，塔筒的弯曲会引起塔顶的位移和转动，进而影响机舱和叶片的平衡。这些耦合效应也会导致 `GlobalStiffnessMatrix` 出现非对角线元素。
+
+**`ExternalForceVector` 的形式：**
+
+这是一个列向量，包含了所有作用在系统上的稳态广义力，例如：
+*   **重力**：塔筒自身质量和塔顶质量的重力贡献。
+*   **气动力**：稳态风载荷作用在塔筒和叶片上的气动力和力矩。
+*   **水动力**：如果考虑海上风机，水动力也会有静态分量。
+*   **其他静态载荷**：如预应力、刹车力（如果视为静态载荷）等。
+
+### 总结
+
+在稳态情况下，`fill_steady_augmat` 函数将构建一个：
+
+*   **`GlobalStiffnessMatrix` (`[K]`)**: 这是一个**非对角线矩阵**，包含了所有自由度之间的弹性耦合。
+*   **`ExternalForceVector` (`F_ext`)**: 这是一个列向量，包含了所有作用在系统上的稳态广义力。
+
+然后，通过求解 `[K]q = F_ext`，我们可以得到系统的稳态广义位移 `q`。
+
+因此，你不能假设 `GlobalStiffnessMatrix` 只有对角线元素。它是一个反映系统整体刚度特性的全矩阵。
+[Timestamp: 2025/11/03 15:21:14]

书籍/力学书籍/OpenFast/FASTKinetics/auto/Drivertrain.md

@@ -0,0 +1,310 @@
+
+The inertia of the drivetrain brings about generalized inertia forces and the load in the generator, high-speed shaft brake, gearbox (friction forces resulting from nonzero GBoxEff ) and the flexibility of the low speed shaft bring about generalized active forces.  Note that all of these equations assume that the rotor is spinning about the positive $c_{I}$ axis (they assume that the rotor can’t be forced to rotate in the opposite direction).  This model works for any gearbox arrangement (including no gearbox, single stage, or multi-stage) as long as the generator rotates about the shaft axis (it may not be skewed relative to the shaft, even though it may rotate in the opposite direction of the low-speed shaft due to the gearbox stages).  If there is no gearbox, simply set $G B R a t i o\,=\,G B o x E f f=G e n D i r=1$ (GB R $e\nu e r s e={\mathrm{False}}$ ).  
+传动系惯性带来广义惯性力和发电机负载、高速轴制动器、齿轮箱（由于GBoxEff非零导致的摩擦力）以及低速轴挠曲性，都会产生广义主动力。需要注意的是，这些方程都假设转子绕正$c_{I}$轴旋转（即假设转子不能被强制旋转到相反方向）。该模型适用于任何齿轮箱布置（包括无齿轮箱、单级或多级），只要发电机绕轴旋转即可（即使它可能由于齿轮箱级数而以低速轴相反的方向旋转，但不能相对于轴倾斜）。如果不存在齿轮箱，只需将$G B R a t i o\,=\,G B o x E f f=G e n D i r=1$ (GB R $e\nu e r s e={\mathrm{False}}$ )。
+
+The mechanical torque within the generator is applied to the high speed shaft and equally and oppositely to the structure that furls with the rotor as follows:  
+发电机内部的机械扭矩施加到高速轴上，并以相等且相反的方式作用于与转子一起摆动的结构，具体如下：
+
+$$
+\left.F_{r}\right|_{G e n}=\left(^{E}\omega_{r}^{G}-{^{E}\!\omega_{r}^{R}}\right)\cdot M_{G e n}^{G}\quad\left(r=l,2,...,22\right)
+$$  
+
+Thus,  
+
+$$
+\left.F_{r}\right|_{G e n}=\left\{\!\!\!\begin{array}{l l}{{^{E}\pmb{\omega}_{G e d z}^{G}\cdot M_{G e n}^{G}}}&{{f o r~~r=G e A z}}\\ {{0}}&{{o t h e r w i s e}}\end{array}\ \right.
+$$  
+
+$$
+\boldsymbol{M_{G e n}^{G}}=-G e n D i\boldsymbol{r}\cdot\boldsymbol{T^{G e n}}\left(G B R a t i o\cdot\dot{q}_{G e d z},t\right)\boldsymbol{c}_{I}
+$$  
+
+Note that a positive $T^{G e n}$ represents a load (positive power extracted) whereas a negative $T^{G e n}$ represents a motoring-up situation (negative power extracted, or power input).  Thus,  
+
+$$
+F_{r}\Big|_{G e n}=\left\{\begin{array}{l l}{\left(G e n D i r\cdot G B R a t i o\cdot c_{I}\right)\cdot\left[-G e n D i r\cdot T^{G e n}\left(G B R a t i o\cdot\dot{q}_{G e d z},t\right)c_{I}\right]}&{f o r\;\;\;r=G e n D i r}\\ {0}&{o t h e r w i s e}\end{array}\right.
+$$  
+
+Or since $G e n D i r^{2}=1$ ,  
+
+$$
+
+F_r \Big|_{Gen} = \begin{cases}
+
+-GRRatio \cdot T^{Gen} (GRRatio \cdot \dot{q}_{GeAz}, t) & \text{for } r = GeAz \\
+
+0 & \text{otherwise}
+
+\end{cases}
+
+$$
+ Thus, 
+
+$$
+[C(q,t)]|_{\text{Gen}} = 0
+$$
+
+$$
+\{-f(\dot{q},q,t)\}|_{\text{Gen}} = \left\{
+\begin{array}{l}
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+-GBRatio \cdot T^{\text{Gen}} (GBRatio \cdot \dot{q}_{\text{GeAz}}, t) \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots
+\end{array}
+\right\}
+$$
+
+Similarly, the mechanical torque applied to the high-speed shaft from the high-speed shaft brake is applied equally and oppositely to the structure that furls with the rotor.  Thus,  
+同样地，从高速轴制动器施加到高速轴上的机械扭矩，也会以相等且相反的方式作用于与转子一同回转的结构上。因此，
+
+$F_{r}\big|_{B r a k e}=\left\{\!\!\!\begin{array}{l l}{{^{E}\pmb{\omega}_{G e A z}^{G}\cdot M_{B r a k e}^{G}}}&{{f o r\;\;\;r=G e A z}}\\ {{0}}&{{o t h e r w i s e}}\end{array}\!\!\right.$ where $M_{B r a k e}^{G}=-G e n D i r\cdot T^{B r a k e}\left(t\right)c_{I}\qquad\mathrm{~and~where}\qquad T^{B r a k e}\left(t\right)=H S S B r k T\left(t\right)$  
+
+which is assumed to be positive in value always.  Thus,  
+
+$$
+F_{r}\bigr|_{B r a k e}=\left\{\begin{array}{l l}{-G B R a t i o\cdot T^{B r a k e}\left(t\right)}&{f o r\;\;\;r=G e A z}\\ {O}&{o t h e r w i s e}\end{array}\right.
+$$  
+
+Thus,  
+
+$$
+[C(q,t)]|_{\text{Brake}} = 0
+$$
+
+$$
+\{-f(\dot{q},q,t)\}|_{\text{Brake}} = \left\{
+\begin{array}{l}
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+-GBRatio \cdot T^{\text{Brake}}(t) \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots
+\end{array}
+\right\}
+$$
+
+If the translational inertia of the drivetrain is assumed to be incorporated into that of the structure that furls with the rotor, then the high-speed shaft generator inertia generalized force is as follows: 
+如果假设动能惯性由驱动系传递至随转子摆动的结构部分，则高速轴发电机惯性概括力如下： 
+
+$$
+F_r^*|_{\text{G}} = {}^E \omega_r^{\text{G}} \cdot (-\overline{\overline{{I}}}^{\text{G}} \cdot {}^E {a}^{\text{G}} - {}^E \boldsymbol{\omega}^{\text{G}} \times \overline{\overline{{I}}}^{\text{G}} \cdot {}^E \boldsymbol{\omega}^{\text{G}}) \quad \text{where} \quad \overline{\overline{{I}}}^{\text{G}} = \text{GenInerc}_{1c1} \tag{r=1,2,...,22}
+$$
+
+or, 
+$$
+F_r^*|_{\text{G}} = {}^E \boldsymbol{\omega}_r^{\text{G}} \cdot \left\{ -\overline{\overline{{I}}}^{\text{G}} \cdot \left\{ \left( \sum_{i=4}^{13} {}^E \omega_i^{\text{G}} \ddot{q}_i \right) + \left[ \sum_{i=7}^{13} \frac{d}{dt} \left( {}^E \boldsymbol{\omega}_i^{\text{G}} \right) \dot{q}_i \right] \right\} - {}^E \boldsymbol{\omega}^{\text{G}} \times \overline{\overline{{I}}}^{\text{G}} \cdot {}^E \boldsymbol{\omega}^{\text{G}} \right\} \tag{r=1,2,...,22}
+$$
+
+However, since $c_{1}\cdot\frac{d}{d t}\Big(\{^{E}\omega_{G e A z}^{G}\Big)\propto c_{1}\cdot(^{E}\omega^{R}\times c_{1}\Big)={^{E}\omega^{R}}\cdot\Big(c_{1}\times c_{1}\Big)=0$ (the first $c_{1}$ coming from $\overline{\overline{{I}}}^{\text{G}}$ ), this simplifies as follows:  
+
+$$
+F_r^*|_{\text{G}} = {}^E \boldsymbol{\omega}_r^{\text{G}} \cdot \left\{ -\overline{\overline{{I}}}^{\text{G}} \cdot \left\{ \left( \sum_{i=4}^{13} {}^E \boldsymbol{\omega}_i^{\text{G}} \ddot{q}_i \right) + \left[ \sum_{i=7}^{12} \frac{d}{dt} \left( {}^E \boldsymbol{\omega}_i^{\text{G}} \right) \dot{q}_i \right] \right\} - {}^E \boldsymbol{\omega}^{\text{G}} \times \overline{\overline{{I}}}^{\text{G}} \cdot {}^E \boldsymbol{\omega}^{\text{G}} \right\} \tag{r=1,2,...,22}
+$$Or, 
+$$
+F_r^*|_{G} = \begin{cases}
+- {}^E \boldsymbol{\omega}_r^{\text{R}} \cdot \overline{\overline{{I}}}^{\text{G}} \cdot \left\{ \left( \sum_{i=4}^{13} {}^E \boldsymbol{\omega}_i^{\text{G}} \ddot{q}_i \right) + \left[ \sum_{i=7}^{12} \frac{d}{dt} \left( {}^E \boldsymbol{\omega}_i^{\text{R}} \right) \dot{q}_i \right] \right\} - ({}^E \boldsymbol{\omega}_r^{\text{R}} \cdot ({}^E \boldsymbol{\omega}^{\text{G}} \times \overline{\overline{{I}}}^{\text{G}} \cdot {}^E \boldsymbol{\omega}^{\text{G}})) & \text{for } r = 4,5,...,12 \\
+- \text{GenDir} \cdot \text{GenIner} \cdot \text{GBRatio} \left\{ \left( \left( \sum_{i=4}^{12} {}^E \boldsymbol{\omega}_i^{\text{R}} \ddot{q}_i \right) + \left[ \sum_{i=7}^{12} \frac{d}{dt} \left( {}^E \boldsymbol{\omega}_i^{\text{R}} \right) \dot{q}_i \right] \right) \cdot c_1 - \text{GenIner} \cdot \text{GBRatio}^2 \cdot \ddot{q}_{\text{GeAz}} \right\} - \text{GenDir} \cdot \text{GBRatioc}_1 \cdot ({}^E \boldsymbol{\omega}^{\text{G}} \times \overline{\overline{{I}}}^{\text{G}} \cdot {}^E \boldsymbol{\omega}^{\text{G}}) & \text{for } r = \text{GeAz} \\
+0 & \text{otherwise}
+\end{cases}
+$$
+
+However since $c_{1}\cdot\left({\,^{E}}{\omega^{G}}\times c_{1}\,\right)={^{E}}{\omega^{G}}\cdot\left(c_{1}\times c_{1}\right)=0$ (the first $c_{I}$ coming from $\overline{{\overline{{I}}}}^{G}$ ), this simplifies again as follows:  
+
+$$
+F_r^*|_{G} = \begin{cases}
+- {}^E \boldsymbol{\omega}_r^{\text{R}} \cdot \overline{\overline{{I}}}^{\text{G}} \cdot \left\{ \left( \sum_{i=4}^{13} {}^E \boldsymbol{\omega}_i^{\text{G}} \ddot{q}_i \right) + \left[ \sum_{i=7}^{12} \frac{d}{dt} \left( {}^E \boldsymbol{\omega}_i^{\text{R}} \right) \dot{q}_i \right] \right\} - ({}^E \boldsymbol{\omega}_r^{\text{R}} \cdot ({}^E \boldsymbol{\omega}^{\text{G}} \times \overline{\overline{{I}}}^{\text{G}} \cdot {}^E \boldsymbol{\omega}^{\text{G}})) & \text{for } r = 4,5,...,12 \\
+- \text{GenDir} \cdot \text{GenIner} \cdot \text{GBRatio} \left\{ \left( \left( \sum_{i=4}^{12} {}^E \boldsymbol{\omega}_i^{\text{R}} \ddot{q}_i \right) + \left[ \sum_{i=7}^{12} \frac{d}{dt} \left( {}^E \boldsymbol{\omega}_i^{\text{R}} \right) \dot{q}_i \right] \right) \cdot c_1 - \text{GenIner} \cdot \text{GBRatio}^2 \cdot \ddot{q}_{\text{GeAz}} \right\} - \text{GenDir} \cdot \text{GBRatioc}_1 \cdot ({}^E \boldsymbol{\omega}^{\text{G}} \times \overline{\overline{{I}}}^{\text{G}} \cdot {}^E \boldsymbol{\omega}^{\text{G}}) & \text{for } r = \text{GeAz} \\
+0 & \text{otherwise}
+\end{cases}
+$$
+
+The terms associated with DOFs 4,5,…,12 represent the fact that the rate of change of angular momentum of the generator can be considered as an additional torque on the structure that furls with the rotor (i.e., in addition to the torques on the structure transmitted directly from the low-speed shaft).  
+与DOF 4、5、…、12相关的项，表示发电机角动量变化率可以被视为一个附加于随转子展开的结构上的扭矩（即，除了直接从低速轴传递到结构上的扭矩之外）。
+
+$$
+[C(q,t)]_G = \begin{cases}
+{}^E \boldsymbol{\omega}_{\text{Row}}^{\text{G}} \cdot \overline{\overline{{I}}}^{\text{G}} \cdot {}^E \boldsymbol{\omega}_{\text{Col}}^{\text{G}} & \text{for (Row, Col = 4,5,...,13)} \\
+0 & \text{otherwise}
+\end{cases}
+$$
+
+$$
+\{-f(\dot{q},q,t)\}_G = \begin{cases}
+- {}^E \boldsymbol{\omega}_{\text{Row}}^{\text{G}} \cdot \left\{ \overline{\overline{{I}}}^{\text{G}} \cdot \left[ \sum_{i=7}^{12} \frac{d}{dt} \left( {}^E \boldsymbol{\omega}_i^{\text{R}} \right) \dot{q}_i \right] + {}^E \boldsymbol{\omega}^{\text{G}} \times \overline{\overline{{I}}}^{\text{G}} \cdot {}^E \boldsymbol{\omega}^{\text{G}} \right\} & \text{for (Row = 4,5,...,12)} \\
+- {}^E \boldsymbol{\omega}_{\text{Row}}^{\text{G}} \cdot \overline{\overline{{I}}}^{\text{G}} \cdot \left[ \sum_{i=7}^{12} \frac{d}{dt} \left( {}^E \boldsymbol{\omega}_i^{\text{R}} \right) \dot{q}_i \right] & \text{for (Row = 13)} \\
+0 & \text{otherwise}
+\end{cases}
+$$
+
+$$
+F_r |_{\text{ElasticDrive}} = -\frac{\partial V_{\text{Drive}}}{\partial q_r} \quad (r = 1,2,...,22) \quad \text{where} \quad V_{\text{Drive}} = \frac{1}{2} \text{DTTorSpr} \cdot q_{\text{DrTr}}^2
+$$
+So, 
+$$
+F_r |_{\text{ElasticDrive}} = \begin{cases}
+-DTTorSpr \cdot q_{DrTr} & \text{for } r = DrTr \\
+0 & \text{otherwise}
+\end{cases}
+$$
+and likewise
+$$
+F_r |_{\text{DampDrive}} = \begin{cases}
+-DTTorDmp \cdot \dot{q}_{DrTr} & \text{for } r = DrTr \\
+0 & \text{otherwise}
+\end{cases}
+$$
+
+
+Thus, 
+
+$$
+[C(q,t)]_{\text{ElasticDrive}} = 0
+$$
+
+$$
+\{-f(\dot{q},q,t)\}_{\text{ElasticDrive}} = \begin{Bmatrix}
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+-DTTorSpr \cdot q_{DrTr} \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\end{Bmatrix}
+$$
+
+$$
+[C(q,t)]_{\text{DampDrive}} = 0
+$$
+
+$$
+\{-f(\dot{q},q,t)\}_{\text{DampDrive}} = \begin{Bmatrix}
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+-DTTorDmp \cdot \dot{q}_{DrTr} \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\end{Bmatrix}
+$$
+
+Similar to the generator and HSS brake, the mechanical friction torque applied to the high speed shaft is applied equally and oppositely to the structure that furls with the rotor.  Thus,  
+类似于发电机和HSS制动器，施加到高速轴上的机械摩擦扭矩，以相等且相反的方式作用于与转子一同旋转的结构上。因此，
+
+$$
+F_r |_{\text{GBFric}} = \begin{cases}
+{}^E\omega_{\text{GeAz}}^G \cdot M_{\text{GBFric}}^G & \text{for } r = \text{GeAz} \\
+0 & \text{otherwise}
+\end{cases}
+\quad \text{where} \quad
+M_{\text{GBFric}}^G = \frac{T_{\text{GBFric}}^G(\ddot{q}, \dot{q}, q, t) c_1}{\text{GBRatio} \cdot \text{GenDir}}
+$$
+
+where, from a free-body diagram of the high and low-speed shafts, it is easily seen that the friction torque applied on the LSS upon the gearbox entrance, $T^{G B F r i c}\left(\ddot{q},\dot{q},q,t\right)$ , is always positive in value and equal to: 
+从高速轴和低速轴的受力图上可以清楚地看到，作用于低速轴进入变速箱处的摩擦转矩，$T^{G B F r i c}\left(\ddot{q},\dot{q},q,t\right)$，其值始终为正，且等于：
+
+$$
+T^{\text{GBFric}}(\ddot{q}, \dot{q}, q, t) = \left[ \text{GenIner} \cdot \text{GBRatio}^2 \cdot \ddot{q}_{\text{GeAz}} + \text{GenDir} \cdot \text{GenIner} \cdot \text{GBRatio} \cdot a^E \cdot c_1^R + \text{GBRatio} \cdot T^{\text{Gen}}(\text{GBRatio} \cdot \dot{q}_{\text{GeAz}}, t) + \text{GBRatio} \cdot T^{\text{Brake}}(t) \right] \cdot \left[ \frac{1}{\text{GBBoxEff}^{\text{SIGN(LSShftTq)}}} - 1 \right]
+$$
+Thus, 
+$$
+F_r |_{\text{GBFric}} = \begin{cases}
+-T^{\text{GBFric}}(\ddot{q}, \dot{q}, q, t) & \text{for } r = \text{GeAz} \\
+0 & \text{otherwise}
+\end{cases}
+$$
+
+$$
+F_r |_{\text{GBFric}} = \begin{cases}
+- \left[ \text{GenIner} \cdot \text{GBRatio}^2 \cdot \ddot{q}_{\text{GeAz}} + \text{GenDir} \cdot \text{GenIner} \cdot \text{GBRatio}^E a^R \cdot c_1 + \text{GBRatio} \cdot T^{\text{Gen}}(\text{GBRatio} \cdot \dot{q}_{\text{GeAz}}, t) + \text{GBRatio} \cdot T^{\text{Brake}}(t) \right] \cdot \left[ \frac{1}{\text{GBoxEff}^{\text{SIGN(LSShftTq)}}} - 1 \right] & \text{for } r = \text{GeAz} \\
+0 & \text{otherwise}
+\end{cases}
+$$
+
+Or, 
+$$
+F_r |_{\text{GBFric}} = \begin{cases}
+- \left[ \text{GenIner} \cdot \text{GBRatio}^2 \cdot \ddot{q}_{\text{GeAz}} \right. \\
+\quad \left. + \text{GenDir} \cdot \text{GenIner} \cdot \text{GBRatio} \left\{ \sum_{i=4}^{12} {^E \omega_i^R \ddot{q}_i} + \left[ \sum_{i=7}^{12} \frac{d}{dt} ({^E \omega_i^R}) \dot{q}_i \right] \right\} \cdot c_1 \right. \\
+\quad \left. + \text{GBRatio} \cdot T^{\text{Gen}}(\text{GBRatio} \cdot \dot{q}_{\text{GeAz}}, t) + \text{GBRatio} \cdot T^{\text{Brake}}(t) \right] \cdot \left[ \frac{1}{\text{GBoxEff}^{\text{SIGN(LSShftTq)}}} - 1 \right] & \text{for } r = \text{GeAz} \\
+0 & \text{otherwise}
+\end{cases}
+$$
+
+Thus, 
+$$
+[C(q,t)]_{\text{GBFric}} = \begin{cases}
+{^E\omega_{\text{Row}}^G \cdot \overline{\overline{{I}}}^{\text{G}} \cdot {^E\omega_{\text{Col}}^G}} \left[ \frac{1}{\text{GBoxEff}^{\text{SIGN(LSShftTq)}}} - 1 \right] & \text{for } (\text{Row}=13, \text{Col}=4,5,...,13) \\
+0 & \text{otherwise}
+\end{cases}
+$$
+
+$$
+\{-f(\dot{q},q,t)\}_{\text{GBFric}} = \begin{cases}
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+- \left\{ {^E\omega_{\text{GeAz}}^G \cdot \overline{\overline{{I}}}^{\text{G}} \cdot \left[ \sum_{i=7}^{12} \frac{d}{dt} ({^E\omega_i^R}) \dot{q}_i \right]} + \text{GBRatio} \cdot T^{\text{Gen}}(\text{GBRatio} \cdot \dot{q}_{\text{GeAz}}, t) + \text{GBRatio} \cdot T^{\text{Brake}}(t) \right\} \cdot \left[ \frac{1}{\text{GBoxEff}^{\text{SIGN(LSShftTq)}}} - 1 \right] \\
+\ldots \\
+\ldots \\
+\ldots \\
+\ldots \\
+\end{cases}
+$$
+
+It is noted that the gearbox friction generalized active force effects both the mass matrix and the forcing function.  Its effect on the mass matrix can be thought of as an apparent mass in the system (i.e., an active friction force which is a function of the generator acceleration).  The gearbox friction causes the mass matrix to become unsymmetric.  Note that the equation for DOF GeAz is an implicit equation (since the gearbox friction is a function of DOF GeAz), which has no closed-form solution.  To avoid the complications resulting from this implicit characterization, the value of the LSShftTq  from the previous time step is used in the SIGN() function in place of the value of the LSShftTq  in the current time step.  This may be done since it may be assumed that LSShftTq  will not change much between time steps if the time step is considered small enough.  Note that gearbox friction is the only term in the equations of motion that cause the mass matrix to be unsymmetrical.  The mass matrix will only be unsymmetric if GBoxEff  is not $100\%$ —if GBoxEff  is $100\%$ , then the mass matrix is completely symmetric.  
+需要注意的是，齿箱摩擦产生的广义主动力效应既影响质量矩阵，也影响强制函数。它对质量矩阵的影响可以被认为是系统中的视质量（即，一个随发电机加速度变化的活动摩擦力）。齿箱摩擦导致质量矩阵变得非对称。需要注意的是，DOF GeAz 的方程是一个隐式方程（因为齿箱摩擦是 DOF GeAz 的函数），因此没有闭合解。为了避免由此隐式特性带来的复杂性，在 SIGN() 函数中，使用上一个时间步长的 LSShftTq 值代替当前时间步长的 LSShftTq 值。这是可以做到的，因为如果时间步长足够小，可以假设 LSShftTq 在时间步长之间不会发生太大变化。需要注意的是，齿箱摩擦是导致质量矩阵非对称的唯一项。只有当 GBoxEff 不是 100% 时，质量矩阵才会是非对称的——如果 GBoxEff 是 100%，那么质量矩阵将完全对称。

书籍/力学书籍/OpenFast/FASTKinetics/auto/FillAugmat.md

@@ -0,0 +1,429 @@
+
+```
+SUBROUTINE FillAugMat( p, x, CoordSys, u, HSSBrTrq, RtHSdat, AugMat )
+!..................................................................................................................................
+
+      ! Passed variables
+   TYPE(ED_ParameterType),       INTENT(IN   )  :: p           !< Parameters
+   TYPE(ED_ContinuousStateType), INTENT(IN   )  :: x           !< Continuous states at Time
+   TYPE(ED_CoordSys),            INTENT(IN   )  :: CoordSys    !< The coordinate systems that have been set for these states/time
+   TYPE(ED_InputType),           INTENT(IN   )  :: u           !< The aero blade forces/moments
+   TYPE(ED_RtHndSide),           INTENT(IN   )  :: RtHSdat     !< data from the RtHndSid module (contains positions to be set)
+   REAL(ReKi),                   INTENT(IN )    :: HSSBrTrq    !<  SIGN( u%HSSBrTrqC, x%QDT(DOF_GeAz) ) or corrected value from FixHSS
+   REAL(R8Ki),                   INTENT(OUT)    :: AugMat(:,:) !< the return matrix 
+   
+      ! Local variables
+   REAL(ReKi)                   :: TmpVec    (3)                                   ! A temporary vector used in various computations.
+   REAL(ReKi)                   :: TmpVec1   (3)                                   ! A temporary vector used in various computations.
+   REAL(ReKi)                   :: TmpVec3   (3)                                   ! A temporary vector used in various computations.
+   REAL(ReKi)                   :: GBoxTrq                                         ! Gearbox torque on the LSS side in N-m (calculated from inputs and parameters).
+   REAL(ReKi)                   :: GBoxEffFac2                                     ! A second gearbox efficiency factor = ( 1 / GBoxEff^SgnPrvLSTQ - 1 )
+
+   INTEGER(IntKi)               :: I                                               ! Loops through some or all of the DOFs
+   INTEGER(IntKi)               :: J                                               ! Counter for elements
+   INTEGER(IntKi)               :: K                                               ! Counter for blades
+   INTEGER(IntKi)               :: L                                               ! Generic index
+
+   
+      ! Initialize the matrix:
+      
+   AugMat      = 0.0
+   GBoxTrq    = ( u%GenTrq + HSSBrTrq )*ABS(p%GBRatio) ! bjj: do we use HSSBrTrqC or HSSBrTrq?
+   
+   DO K = 1,p%NumBl ! Loop through all blades
+   
+
+      ! Initialize the portions of the mass matrix on and below the diagonal associated with purely blade DOFs (these portions can't
+      !   be calculated using partial loads) using the tip mass effects. 
+      ! Also, initialize the portions of the forcing vector associated with purely blade DOFs (these portions can't be calculated 
+      !   using partial loads) using the tip mass effects:
+      ! NOTE: The vector subscript array, PSBE(), used in the following loops must be sorted from smallest to largest DOF index in 
+      !       order for the loops to work to enter values only on and below the diagonal of AugMat():
+
+   
+      DO L = 1,p%DOFs%NPSBE(K)    ! Loop through all active (enabled) blade DOFs that contribute to the QD2T-related linear accelerations of the tip of blade K (point S(p%BldFlexL))
+         DO I = L,p%DOFs%NPSBE(K) ! Loop through all active (enabled) blade DOFs greater than or equal to L
+            AugMat(p%DOFs%PSBE(K,I),p%DOFs%PSBE(K,L)) = p%TipMass(K)*&
+                                        DOT_PRODUCT( RtHSdat%PLinVelES(K, p%TipNode, p%DOFs%PSBE(K,I),0,:), &   ! [C(q,t)]B
+                                                     RtHSdat%PLinVelES(K, p%TipNode, p%DOFs%PSBE(K,L),0,:)    )
+         ENDDO             ! I - All active (enabled) blade DOFs greater than or equal to L
+      ENDDO                ! L - All active (enabled) blade DOFs that contribute to the QD2T-related linear accelerations of the tip of blade K (point S(p%BldFlexL))
+
+      TmpVec1 = RtHSdat%FSTipDrag(:,K) - p%TipMass(K)*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccESt(:,K,p%TipNode) ) ! The portion of FrcS0Bt associated with the tip brake
+      DO I = 1,p%DOFs%NPSBE(K)    ! Loop through all active (enabled) blade DOFs that contribute to the QD2T-related linear accelerations of the tip of blade K (point S(p%BldFlexL))
+            AugMat(p%DOFs%PSBE(K,I), p%NAug) = DOT_PRODUCT( RtHSdat%PLinVelES(K,p%TipNode,p%DOFs%PSBE(K,I),0,:), &   ! {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB
+                                                              TmpVec1                               ) ! NOTE: TmpVec1 is still the portion of FrcS0Bt associated with the tip brake
+      ENDDO                ! I - All active (enabled) blade DOFs that contribute to the QD2T-related linear accelerations of the tip of blade K (point S(p%BldFlexL))
+   
+      
+
+      DO J = 1,p%BldNodes ! Loop through the blade nodes / elements
+
+
+      ! Integrate to find the portions of the mass matrix on and below the diagonal associated with purely blade DOFs (these portions
+      !   can't be calculated using partial loads).  Also, integrate to find the portions of the forcing vector associated with
+      !    purely blade DOFs (these portions can't be calculated using partial loads):
+      ! NOTE: The vector subscript array, PSBE(), used in the following loops must
+      !       be sorted from smallest to largest DOF index in order for the loops
+      !       to work to enter values only on and below the diagonal of AugMat():
+  
+         DO L = 1,p%DOFs%NPSBE(K)    ! Loop through all active (enabled) blade DOFs that contribute to the QD2T-related linear accelerations of the blade
+            DO I = L,p%DOFs%NPSBE(K) ! Loop through all active (enabled) blade DOFs greater than or equal to L
+               AugMat(p%DOFs%PSBE(K,I),p%DOFs%PSBE(K,L)) = AugMat(p%DOFs%PSBE(K,I),p%DOFs%PSBE(K,L)) + p%BElmntMass(J,K)*&
+                                             DOT_PRODUCT( RtHSdat%PLinVelES(K,J,p%DOFs%PSBE(K,I),0,:), &           ! [C(q,t)]B
+                                                          RtHSdat%PLinVelES(K,J,p%DOFs%PSBE(K,L),0,:)   )
+            ENDDO             ! I - All active (enabled) blade DOFs greater than or equal to L
+         ENDDO                ! L - All active (enabled) blade DOFs that contribute to the QD2T-related linear accelerations of the blade
+      
+         TmpVec1 = RtHSdat%FSAero(:,K,J)*p%DRNodes(J) - p%BElmntMass(J,K)*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccESt(:,K,J) ) ! The portion of FrcS0Bt associated with blade element J
+         TmpVec3 = RtHSdat%MMAero(:,K,J)*p%DRNodes(J)                                               ! The total external moment applied to blade element J
+         DO I = 1,p%DOFs%NPSBE(K)    ! Loop through all active (enabled) blade DOFs that contribute to the QD2T-related linear accelerations of the blade
+               AugMat(p%DOFs%PSBE(K,I), p%NAug) = AugMat(p%DOFs%PSBE(K,I),     p%NAug)                      & ! {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB
+                                           + DOT_PRODUCT( RtHSdat%PLinVelES(K,J,p%DOFs%PSBE(K,I),0,:), TmpVec1 ) & ! NOTE: TmpVec1 is still the portion of FrcS0Bt associated with blade element J
+                                           + DOT_PRODUCT( RtHSdat%PAngVelEM(K,J,p%DOFs%PSBE(K,I),0,:), TmpVec3 )   !       and TmpVec3 is still the total external moment applied to blade element J
+         ENDDO                ! I - All active (enabled) blade DOFs that contribute to the QD2T-related linear accelerations of the blade
+
+
+      ENDDO ! J - Blade nodes / elements      
+      
+      
+      
+
+      ! Initialize the portions of the mass matrix below the diagonal associated
+      !   with the teeter and pure blade DOFs using the partial loads at the teeter pin; only do this if necessary:
+
+      IF ( ( p%NumBl == 2 ) ) THEN
+         IF ( p%DOF_Flag(DOF_Teet) ) THEN ! NOTE: two "ifs" since DOF_Teet might be out of bound
+            DO L = 1,p%DOFs%NPSBE(K) ! Loop through all active (enabled) blade DOFs that contribute to the QD2T-related linear accelerations of the blade
+               AugMat(DOF_Teet,p%DOFs%PSBE(K,L)) = -DOT_PRODUCT( RtHSdat%PAngVelEH(DOF_Teet,0,:), &
+                                                                 RtHSdat%PMomLPRot(:,p%DOFs%PSBE(K,L)) )  ! [C(q,t)]B
+            ENDDO             ! L - All active (enabled) blade DOFs that contribute to the QD2T-related linear accelerations of the blade
+         ENDIF
+      ENDIF
+
+
+
+      ! If the associated DOFs are enabled, add the blade elasticity and damping
+      !   forces to the forcing vector (these portions can't be calculated using
+      !   partial loads):
+
+      IF ( p%DOF_Flag(DOF_BF(K,1)) )  THEN
+         AugMat(    DOF_BF(K,1),p%NAug) = AugMat(DOF_BF(K,1),p%NAug)      & !
+                                        - p%KBF(K,1,1)*x%QT( DOF_BF(K,1)) &
+                                        - p%KBF(K,1,2)*x%QT( DOF_BF(K,2)) &
+                                        - p%CBF(K,1,1)*x%QDT(DOF_BF(K,1)) &
+                                        - p%CBF(K,1,2)*x%QDT(DOF_BF(K,2))
+      ENDIF
+      IF ( p%DOF_Flag(DOF_BF(K,2)) )  THEN
+         AugMat(    DOF_BF(K,2),p%NAug) = AugMat(DOF_BF(K,2),p%NAug)      & ! {-f(qd,q,t)}ElasticB + {-f(qd,q,t)}DampB
+                                        - p%KBF(K,2,1)*x%QT( DOF_BF(K,1)) &
+                                        - p%KBF(K,2,2)*x%QT( DOF_BF(K,2)) &
+                                        - p%CBF(K,2,1)*x%QDT(DOF_BF(K,1)) &
+                                        - p%CBF(K,2,2)*x%QDT(DOF_BF(K,2))
+      ENDIF
+      IF ( p%DOF_Flag(DOF_BE(K,1)) )  THEN
+         AugMat(    DOF_BE(K,1),p%NAug) = AugMat(DOF_BE(K,1),p%NAug)      & !
+                                        - p%KBE(K,1,1)*x%QT( DOF_BE(K,1)) &
+                                        - p%CBE(K,1,1)*x%QDT(DOF_BE(K,1))
+      ENDIF
+                  
+      
+   END DO !k
+         
+   
+      ! Initialize the portions of the mass matrix on and below the diagonal
+      !   associated with purely tower DOFs (these portions can't be calculated
+      !   using partial loads) using the yaw bearing mass effects.
+      !   Also, initialize the portions of the forcing vector associated with
+      !   purely blade DOFs (these portions can't be calculated using partial
+      !   loads) using the yaw bearing mass effects:
+      ! NOTE: The vector subscript array, PTTE(), used in the following loops must
+      !       be sorted from smallest to largest DOF index in order for the loops
+      !       to work to enter values only on and below the diagonal of AugMat():
+
+   DO L = 1,p%DOFs%NPTTE    ! Loop through all active (enabled) tower DOFs that contribute to the QD2T-related linear accelerations of the yaw bearing (point O)
+      DO I = L,p%DOFs%NPTTE ! Loop through all active (enabled) tower DOFs greater than or equal to L
+         AugMat(p%DOFs%PTTE(I),p%DOFs%PTTE(L)) = p%YawBrMass*DOT_PRODUCT( RtHSdat%PLinVelEO(p%DOFs%PTTE(I),0,:), &     ! [C(q,t)]T of YawBrMass
+                                                      RtHSdat%PLinVelEO(p%DOFs%PTTE(L),0,:)    )
+      ENDDO          ! I - All active (enabled) tower DOFs greater than or equal to L
+   ENDDO             ! L - All active (enabled) tower DOFs that contribute to the QD2T-related linear accelerations of the yaw bearing (point O)
+
+   TmpVec1 = -p%YawBrMass*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccEOt ) ! The portion of FrcT0Trbt associated with the YawBrMass
+   DO I = 1,p%DOFs%NPTTE    ! Loop through all active (enabled) tower DOFs that contribute to the QD2T-related linear accelerations of the yaw bearing (point O)
+         AugMat(p%DOFs%PTTE(I),   p%NAug) =           DOT_PRODUCT( RtHSdat%PLinVelEO(p%DOFs%PTTE(I),0,:), &     ! {-f(qd,q,t)}T + {-f(qd,q,t)}GravT of YawBrMass
+                                                      TmpVec1                   )   ! NOTE: TmpVec1 is still the portion of FrcT0Trbt associated with YawBrMass
+   ENDDO             ! I - All active (enabled) tower DOFs that contribute to the QD2T-related linear accelerations of the yaw bearing (point O)
+   
+   
+   
+   
+
+   DO J = 1,p%TwrNodes
+
+   !..................................................................................................................................
+   ! Integrate to find the portions of the mass matrix on and below the diagonal associated with purely tower DOFs (these portions
+   !   can't be calculated using partial loads).  Also, integrate to find the portions of the forcing vector associated with purely
+   !   tower DOFs (these portions can't be calculated using partial loads).
+   ! NOTE: The vector subscript array, PTTE(), used in the following loops must be sorted from smallest to largest DOF index in order
+   !   for the loops to work to enter values only on and below the diagonal of AugMat():
+   !..................................................................................................................................
+
+      DO L = 1,p%DOFs%NPTTE    ! Loop through all active (enabled) tower DOFs that contribute to the QD2T-related linear accelerations of the tower
+         DO I = L,p%DOFs%NPTTE ! Loop through all active (enabled) tower DOFs greater than or equal to L
+            AugMat(p%DOFs%PTTE(I),p%DOFs%PTTE(L)) = AugMat(p%DOFs%PTTE(I),p%DOFs%PTTE(L))  &
+                                                  + p%TElmntMass(J)   *DOT_PRODUCT( RtHSdat%PLinVelET(J,p%DOFs%PTTE(I),0,:),  &
+                                                                              RtHSdat%PLinVelET(J,p%DOFs%PTTE(L),0,:) ) &   ! [C(q,t)]T + [C(q,t)]HydroT
+                                                  - abs(p%DHNodes(J))*DOT_PRODUCT( RtHSdat%PLinVelET(J,p%DOFs%PTTE(I),0,:),  &
+                                                                              RtHSdat%PFTHydro (:,J,p%DOFs%PTTE(L)  ) ) &
+                                                  - abs(p%DHNodes(J))*DOT_PRODUCT( RtHSdat%PAngVelEF(J,p%DOFs%PTTE(I),0,:),  &
+                                                                              RtHSdat%PMFHydro (:,J,p%DOFs%PTTE(L)  ) )
+         ENDDO                 ! I - All active (enabled) tower DOFs greater than or equal to L
+      ENDDO                    ! L - All active (enabled) tower DOFs that contribute to the QD2T-related linear accelerations of the tower
+
+      TmpVec1 = ( RtHSdat%FTHydrot(:,J) )*abs(p%DHNodes(J)) &
+              - p%TElmntMass(J)*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccETt(:,J) )          ! The portion of FrcT0Trbt associated with tower element J
+      TmpVec3 = ( RtHSdat%MFHydrot(:,J) )*abs(p%DHNodes(J))             ! The external moment applied to tower element J
+      DO I = 1,p%DOFs%NPTTE    ! Loop through all active (enabled) tower DOFs that contribute to the QD2T-related linear accelerations of the tower
+            AugMat(p%DOFs%PTTE(I),        p%NAug) = AugMat(p%DOFs%PTTE(I),   p%NAug)                         &                 ! {-f(qd,q,t)}T + {-f(qd,q,t)}GravT + {-f(qd,q,t)}AeroT + {-f(qd,q,t)}HydroT
+                                                  +  DOT_PRODUCT( RtHSdat%PLinVelET(J,p%DOFs%PTTE(I),0,:), TmpVec1        ) &  ! NOTE: TmpVec1 is still the portion of FrcT0Trbt associated with tower element J
+                                                  +  DOT_PRODUCT( RtHSdat%PAngVelEF(J,p%DOFs%PTTE(I),0,:), TmpVec3        )    !       and TmpVec3 is still the total external moment to tower element J
+      ENDDO                    ! I - All active (enabled) tower DOFs that contribute to the QD2T-related linear accelerations of the tower
+
+   ENDDO ! J - Tower nodes / elements   
+   
+   !..................................................................................................................................
+   ! If the associated DOFs are enabled, add the tower elasticity and damping forces to the forcing vector (these portions can't be
+   !   calculated using partial loads):
+   !..................................................................................................................................
+
+   IF ( p%DOF_Flag(DOF_TFA1) )  THEN
+      AugMat(    DOF_TFA1,p%NAug) = AugMat(DOF_TFA1,p%NAug)                                   &
+                                  - p%KTFA(1,1)*x%QT( DOF_TFA1) - p%KTFA(1,2)*x%QT( DOF_TFA2) &                                     !
+                                  - p%CTFA(1,1)*x%QDT(DOF_TFA1) - p%CTFA(1,2)*x%QDT(DOF_TFA2)
+   ENDIF
+   IF ( p%DOF_Flag(DOF_TSS1) )  THEN
+      AugMat(    DOF_TSS1,p%NAug) = AugMat(DOF_TSS1,p%NAug)                                   &
+                                  - p%KTSS(1,1)*x%QT( DOF_TSS1) - p%KTSS(1,2)*x%QT( DOF_TSS2) &                                     ! {-f(qd,q,t)}ElasticT + {-f(qd,q,t)}DampT
+                                  - p%CTSS(1,1)*x%QDT(DOF_TSS1) - p%CTSS(1,2)*x%QDT(DOF_TSS2)
+   ENDIF
+   IF ( p%DOF_Flag(DOF_TFA2) )  THEN
+      AugMat(    DOF_TFA2,p%NAug) = AugMat(DOF_TFA2,p%NAug)                                   &
+                                  - p%KTFA(2,1)*x%QT( DOF_TFA1) - p%KTFA(2,2)*x%QT( DOF_TFA2) &                                     !
+                                  - p%CTFA(2,1)*x%QDT(DOF_TFA1) - p%CTFA(2,2)*x%QDT(DOF_TFA2)
+   ENDIF
+   IF ( p%DOF_Flag(DOF_TSS2) )  THEN
+      AugMat(    DOF_TSS2,p%NAug) = AugMat(DOF_TSS2,p%NAug)                                   &
+                                  - p%KTSS(2,1)*x%QT( DOF_TSS1) - p%KTSS(2,2)*x%QT( DOF_TSS2) &                                     !
+                                  - p%CTSS(2,1)*x%QDT(DOF_TSS1) - p%CTSS(2,2)*x%QDT(DOF_TSS2)
+   ENDIF
+   
+   
+   
+!..................................................................................................................................
+! Now that all of the partial loads have been found, let's fill in the portions of the mass matrix on and below the diagonal that
+! may be calculated with the help of the partial loads.
+! Also, let's fill in the portions of the forcing vector that may be calculated with the help of the partial loads.
+! Also let's add in additional terms to the forcing function that can't be added with the help of the partial loads.
+!
+! NOTE: The vector subscript array, SrtPS(), used in the following loops must be sorted from smallest to largest DOF index in order
+!   for the loops to work to enter values only on and below the diagonal of AugMat():
+!..................................................................................................................................
+
+   IF ( p%DOF_Flag (DOF_Sg  ) )  THEN
+      DO I = p%DOFs%Diag(DOF_Sg  ),p%DOFs%NActvDOF   ! Loop through all active (enabled) DOFs on or below the diagonal
+         AugMat(p%DOFs%SrtPS(I),DOF_Sg  ) = -1.*DOT_PRODUCT( RtHSdat%PLinVelEZ(DOF_Sg ,0,:), RtHSdat%PFrcZAll (:,p%DOFs%SrtPS(I)) ) ! [C(q,t)]X + [C(q,t)]HydroX + [C(q,t)]T + [C(q,t)]HydroT + [C(q,t)]N + [C(q,t)]R + [C(q,t)]H + [C(q,t)]B + [C(q,t)]A
+      ENDDO                            ! I - All active (enabled) DOFs on or below the diagonal
+         AugMat(DOF_Sg  ,         p%NAug) =     DOT_PRODUCT( RtHSdat%PLinVelEZ(DOF_Sg ,0,:), RtHSdat%FrcZAllt              )        ! {-f(qd,q,t)}X + {-f(qd,q,t)}HydroX + {-f(qd,q,t)}T + {-f(qd,q,t)}AeroT + {-f(qd,q,t)}HydroT + {-f(qd,q,t)}N + {-f(qd,q,t)}R + {-f(qd,q,t)}H + {-f(qd,q,t)}B + {-f(qd,q,t)}AeroB + {-f(qd,q,t)}A + {-f(qd,q,t)}AeroA
+   ENDIF
+
+   IF ( p%DOF_Flag (DOF_Sw  ) )  THEN
+      DO I = p%DOFs%Diag(DOF_Sw  ),p%DOFs%NActvDOF   ! Loop through all active (enabled) DOFs on or below the diagonal
+         AugMat(p%DOFs%SrtPS(I),DOF_Sw  ) = -1.*DOT_PRODUCT( RtHSdat%PLinVelEZ(DOF_Sw ,0,:), RtHSdat%PFrcZAll (:,p%DOFs%SrtPS(I)) ) ! [C(q,t)]X + [C(q,t)]HydroX + [C(q,t)]T + [C(q,t)]HydroT + [C(q,t)]N + [C(q,t)]R + [C(q,t)]H + [C(q,t)]B + [C(q,t)]A
+      ENDDO                            ! I - All active (enabled) DOFs on or below the diagonal
+         AugMat(DOF_Sw  ,         p%NAug) =     DOT_PRODUCT( RtHSdat%PLinVelEZ(DOF_Sw ,0,:), RtHSdat%FrcZAllt              )        ! {-f(qd,q,t)}X + {-f(qd,q,t)}HydroX + {-f(qd,q,t)}T + {-f(qd,q,t)}AeroT + {-f(qd,q,t)}HydroT + {-f(qd,q,t)}N + {-f(qd,q,t)}R + {-f(qd,q,t)}H + {-f(qd,q,t)}B + {-f(qd,q,t)}AeroB + {-f(qd,q,t)}A + {-f(qd,q,t)}AeroA
+   ENDIF
+
+   IF ( p%DOF_Flag (DOF_Hv  ) )  THEN
+      DO I = p%DOFs%Diag(DOF_Hv  ),p%DOFs%NActvDOF   ! Loop through all active (enabled) DOFs on or below the diagonal
+         AugMat(p%DOFs%SrtPS(I),DOF_Hv  ) = -1.*DOT_PRODUCT( RtHSdat%PLinVelEZ(DOF_Hv ,0,:), RtHSdat%PFrcZAll (:,p%DOFs%SrtPS(I)) ) ! [C(q,t)]X + [C(q,t)]HydroX + [C(q,t)]T + [C(q,t)]HydroT + [C(q,t)]N + [C(q,t)]R + [C(q,t)]H + [C(q,t)]B + [C(q,t)]A
+      ENDDO                            ! I - All active (enabled) DOFs on or below the diagonal
+         AugMat(DOF_Hv  ,         p%NAug) =     DOT_PRODUCT( RtHSdat%PLinVelEZ(DOF_Hv ,0,:), RtHSdat%FrcZAllt              )        ! {-f(qd,q,t)}X + {-f(qd,q,t)}GravX + {-f(qd,q,t)}HydroX + {-f(qd,q,t)}T + {-f(qd,q,t)}GravT + {-f(qd,q,t)}AeroT + {-f(qd,q,t)}HydroT + {-f(qd,q,t)}N + {-f(qd,q,t)}GravN + {-f(qd,q,t)}R + {-f(qd,q,t)}GravR + {-f(qd,q,t)}H + {-f(qd,q,t)}GravH + {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB + {-f(qd,q,t)}A + {-f(qd,q,t)}GravA + {-f(qd,q,t)}AeroA
+   ENDIF
+
+   IF ( p%DOF_Flag (DOF_R   ) )  THEN
+      DO I = p%DOFs%Diag(DOF_R   ),p%DOFs%NActvDOF   ! Loop through all active (enabled) DOFs on or below the diagonal
+
+         AugMat(p%DOFs%SrtPS(I),DOF_R   ) = -1.*DOT_PRODUCT( RtHSdat%PAngVelEX(DOF_R  ,0,:), RtHSdat%PMomXAll (:,p%DOFs%SrtPS(I)) ) ! [C(q,t)]X + [C(q,t)]HydroX + [C(q,t)]T + [C(q,t)]HydroT + [C(q,t)]N + [C(q,t)]R + [C(q,t)]G + [C(q,t)]H + [C(q,t)]B + [C(q,t)]A
+      ENDDO                            ! I - All active (enabled) DOFs on or below the diagonal
+         AugMat(DOF_R   ,         p%NAug) =     DOT_PRODUCT( RtHSdat%PAngVelEX(DOF_R  ,0,:), RtHSdat%MomXAllt              )        ! {-f(qd,q,t)}X + {-f(qd,q,t)}GravX + {-f(qd,q,t)}HydroX + {-f(qd,q,t)}T + {-f(qd,q,t)}GravT + {-f(qd,q,t)}AeroT + {-f(qd,q,t)}HydroT + {-f(qd,q,t)}N + {-f(qd,q,t)}GravN + {-f(qd,q,t)}R + {-f(qd,q,t)}GravR + {-f(qd,q,t)}G + {-f(qd,q,t)}H + {-f(qd,q,t)}GravH + {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB + {-f(qd,q,t)}A + {-f(qd,q,t)}GravA + {-f(qd,q,t)}AeroA
+   ENDIF
+
+   IF ( p%DOF_Flag (DOF_P   ) )  THEN
+      DO I = p%DOFs%Diag(DOF_P   ),p%DOFs%NActvDOF    ! Loop through all active (enabled) DOFs on or below the diagonal
+         AugMat(p%DOFs%SrtPS(I),DOF_P   ) = -1.*DOT_PRODUCT( RtHSdat%PAngVelEX(DOF_P  ,0,:), RtHSdat%PMomXAll (:,p%DOFs%SrtPS(I)) ) ! [C(q,t)]X + [C(q,t)]HydroX + [C(q,t)]T + [C(q,t)]HydroT + [C(q,t)]N + [C(q,t)]R + [C(q,t)]G + [C(q,t)]H + [C(q,t)]B + [C(q,t)]A
+      ENDDO                                                             ! I - All active (enabled) DOFs on or below the diagonal
+         AugMat(DOF_P            ,p%NAug) =     DOT_PRODUCT( RtHSdat%PAngVelEX(DOF_P  ,0,:), RtHSdat%MomXAllt              )        ! {-f(qd,q,t)}X + {-f(qd,q,t)}GravX + {-f(qd,q,t)}HydroX + {-f(qd,q,t)}T + {-f(qd,q,t)}GravT + {-f(qd,q,t)}AeroT + {-f(qd,q,t)}HydroT + {-f(qd,q,t)}N + {-f(qd,q,t)}GravN + {-f(qd,q,t)}R + {-f(qd,q,t)}GravR + {-f(qd,q,t)}G + {-f(qd,q,t)}H + {-f(qd,q,t)}GravH + {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB + {-f(qd,q,t)}A + {-f(qd,q,t)}GravA + {-f(qd,q,t)}AeroA
+   END IF
+
+   IF ( p%DOF_Flag (DOF_Y   ) )  THEN
+      DO I = p%DOFs%Diag(DOF_Y   ),p%DOFs%NActvDOF    ! Loop through all active (enabled) DOFs on or below the diagonal
+         AugMat(p%DOFs%SrtPS(I),DOF_Y   ) = -1.*DOT_PRODUCT( RtHSdat%PAngVelEX(DOF_Y  ,0,:), RtHSdat%PMomXAll (:,p%DOFs%SrtPS(I)) ) ! [C(q,t)]X + [C(q,t)]HydroX + [C(q,t)]T + [C(q,t)]HydroT + [C(q,t)]N + [C(q,t)]R + [C(q,t)]G + [C(q,t)]H + [C(q,t)]B + [C(q,t)]A
+      ENDDO                                                             ! I - All active (enabled) DOFs on or below the diagonal
+         AugMat(DOF_Y   ,         p%NAug) =     DOT_PRODUCT( RtHSdat%PAngVelEX(DOF_Y  ,0,:), RtHSdat%MomXAllt              )        ! {-f(qd,q,t)}X + {-f(qd,q,t)}GravX + {-f(qd,q,t)}HydroX + {-f(qd,q,t)}T + {-f(qd,q,t)}GravT + {-f(qd,q,t)}AeroT + {-f(qd,q,t)}HydroT + {-f(qd,q,t)}N + {-f(qd,q,t)}GravN + {-f(qd,q,t)}R + {-f(qd,q,t)}GravR + {-f(qd,q,t)}G + {-f(qd,q,t)}H + {-f(qd,q,t)}GravH + {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB + {-f(qd,q,t)}A + {-f(qd,q,t)}GravA + {-f(qd,q,t)}AeroA
+   ENDIF
+
+   IF ( p%DOF_Flag (DOF_TFA1) )  THEN
+      DO I = p%DOFs%Diag(DOF_TFA1),p%DOFs%NActvDOF   ! Loop through all active (enabled) DOFs on or below the diagonal
+         AugMat(p%DOFs%SrtPS(I),DOF_TFA1) = AugMat(p%DOFs%SrtPS(I),DOF_TFA1)                             &
+                                          -  DOT_PRODUCT( RtHSdat%PLinVelEO(DOF_TFA1,0,:),       &
+                                                          RtHSdat%PFrcONcRt(:,p%DOFs%SrtPS(I)) ) &                          ! [C(q,t)]N + [C(q,t)]R + [C(q,t)]G + [C(q,t)]H + [C(q,t)]B + [C(q,t)]A
+                                          -  DOT_PRODUCT( RtHSdat%PAngVelEB(DOF_TFA1,0,:),       &
+                                                          RtHSdat%PMomBNcRt(:,p%DOFs%SrtPS(I)) )
+      ENDDO                            ! I - All active (enabled) DOFs on or below the diagonal
+         AugMat(DOF_TFA1,         p%NAug) = AugMat(DOF_TFA1,    p%NAug)                                  &
+                                          +  DOT_PRODUCT( RtHSdat%PLinVelEO(DOF_TFA1,0,:), RtHSdat%FrcONcRtt  ) &   ! {-f(qd,q,t)}N + {-f(qd,q,t)}GravN + {-f(qd,q,t)}R + {-f(qd,q,t)}GravR + {-f(qd,q,t)}G + {-f(qd,q,t)}H + {-f(qd,q,t)}GravH + {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB + {-f(qd,q,t)}A + {-f(qd,q,t)}GravA + {-f(qd,q,t)}AeroA
+                                          +  DOT_PRODUCT( RtHSdat%PAngVelEB(DOF_TFA1,0,:), RtHSdat%MomBNcRtt  )
+   ENDIF
+
+   IF ( p%DOF_Flag (DOF_TSS1) )  THEN
+      DO I = p%DOFs%Diag(DOF_TSS1),p%DOFs%NActvDOF   ! Loop through all active (enabled) DOFs on or below the diagonal
+         AugMat(p%DOFs%SrtPS(I),DOF_TSS1) = AugMat(p%DOFs%SrtPS(I),DOF_TSS1)                             &
+                                          -  DOT_PRODUCT( RtHSdat%PLinVelEO(DOF_TSS1,0,:),       &
+                                                          RtHSdat%PFrcONcRt(:,p%DOFs%SrtPS(I)) ) &                          ! [C(q,t)]N + [C(q,t)]R + [C(q,t)]G + [C(q,t)]H + [C(q,t)]B + [C(q,t)]A
+                                          -  DOT_PRODUCT( RtHSdat%PAngVelEB(DOF_TSS1,0,:),       &
+                                                          RtHSdat%PMomBNcRt(:,p%DOFs%SrtPS(I)) )
+      ENDDO                            ! I - All active (enabled) DOFs on or below the diagonal
+         AugMat(DOF_TSS1,         p%NAug) = AugMat(DOF_TSS1,    p%NAug)                                  &
+                                          +  DOT_PRODUCT( RtHSdat%PLinVelEO(DOF_TSS1,0,:), RtHSdat%FrcONcRtt  ) &   ! {-f(qd,q,t)}N + {-f(qd,q,t)}GravN + {-f(qd,q,t)}R + {-f(qd,q,t)}GravR + {-f(qd,q,t)}G + {-f(qd,q,t)}H + {-f(qd,q,t)}GravH + {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB + {-f(qd,q,t)}A + {-f(qd,q,t)}GravA + {-f(qd,q,t)}AeroA
+                                          +  DOT_PRODUCT( RtHSdat%PAngVelEB(DOF_TSS1,0,:), RtHSdat%MomBNcRtt  )
+   ENDIF
+
+   IF ( p%DOF_Flag (DOF_TFA2) )  THEN
+      DO I = p%DOFs%Diag(DOF_TFA2),p%DOFs%NActvDOF   ! Loop through all active (enabled) DOFs on or below the diagonal
+         AugMat(p%DOFs%SrtPS(I),DOF_TFA2) = AugMat(p%DOFs%SrtPS(I),DOF_TFA2)                             &
+                                          -  DOT_PRODUCT( RtHSdat%PLinVelEO(DOF_TFA2,0,:),       &
+                                                          RtHSdat%PFrcONcRt(:,p%DOFs%SrtPS(I)) ) &                          ! [C(q,t)]N + [C(q,t)]R + [C(q,t)]G + [C(q,t)]H + [C(q,t)]B + [C(q,t)]A
+                                          -  DOT_PRODUCT( RtHSdat%PAngVelEB(DOF_TFA2,0,:),       &
+                                                          RtHSdat%PMomBNcRt(:,p%DOFs%SrtPS(I)) )
+      ENDDO                            ! I - All active (enabled) DOFs on or below the diagonal
+         AugMat(DOF_TFA2,         p%NAug) = AugMat(DOF_TFA2,    p%NAug)                                  &
+                                          +  DOT_PRODUCT( RtHSdat%PLinVelEO(DOF_TFA2,0,:), RtHSdat%FrcONcRtt  ) &   ! {-f(qd,q,t)}N + {-f(qd,q,t)}GravN + {-f(qd,q,t)}R + {-f(qd,q,t)}GravR + {-f(qd,q,t)}G + {-f(qd,q,t)}H + {-f(qd,q,t)}GravH + {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB + {-f(qd,q,t)}A + {-f(qd,q,t)}GravA + {-f(qd,q,t)}AeroA
+                                          +  DOT_PRODUCT( RtHSdat%PAngVelEB(DOF_TFA2,0,:), RtHSdat%MomBNcRtt  )
+   ENDIF
+
+   IF ( p%DOF_Flag (DOF_TSS2) )  THEN
+      DO I = p%DOFs%Diag(DOF_TSS2),p%DOFs%NActvDOF   ! Loop through all active (enabled) DOFs on or below the diagonal
+         AugMat(p%DOFs%SrtPS(I),DOF_TSS2) = AugMat(p%DOFs%SrtPS(I),DOF_TSS2)                             &
+                                          -  DOT_PRODUCT( RtHSdat%PLinVelEO(DOF_TSS2,0,:),       &
+                                                          RtHSdat%PFrcONcRt(:,p%DOFs%SrtPS(I)) ) &                          ! [C(q,t)]N + [C(q,t)]R + [C(q,t)]G + [C(q,t)]H + [C(q,t)]B + [C(q,t)]A
+                                          -  DOT_PRODUCT( RtHSdat%PAngVelEB(DOF_TSS2,0,:),       &
+                                                          RtHSdat%PMomBNcRt(:,p%DOFs%SrtPS(I)) )
+      ENDDO                            ! I - All active (enabled) DOFs on or below the diagonal
+         AugMat(DOF_TSS2,         p%NAug) = AugMat(DOF_TSS2,    p%NAug)                                  &
+                                          +  DOT_PRODUCT( RtHSdat%PLinVelEO(DOF_TSS2,0,:), RtHSdat%FrcONcRtt  ) &   ! {-f(qd,q,t)}N + {-f(qd,q,t)}GravN + {-f(qd,q,t)}R + {-f(qd,q,t)}GravR + {-f(qd,q,t)}G + {-f(qd,q,t)}H + {-f(qd,q,t)}GravH + {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB + {-f(qd,q,t)}A + {-f(qd,q,t)}GravA + {-f(qd,q,t)}AeroA
+                                          +  DOT_PRODUCT( RtHSdat%PAngVelEB(DOF_TSS2,0,:), RtHSdat%MomBNcRtt  )
+   ENDIF
+   
+   IF ( p%DOF_Flag (DOF_Yaw ) )  THEN
+      DO I = p%DOFs%Diag(DOF_Yaw ),p%DOFs%NActvDOF   ! Loop through all active (enabled) DOFs on or below the diagonal
+         AugMat(p%DOFs%SrtPS(I),DOF_Yaw ) = -DOT_PRODUCT( RtHSdat%PAngVelEN(DOF_Yaw ,0,:), RtHSdat%PMomBNcRt(:,p%DOFs%SrtPS(I)) )   ! [C(q,t)]N + [C(q,t)]R + [C(q,t)]G + [C(q,t)]H + [C(q,t)]B + [C(q,t)]A
+      ENDDO                            ! I - All active (enabled) DOFs on or below the diagonal
+         AugMat(DOF_Yaw ,         p%NAug) =  DOT_PRODUCT( RtHSdat%PAngVelEN(DOF_Yaw ,0,:), RtHSdat%MomBNcRtt             ) &        ! {-f(qd,q,t)}N + {-f(qd,q,t)}GravN + {-f(qd,q,t)}R + {-f(qd,q,t)}GravR + {-f(qd,q,t)}G + {-f(qd,q,t)}H + {-f(qd,q,t)}GravH + {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB + {-f(qd,q,t)}A + {-f(qd,q,t)}GravA + {-f(qd,q,t)}AeroA
+                                                              + u%YawMom                                                            ! + {-f(qd,q,t)}SpringYaw  + {-f(qd,q,t)}DampYaw; NOTE: The neutral yaw rate, YawRateNeut, defaults to zero.  It is only used for yaw control.
+         !PRINT '(A, F30.20)', 'AugMat[11, 25] = ', AugMat(DOF_Yaw ,p%NAug)
+   ENDIF
+   
+   
+   IF ( p%DOF_Flag (DOF_RFrl) )  THEN
+      DO I = p%DOFs%Diag(DOF_RFrl),p%DOFs%NActvDOF   ! Loop through all active (enabled) DOFs on or below the diagonal
+         AugMat(p%DOFs%SrtPS(I),DOF_RFrl) = -DOT_PRODUCT( RtHSdat%PAngVelER(DOF_RFrl,0,:),       &
+                                                          RtHSdat%PMomNGnRt(:,p%DOFs%SrtPS(I)) )                            ! [C(q,t)]R + [C(q,t)]G + [C(q,t)]H + [C(q,t)]B
+      ENDDO                            ! I - All active (enabled) DOFs on or below the diagonal
+         AugMat(DOF_RFrl,         p%NAug) =  DOT_PRODUCT( RtHSdat%PAngVelER(DOF_RFrl,0,:), RtHSdat%MomNGnRtt  ) &   ! {-f(qd,q,t)}R + {-f(qd,q,t)}GravR + {-f(qd,q,t)}G + {-f(qd,q,t)}H + {-f(qd,q,t)}GravH + {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB
+                                                                              +  RtHSdat%RFrlMom                            ! + {-f(qd,q,t)}SpringRF + {-f(qd,q,t)}DampRF
+   ENDIF
+
+   TmpVec = p%GenIner*CoordSys%c1*DOT_PRODUCT( CoordSys%c1, RtHSdat%PAngVelEG(DOF_GeAz,0,:) )  ! = ( generator inertia dyadic ) Dot ( partial angular velocity of G in E for DOF_GeAz )
+
+   IF ( p%DOF_Flag (DOF_GeAz) )  THEN
+      DO I = p%DOFs%Diag(DOF_GeAz),p%DOFs%NActvDOF   ! Loop through all active (enabled) DOFs on or below the diagonal
+         AugMat(p%DOFs%SrtPS(I),DOF_GeAz) = -DOT_PRODUCT( RtHSdat%PAngVelEL(DOF_GeAz,0,:), RtHSdat%PMomLPRot(:,p%DOFs%SrtPS(I)) )! [C(q,t)]H + [C(q,t)]B
+      ENDDO                            ! I - All active (enabled) DOFs on or below the diagonal
+         AugMat(DOF_GeAz,         p%NAug) =  DOT_PRODUCT( RtHSdat%PAngVelEL(DOF_GeAz,0,:), RtHSdat%MomLPRott             ) &     ! {-f(qd,q,t)}H + {-f(qd,q,t)}GravH + {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB
+                                                              -  GBoxTrq                                                         ! + {-f(qd,q,t)}Gen + {-f(qd,q,t)}Brake
+
+
+      ! The previous loop (DO I = p%DOFs%Diag(DOF_GeAz),p%DOFs%NActvDOF) misses the
+      !   generator inertia-contribution to the mass matrix and forcing function.
+      !   Thus, add these in as well:
+
+
+         AugMat(DOF_GeAz,       DOF_GeAz) = AugMat(DOF_GeAz,DOF_GeAz)                                    &
+                                            +  DOT_PRODUCT( RtHSdat%PAngVelEG(DOF_GeAz,0,:), TmpVec                )             ! [C(q,t)]G
+         AugMat(DOF_GeAz,         p%NAug) = AugMat(DOF_GeAz,  p%NAug)                                    &
+                                            -  DOT_PRODUCT( RtHSdat%AngAccEGt              , TmpVec                )             ! {-f(qd,q,t)}G
+
+
+   ENDIF
+
+   IF ( p%DOF_Flag (DOF_DrTr) )  THEN
+      DO I = p%DOFs%Diag(DOF_DrTr),p%DOFs%NActvDOF   ! Loop through all active (enabled) DOFs on or below the diagonal
+         AugMat(p%DOFs%SrtPS(I),DOF_DrTr) = -DOT_PRODUCT( RtHSdat%PAngVelEL(DOF_DrTr,0,:), RtHSdat%PMomLPRot(:,p%DOFs%SrtPS(I)) ) ! [C(q,t)]H + [C(q,t)]B
+      ENDDO                            ! I - All active (enabled) DOFs on or below the diagonal
+         AugMat(DOF_DrTr,         p%NAug) =  DOT_PRODUCT( RtHSdat%PAngVelEL(DOF_DrTr,0,:), RtHSdat%MomLPRott             ) &      ! {-f(qd,q,t)}H + {-f(qd,q,t)}GravH + {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB
+                                                          -  p%DTTorSpr*x%QT (DOF_DrTr)                                    &      ! + {-f(qd,q,t)}ElasticDrive
+                                                          -  p%DTTorDmp*x%QDT(DOF_DrTr)                                           ! + {-f(qd,q,t)}DampDrive
+   ENDIF
+
+   IF ( p%DOF_Flag (DOF_TFrl) )  THEN
+      ! The tail-furl DOF does not affect any DOF index larger than DOF_TFrl.  Therefore, there is no need to perform the loop: DO I = Diag(DOF_TFrl),NActvDOF
+         AugMat(DOF_TFrl,       DOF_TFrl) = -DOT_PRODUCT( RtHSdat%PAngVelEA(DOF_TFrl,0,:), RtHSdat%PMomNTail(:,DOF_TFrl) )        ! [C(q,t)]A
+         AugMat(DOF_TFrl,         p%NAug) =  DOT_PRODUCT( RtHSdat%PAngVelEA(DOF_TFrl,0,:), RtHSdat%MomNTailt             ) &      ! {-f(qd,q,t)}A + {-f(qd,q,t)}GravA + {-f(qd,q,t)}AeroA
+                                                              +  RtHSdat%TFrlMom                                                  ! + {-f(qd,q,t)}SpringTF + {-f(qd,q,t)}DampTF
+   ENDIF
+
+   IF ( ( p%NumBl == 2 ) ) THEN 
+      IF ( p%DOF_Flag(DOF_Teet) )  THEN  ! NOTE: two "ifs" since DOF_Teet may be out of bound
+         ! The teeter DOF does not affect any DOF index larger than DOF_Teet.  Therefore, there is no need to perform the loop: DO I = Diag(DOF_Teet),NActvDOF
+         AugMat(DOF_Teet,       DOF_Teet) = -DOT_PRODUCT( RtHSdat%PAngVelEH(DOF_Teet,0,:), RtHSdat%PMomLPRot(:,DOF_Teet) )        ! [C(q,t)]H + [C(q,t)]B
+         AugMat(DOF_Teet,         p%NAug) =  DOT_PRODUCT( RtHSdat%PAngVelEH(DOF_Teet,0,:), RtHSdat%MomLPRott             ) &      ! {-f(qd,q,t)}H + {-f(qd,q,t)}GravH + {-f(qd,q,t)}B + {-f(qd,q,t)}GravB + {-f(qd,q,t)}AeroB
+                                                              +  RtHSdat%TeetMom                                                  ! + {-f(qd,q,t)}SpringTeet + {-f(qd,q,t)}DampTeet
+      ENDIF
+   ENDIF
+   !..................................................................................................................................
+   ! So far, we have only filled in the portions of the mass matrix on and below the diagonal.  Because the mass matrix is symmetric
+   !   up to this point, let's fill in the portion above the diagonal by mirroring the values from below:
+   ! NOTE: The vector subscript array, SrtPS(), used in the following loops must be sorted from smallest to largest DOF index in order
+   !   for the loops to work to enter values only on and below the diagonal of AugMat():
+   !..................................................................................................................................
+
+      DO L = 2,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs above the diagonal (columns)
+         DO I = 1,L-1   ! Loop through all active (enabled) DOFs above the diagonal (rows)
+            AugMat(p%DOFs%SrtPS(I),p%DOFs%SrtPS(L)) = AugMat(p%DOFs%SrtPS(L),p%DOFs%SrtPS(I))
+         ENDDO          ! I - All active (enabled) DOFs above the diagonal (rows)
+      ENDDO             ! L - All active (enabled) DOFs above the diagonal (columns)
+
+   ! Let's add the gearbox friction terms to the mass matrix and forcing
+   !   function.  These only effect the equation for the generator azimuth DOF.
+   ! NOTE: the MASS MATRIX WILL NO LONGER BE SYMMETRIC after adding these
+   !       terms, unless the gearbox efficiency, GBoxEff, was set to 100%:
+   
+   
+   GBoxEffFac2 = ( 1.0/RtHSdat%GBoxEffFac - 1.0 ) ! = ( 1 / GBoxEff^SgnPrvLSTQ - 1 )
+   !TmpVec = p%GenIner*CoordSys%c1*DOT_PRODUCT( CoordSys%c1, RtHSdat%PAngVelEG(DOF_GeAz,0,:) )  ! = ( generator inertia dyadic ) Dot ( partial angular velocity of G in E for DOF_GeAz )
+
+   DO I = 1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs
+
+      AugMat(DOF_GeAz,p%DOFs%SrtPS(I)) = AugMat(DOF_GeAz,p%DOFs%SrtPS(I)) &                                          ! NOTE: TmpVec is still = ( generator inertia dyadic ) Dot ( partial angular velocity of G in E for DOF_GeAz ) in the following equation
+                                + GBoxEffFac2*  DOT_PRODUCT( RtHSdat%PAngVelEG(p%DOFs%SrtPS(I),0,:), TmpVec )        ! [C(q,t)]GBFric
+
+   ENDDO             ! I - All active (enabled) DOFs
+
+   AugMat(   DOF_GeAz,    p%NAug) = AugMat(DOF_GeAz,    p%NAug) &                                                    ! NOTE: TmpVec is still = ( generator inertia dyadic ) Dot ( partial angular velocity of G in E for DOF_GeAz ) in the following equation
+                                - GBoxEffFac2*( DOT_PRODUCT( RtHSdat%AngAccEGt              , TmpVec ) + GBoxTrq )   ! {-f(qd,q,t)}GBFric
+
+   
+   
+   
+END SUBROUTINE FillAugMat
+```

书籍/力学书籍/OpenFast/FASTKinetics/auto/Tower.md

@@ -0,0 +1,155 @@
+
+The distributed properties of the tower bring about generalized inertia forces and generalized active forces associated with tower elasticity, tower damping, tower aerodynamics, and tower weight.  Note that I eliminated the tower mass tuners, since it is redundant to have both mass and stiffness tuners when trying to tune tower frequencies (to tune the frequencies for individual modes, all that is needed is to tune the mass or the stiffness for the individual modes, but not both).  Note also that I eliminated the tower stiffness tuner’s effects on the gravitational destiffening loads.  It is also beneficial to eliminate the tower mass tuners because the tower mass density is needed to compute the tower base loads and thus these tuners affect the tower base loads directly—this makes the form of the tower base load equations considerably more complex and considerably less intuitive. Since the tower elastic stiffness does not directly influence the tower base loads in a fundamental way, the retention of the tower stiffness tuners is much more favorable than the retention of the tower mass tuners (recall that only one set of tuners needs to be retained in order to permit the user to match natural frequencies).  The elimination of the tower stiffness tuner’s effects on the gravitational destiffening was done for the same reason (i.e., the gravity loads directly affect the tower base loads, and thus, tower stiffness tuners make the form of the tower base load equations considerable more complex and considerably less intuitive).  The fact that the gravitational destiffening of the tower is small compared to the overall stiffness of the tower is another reason this elimination of stiffness tuning effects should not be of significant concern.  
+
+塔架的分布式特性带来了与塔架弹性、塔架阻尼、塔架气动特性和塔架重量相关的广义惯性力和广义主动力。值得注意的是，我消除了塔架质量调谐器，因为在尝试调谐塔架频率时，同时拥有质量和刚度调谐器是多余的（为了调谐单个模态的频率，所需的只是调谐单个模态的质量或刚度，而不是两者都调谐）。另请注意，我消除了塔架刚度调谐器对重力去刚化载荷的影响。消除塔架质量调谐器也是有益的，因为计算塔基载荷需要塔架质量密度，因此这些调谐器直接影响塔基载荷——这使得塔基载荷方程的形式变得相当复杂且直观性大大降低。由于塔架弹性刚度从根本上不直接影响塔基载荷，因此保留塔架刚度调谐器比保留塔架质量调谐器更有利（回想一下，为了允许用户匹配固有频率，只需保留一套调谐器）。消除塔架刚度调谐器对重力去刚化的影响也是出于同样的原因（即，重力载荷直接影响塔基载荷，因此，塔架刚度调谐器使得塔基载荷方程的形式变得相当复杂且直观性大大降低）。塔架的重力去刚化相对于塔架的整体刚度而言很小，这是这种消除刚度调谐影响的做法不应引起重大关注的另一个原因。
+
+$$
+F_r^* \Big|_T = - \int_0^{TwFlexL} \mu^T (h) \cdot {}^E v_r^T (h) \cdot {}^E a^T (h)dh - YawBrMass {}^E v_r^O \cdot {}^E a^O \quad (r = 1, 2, \ldots, 22) \\
+\text{where} \quad \mu^T (h) = AdjTwMa \cdot TMassDen(h)\tag{1}
+$$
+
+or, 
+
+$$
+F_r^* \Big|_T = - \int_0^{TwrFlexL} \mu^T (h) \cdot {}^E \mathbf{v}_r^T (h) \cdot \left\{ \left( \sum_{i=1}^{10} {}^E \mathbf{v}_i^T (h) \ddot{q}_i \right) + \left[ \sum_{i=4}^{10} \frac{d}{dt} \left( {}^E \mathbf{v}_i^T (h) \right) \dot{q}_i \right] \right\} dh - YawBrMass {}^E \mathbf{v}_r^O \cdot \left\{ \left( \sum_{i=1}^{10} {}^E \mathbf{v}_i^O \ddot{q}_i \right) + \left[ \sum_{i=4}^{10} \frac{d}{dt} \left( {}^E \mathbf{v}_i^O \right) \dot{q}_i \right] \right\} \quad (r=1,2,\ldots,10)\tag{2}
+$$
+
+Thus, 
+$$
+[C(\mathbf{q},t)]_T (Row, Col) = \int_0^{TwrFlexL} \mu^T (h) {}^E \mathbf{v}_{Row}^T (h) \cdot {}^E \mathbf{v}_{Col}^T (h)dh - YawBrMass {}^E \mathbf{v}_{Row}^O \cdot {}^E \mathbf{v}_{Col}^O \quad (Row, Col = 1, 2, \ldots, 10)\tag{3}
+$$
+$$
+\left\{ - f(\dot{q},q,t) \right\} \Big|_T (Row) = - \int_0^{TwrFlexL} \mu^T (h) {}^E \mathbf{v}_{Row}^T (h) \cdot \left( \sum_{i=4}^{10} \frac{d}{dt} \left[ {}^E \mathbf{v}_i^T (h) \right] \dot{q}_i \right) dh - YawBrMass {}^E \mathbf{v}_{Row}^O \cdot \left[ \sum_{i=4}^{10} \frac{d}{dt} \left( {}^E \mathbf{v}_i^O \right) \dot{q}_i \right] \quad (Row = 1,2,\ldots,10)\tag{4}
+$$
+
+$$
+F_{r \Big| ElasticT} = - \frac{\partial V'^T}{\partial q_r} \quad (r = 1,2,\ldots,22)\tag{5}
+$$
+
+so, 
+$$
+F_{r \Big| ElasticT} = \left\{
+\begin{array}{ll}
+-k_{11}^{'TFA} q_{TFA1} - k_{12}^{'TFA} q_{TFA2} & \text{for } r=TFA1 \\
+-k_{11}^{'TSS} q_{TSS1} - k_{12}^{'TSS} q_{TSS2} & \text{for } r=TSS1 \\
+-k_{21}^{'TFA} q_{TFA1} - k_{22}^{'TFA} q_{TFA2} & \text{for } r=TFA2 \\
+-k_{21}^{'TSS} q_{TSS1} - k_{22}^{'TSS} q_{TSS2} & \text{for } r=TSS2 \\
+0 & \text{otherwise}
+\end{array}
+\right.\tag{6}
+$$
+
+where $k_{\ i j}^{\,\prime T F A}$ and $k_{{i j}}^{\,\prime{T S S}}$ are the generalized stiffnesses of the tower in the fore-aft and side-to-side directions respectively when gravitational destiffening effects are not included as follows: 
+其中，$k_{\ i j}^{\,\prime T F A}$ 和 $k_{{i j}}^{\,\prime{T S S}}$ 分别是塔筒在前后和侧向方向上的广义刚度，此时不包括重力去刚化效应，如下所示： 
+$$
+k_{ij}^{'TFA} = \sqrt{FAStTunr(i) FAStTunr(j)} \int_0^{TwrFlexL} EI^{TFA}(h) \frac{d^2 \phi_i^{TFA}(h)}{dh^2} \frac{d^2 \phi_j^{TFA}(h)}{dh^2} dh \quad (i, j = 1,2) \text{ (which is symmetric)}\tag{7}
+$$
+where 
+$$
+EI^{TFA}(h) = AdjFASt \cdot TwFAStif(h)\tag{8}
+$$
+
+$$
+k_{ij}^{'TSS} = \sqrt{SSStTunr(i) SSStTunr(j)} \int_0^{TwrFlexL} EI^{TSS}(h) \frac{d^2 \phi_i^{TSS}(h)}{dh^2} \frac{d^2 \phi_j^{TSS}(h)}{dh^2} dh \quad (i, j = 1,2) \text{ (which is symmetric)}\tag{9}
+$$
+
+where 
+$$
+EI^{TSS}(h) = AdjSSSt \cdot TwSSStif(h)\tag{10}
+$$
+
+The coefficient in front of the integral in these generalized stiffnesses represents the individual modal stiffness tuning, which allows the user to vary the stiffness of the tower between the individual modes to permit better matching of tower frequencies.  To be precise, the tuner coefficient only really makes sense when working with a generalized stiffness of a single mode (i.e., $k_{\;\;I I}^{\;\prime T F A},\;k_{\;\;22}^{\;\prime T F A},\;k_{\;\;I I}^{\;\prime T S S}$ , or $k\,_{\,22}^{\prime T S S}$ ), in which case the coefficient for mode $i$ is simply $F A S t T u n r(i)$ or SSStTunr $(i)$ .  However, since the cross-correlation elements of the generalized stiffness matrix will, in general, not vanish, the coefficient in the form above permits the tuning to apply to these terms in a consistent fashion. 
+这些广义刚度中的积分系数代表了各个模态刚度调谐，允许用户在各个模态之间改变塔的刚度，从而更好地匹配塔的固有频率。 值得注意的是，调谐系数的意义主要体现在针对单一模态的广义刚度时（例如，$k_{\;\;I I}^{\;\prime T F A},\;k_{\;\;22}^{\;\prime T F A},\;k_{\;\;I I}^{\;\prime T S S}$ , 或 $k\,_{\,22}^{\prime T S S}$ ），此时第 $i$ 阶模态的系数仅仅是 $F A S t T u n r(i)$ 或 SSStTunr $(i)$ 。 然而，由于广义刚度矩阵的互相关元素通常不会消失，因此上述形式的系数允许以一致的方式将调谐应用于这些项。 
+
+Similarly, when using the Rayleigh damping technique where the damping is assumed proportional to the stiffness, then
+同样地，当使用瑞利阻尼技术（其中阻尼被假定与刚度成比例）时，则
+$$
+F_r |_{\text{Damp}T} = \begin{cases}
+- \frac{\zeta_1^{TFA} k_{11}^{'TFA}}{\pi f_1^{'TFA}} \dot{q}_{TFA1} - \frac{\zeta_2^{TFA} k_{12}^{'TFA}}{\pi f_2^{'TFA}} \dot{q}_{TFA2} & \text{for } r=TFA1 \\
+- \frac{\zeta_1^{TSS} k_{11}^{'TSS}}{\pi f_1^{'TSS}} \dot{q}_{TSS1} - \frac{\zeta_2^{TSS} k_{12}^{'TSS}}{\pi f_2^{'TSS}} \dot{q}_{TSS2} & \text{for } r=TSS1 \\
+- \frac{\zeta_1^{TFA} k_{21}^{'TFA}}{\pi f_1^{'TFA}} \dot{q}_{TFA1} - \frac{\zeta_2^{TFA} k_{22}^{'TFA}}{\pi f_2^{'TFA}} \dot{q}_{TFA2} & \text{for } r=TFA2 \\
+- \frac{\zeta_1^{TSS} k_{21}^{'TSS}}{\pi f_1^{'TSS}} \dot{q}_{TSS1} - \frac{\zeta_2^{TSS} k_{22}^{'TSS}}{\pi f_2^{'TSS}} \dot{q}_{TSS2} & \text{for } r=TSS2 \\
+0 & \text{otherwise}
+\end{cases}\tag{11}
+$$
+
+where $\zeta_{i}^{T F A}$ and $\zeta_{i}^{T S S}$ represent  the  structural  damping  ratio  of  the  tower  for  the $i^{\mathrm{th}}$ mode  in  the  fore-aft  and  side-to-side  directions, $T w r F A D m p(i)/l O O$ and $T w r S S D m p(i)/l O O$ respectively.  Also, $f_{\ i}^{\,\prime T F A}$ and $f_{\ i}^{\prime T S S}$ represent the natural frequency of the tower for the $i^{\mathrm{th}}$ mode in the fore-aft and side-to-side directions respectively without tower-top mass or gravitational destiffening effects.  That is,  
+其中，$\zeta_{i}^{T F A}$ 和 $\zeta_{i}^{T S S}$ 分别表示塔针对第 $i$ 阶模态在前后向和侧向上的结构阻尼比，分别记为 $T w r F A D m p(i)/100$ 和 $T w r S S D m p(i)/100$。 此外，$f_{\ i}^{\,\prime T F A}$ 和 $f_{\ i}^{\prime T S S}$ 分别表示塔针对第 $i$ 阶模态在前后向和侧向上的固有频率，不考虑塔顶质量和重力软化效应。 即，
+
+$$
+f_{\;\;i}^{\;\prime^{T F A}}=\frac{1}{2\pi}\sqrt{\frac{k_{\;i i}^{\prime^{T F A}}}{m_{\;i i}^{\prime^{T F A}}}}\;\;\;\;\;\mathrm{and}\;\;\;\;\;\;\;f_{\;i}^{\;\prime^{T S S}}=\frac{1}{2\pi}\sqrt{\frac{k_{\;i i}^{\prime T S S}}{m_{\;i i}^{\prime T S S}}}\tag{12}
+$$  
+
+where $m_{\;i i}^{\;\prime T F A}$ and $m_{\ i i}^{\prime T S S}$ represent the generalized mass of the tower for the $i^{\mathrm{th}}$ mode in the fore-aft and side-to-side directions respectively without tower-top mass effects as follows:  
+其中，$m_{\;i i}^{\;\prime T F A}$ 和 $m_{\ i i}^{\prime T S S}$ 分别表示在塔顶质量效应忽略的情况下，第 $i$ 阶模态在前后方向和左右方向上的广义质量，具体如下：
+
+$$
+m_{\ i j}^{\prime^{T F A}}=\int_{0}^{T w r F l e x L}\mu^{T}\left(h\right)\phi_{i}^{T F A}\left(h\right)\phi_{j}^{T F A}\left(h\right)d h\quad\left(i,j=1,2\right)\tag{13}
+$$  
+
+and  
+
+$$
+m_{\it i j}^{\prime T S S}=\int_{0}^{T w r F l e x L}\mu^{T}\left(h\right)\phi_{i}^{T S S}\left(h\right)\phi_{j}^{T S S}\left(h\right)d h\quad\left(i,j=1,2\right)\tag{14}
+$$  
+Thus, 
+$$
+[C(q,t)]|_{\text{Elastic}T} = 0\tag{15}
+$$
+
+$$
+\{-f(\dot{q},q,t)\}|_{\text{Elastic}T} = \begin{Bmatrix}
+\vdots \\
+\vdots \\
+\vdots \\
+-k_{11}^{'TFA} q_{TFA1} - k_{12}^{'TFA} q_{TFA2} \\
+-k_{11}^{'TSS} q_{TSS1} - k_{12}^{'TSS} q_{TSS2} \\
+-k_{21}^{'TFA} q_{TFA1} - k_{22}^{'TFA} q_{TFA2} \\
+-k_{21}^{'TSS} q_{TSS1} - k_{22}^{'TSS} q_{TSS2} \\
+\vdots \\
+\vdots \\
+\vdots \\
+\end{Bmatrix}\tag{16}
+$$
+
+$$
+[C(q,t)]|_{\text{Damp}T} = 0\tag{17}
+$$
+
+$$
+\{-f(\dot{q},q,t)\}|_{\text{Damp}T} = \begin{Bmatrix}
+\vdots \\
+\vdots \\
+\vdots \\
+- \frac{\zeta_1^{TFA} k_{11}^{'TFA}}{\pi f_1^{'TFA}} \dot{q}_{TFA1} - \frac{\zeta_2^{TFA} k_{12}^{'TFA}}{\pi f_2^{'TFA}} \dot{q}_{TFA2} \\
+- \frac{\zeta_1^{TSS} k_{11}^{'TSS}}{\pi f_1^{'TSS}} \dot{q}_{TSS1} - \frac{\zeta_2^{TSS} k_{12}^{'TSS}}{\pi f_2^{'TSS}} \dot{q}_{TSS2} \\
+- \frac{\zeta_1^{TFA} k_{21}^{'TFA}}{\pi f_1^{'TFA}} \dot{q}_{TFA1} - \frac{\zeta_2^{TFA} k_{22}^{'TFA}}{\pi f_2^{'TFA}} \dot{q}_{TFA2} \\
+- \frac{\zeta_1^{TSS} k_{21}^{'TSS}}{\pi f_1^{'TSS}} \dot{q}_{TSS1} - \frac{\zeta_2^{TSS} k_{22}^{'TSS}}{\pi f_2^{'TSS}} \dot{q}_{TSS2} \\
+\vdots \\
+\vdots \\
+\vdots
+\end{Bmatrix}\tag{18}
+$$
+$$
+F_r |_{\text{Grav}T} = \int_{0}^{\text{TwrFlexL}} \! {}^E \mathbf{v}_r^T(h) \cdot [-\mu^T(h) \mathbf{g}\mathbf{z}_2] \, dh + {}^E \mathbf{v}_r^O \cdot (-\text{YawBrMass} \cdot \mathbf{g}\mathbf{z}_2) \quad (r=3,4,...,10)\tag{19}
+$$
+Thus, 
+$$
+[C(q,t)]|_{\text{Grav}T} = 0\tag{20}
+$$
+
+$$
+\{-f(\dot{q},q,t)\}|_{\text{Grav}T} (\text{Row}) = - \int_{0}^{\text{TwrFlexL}} \! \mu^T(h) g {}^E \mathbf{v}_{\text{Row}}^T(h) \cdot \mathbf{z}_2 \, dh - \text{YawBrMass} \cdot g {}^E \mathbf{v}_{\text{Row}}^O \cdot \mathbf{z}_2 \quad (Row=3,4,...,10)\tag{21}
+$$
+$$
+F_r |_{\text{Aero}T} = \int_{0}^{\text{TwrFlexL}} \! \left[ {}^E \mathbf{v}_r^T(h) \cdot \mathbf{F}_{\text{Aero}T}^T(h) + {}^E \boldsymbol{\omega}_r^F(h) \cdot \mathbf{M}_{\text{Aero}T}^F(h) \right] dh \quad (r=1,2,...,10)\tag{22}
+$$
+
+where $F_{A e r o T}^{T}\left(h\right)$ and $M_{A e r o T}^{F}\left(h\right)$ are aerodynamic forces and moments applied to point $\mathrm{T}$ on the tower respectively expressed per unit height. Thus,  
+$$
+[C(q,t)]|_{\text{Aero}T} = 0\tag{23}
+$$
+
+$$
+\{-f(\dot{q},q,t)\}|_{\text{Aero}T} (\text{Row}) = \int_{0}^{\text{TwrFlexL}} \! \left[ {}^E \mathbf{v}_{\text{Row}}^T(h) \cdot \mathbf{F}_{\text{Aero}T}^T(h) + {}^E \boldsymbol{\omega}_{\text{Row}}^F(h) \cdot \mathbf{M}_{\text{Aero}T}^F(h) \right] dh \quad (Row=1,2,...,10)\tag{24}
+$$

工作OKRs/25.9-11 OKR.canvas

@@ -6,7 +6,7 @@
 		{"id":"82708a439812fdc7","type":"text","text":"# 10月已完成\n\n- 编写输出量\n\t- 方法测试 done\n\t- 批量输出编写 done\n\nP2 湍流 气动 多体 控制联调 \n- 更换600s湍流风 done\n\n- 低风速震荡问题解决 不发电？bug解决\n- 测试额定风速以下，附近，额定风速上稳态风及湍流风工况\n\nP1 线性化方法编写 ing\n- 再把握理论 开启哪些自由度 done\n\n","x":-220,"y":134,"width":440,"height":560},
 		{"id":"505acb3e6b119076","type":"text","text":"# 9月已完成\n\nP1 湍流 气动 多体 控制联调  done\n- 5mw 通了\n\t- 纯叶片变形\n\t- 纯塔架变形\n\t- 叶片+塔架变形 \n\nP1 bladed对比--稳态运行载荷，产出报告\n- 气动新版本稳态跑通 done\n- 多体模块输出方法更新\n\t- 增加输出结构体\n\t- 新增输出宏，自动编写变量头和输出\n\t- 输出位置从ed_calcoutput拆分\n- 稳态变形量对比\n- 所有输出量梳理，对比","x":-700,"y":134,"width":440,"height":560},
 		{"id":"30cb7486dc4e224c","type":"text","text":"# 11月已完成\n\n\n\n","x":260,"y":134,"width":440,"height":560},
-		{"id":"c18d25521d773705","type":"text","text":"# 计划\n这周要做的3~5件重要的事情，这些事情能有效推进实现OKR。\n\nP1 必须做。P2 应该做\n\n\nP2 柔性部件 叶片、塔架变形算法 主线\n- 变形体动力学 简略看看ing\n- 柔性梁弯曲变形振动学习，主线 \n\t- 广义质量 刚度矩阵及含义\n\t\n- 梳理bladed动力学框架\n\t- 子结构文献阅读\n\t- 叶片模型建模 done\n- 共旋方法学习\n- DTU 变形量计算方法学习\n\n\nP1 线性化方法编写 ing\n\n- 开始编写扰动代码\n- 形成系统矩阵\n- 分析矩阵\n\n\nP2 气动、多体、控制、水动联调\n\nP1 启机工况等调试\n\nP1 稳态工况前端对接-董\n- 是否拆分成单独的bin，等待气动完成后开始\n- 如何接收参数  配置文件 \n\nP1 专利\n- 做出solidworks模型，写专利\n\nP2 bladed对比--稳态运行载荷，产出报告\n- 气动参与模块对比\n- 模态对比   两种描述方法不同，bladed方向更多，x y z deflection, x y z rotation，不好对比\n- 气动对比 aerodynamic info 轴向切向诱导因子，根部，尖部差距较大\n\nP2 yaw 自由度再bug确认 已知原理了\n","x":-597,"y":-803,"width":453,"height":457},
+		{"id":"c18d25521d773705","type":"text","text":"# 计划\n这周要做的3~5件重要的事情，这些事情能有效推进实现OKR。\n\nP1 必须做。P2 应该做\n\n\nP2 柔性部件 叶片、塔架变形算法 主线\n- 变形体动力学 简略看看ing\n- 柔性梁弯曲变形振动学习，主线 \n\t- 广义质量 刚度矩阵及含义\n\t\n- 梳理bladed动力学框架\n\t- 子结构文献阅读\n\t- 叶片模型建模 done\n- 共旋方法学习\n- DTU 变形量计算方法学习\n\n\nP1 线性化方法编写 ing\n\n- 其他自由度稳态计算方法公式推导和算法编写\n\t- drtr\n\t- tower\n- 开始编写扰动代码\n- 形成系统矩阵\n- 分析矩阵\n\n\nP2 气动、多体、控制、水动联调\n\nP1 启机工况等调试\n- 兰总进度\n\nP1 稳态工况前端对接-董\n- 是否拆分成单独的bin，等待气动完成后开始\n- 如何接收参数  配置文件 \n\nP1 专利\n- 做出solidworks模型，写专利\n\nP2 bladed对比--稳态运行载荷，产出报告\n- 气动参与模块对比\n- 模态对比   两种描述方法不同，bladed方向更多，x y z deflection, x y z rotation，不好对比\n- 气动对比 aerodynamic info 轴向切向诱导因子，根部，尖部差距较大\n\nP2 yaw 自由度再bug确认 已知原理了\n","x":-597,"y":-803,"width":453,"height":457},
 		{"id":"86ab96a25a3bf82e","type":"text","text":" 湍流风+ 控制的联调，bladed也算一个算例\n- 加水动的联调\n- 8月份底完成这两个\n- 9月份完成停机等工况测试\n- 10月份明阳实际机型测试","x":580,"y":-803,"width":480,"height":220},
 		{"id":"e355f33c92cf18ea","type":"text","text":"9月份定常计算对接前端\n非定常测试完也对接前端","x":580,"y":-500,"width":480,"height":100},
 		{"id":"859e6853b7f1b92b","type":"text","text":"年底考核：\n专利\n线性化模块","x":1200,"y":-803,"width":320,"height":110}