From 7d12a4abc045f9b4b03a56fc17fd949ddb071cb0 Mon Sep 17 00:00:00 2001 From: yize Date: Wed, 7 Jan 2026 17:42:48 +0800 Subject: [PATCH] vault backup: 2026-01-07 17:42:46 --- .obsidian/plugins/copilot/data.json | 2 +- ...7_103211__[[Blade]]详细展示公式推导过程.md | 799 ++++++++++++++++++ .../OpenFast/FASTKinetics/auto/Blade.md | 94 ++- .../FASTKinetics/auto/CalculateAngularPosVelPAcc.md | 222 +++++ .../FASTKinetics/auto/CalculateForcesMoments.md | 634 ++++++++++++++ .../FASTKinetics/auto/CalculatePositions.md | 143 ++++ .../OpenFast/FASTKinetics/auto/Coeff.md | 678 +++++++++++++++ .../明阳机型转速发散 叶片变形量为负.md | 14 +- 未命名 1.base | 3 + 9 files changed, 2550 insertions(+), 39 deletions(-) create mode 100644 copilot-conversations/20260107_103211__[[Blade]]详细展示公式推导过程.md create mode 100644 书籍/力学书籍/OpenFast/FASTKinetics/auto/CalculateAngularPosVelPAcc.md create mode 100644 书籍/力学书籍/OpenFast/FASTKinetics/auto/CalculateForcesMoments.md create mode 100644 书籍/力学书籍/OpenFast/FASTKinetics/auto/CalculatePositions.md create mode 100644 书籍/力学书籍/OpenFast/FASTKinetics/auto/Coeff.md create mode 100644 未命名 1.base diff --git a/.obsidian/plugins/copilot/data.json b/.obsidian/plugins/copilot/data.json index 090a3ea..752fffe 100644 --- a/.obsidian/plugins/copilot/data.json +++ b/.obsidian/plugins/copilot/data.json @@ -316,7 +316,7 @@ "name": "Translate to Chinese", "prompt": "Translate the text below into Chinese:\n 1. Preserve the meaning and tone\n 2. Maintain appropriate cultural context\n 3. Keep formatting and structure\n \n \n\n{copilot-selection}\n\n1. Blade翻译为叶片，flapwise翻译为挥舞，edgewise翻译为摆振，pitch angle翻译成变桨角度，twist angle翻译为扭角，rotor翻译为风轮，turbine、wind turbine翻译为机组、风电机组，span翻译为展向，deflection翻译为变形，mode翻译为模态，normal mode翻译为简正模态，jacket 翻译为导管架，superelement翻译为超单元，shaft翻译为主轴，azimuth、azimuth angle翻译为方位角，neutral axes 翻译为中性轴\n2. Return only the translated text.\n", "showInContextMenu": true, - "modelKey": "gemini-2.5-flash|google" + "modelKey": "ministral-3:14b|ollama" }, { "name": "Summarize", diff --git a/copilot-conversations/20260107_103211__[[Blade]]详细展示公式推导过程.md b/copilot-conversations/20260107_103211__[[Blade]]详细展示公式推导过程.md new file mode 100644 index 0000000..af4ff24 --- /dev/null +++ b/copilot-conversations/20260107_103211__[[Blade]]详细展示公式推导过程.md @@ -0,0 +1,799 @@ +--- +epoch: 1767753131620 +modelKey: gemini-2.5-flash|google +tags: + - copilot-conversation +--- + +**user**: [[Blade]]详细展示公式推导过程 +[Timestamp: 2026/01/07 10:32:11] + +**ai**: 该笔记 [[Blade]] 详细展示了叶片广义力（包括惯性力、弹性力、阻尼力、重力、气动力）的最终公式形式，但并未提供这些公式从基本原理（如拉格朗日方程或虚功原理）出发的详细推导过程。 + +以下是对笔记中各项公式的解释及其背后的物理原理： + +### 1. 广义惯性力 ($F_{r}^{*}\big|_{B1}$) + +广义惯性力是由于物体运动（包括平动和转动）产生的惯性效应，通过虚功原理或拉格朗日方程推导而来。它通常与系统的质量和加速度有关。 + +* **公式 (1):** + $$ + F_{r}^{*}\big|_{B1} = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, {}^E \, \boldsymbol{v}_{r}^{S1}(r) \cdot {}^E \, \boldsymbol{a}^{S1}(r) \, dr - m^{BldTip} \, {}^E \, \boldsymbol{v}_{r}^{S1}(BldFlexL) \cdot {}^E \, \boldsymbol{a}^{S1}(BldFlexL) + $$ + 这个公式表示叶片1的广义惯性力，它由两部分组成： + 1. **积分项:** 沿叶片柔性段长度 $BldFlexL$ 的积分，表示叶片分布式质量的惯性力贡献。其中 $\mu^{B1}(r)$ 是单位长度质量，${}^E\boldsymbol{v}_{r}^{S1}(r)$ 是广义坐标 $q_r$ 对应的虚速度（或模态速度），${}^E\boldsymbol{a}^{S1}(r)$ 是叶片上点 $S1$ 的绝对加速度。点积表示虚功。 + 2. **叶尖质量项:** $m^{BldTip}$ 是叶尖质量，${}^E\boldsymbol{v}_{r}^{S1}(BldFlexL)$ 和 ${}^E\boldsymbol{a}^{S1}(BldFlexL)$ 分别是叶尖处的虚速度和加速度。 + 负号表示惯性力与加速度方向相反（D'Alembert 原理）。 + +* **公式 (2) 和 (3):** + $$ + \mu^{B1}(r) = AdjBlMs^{B1} \cdot BMassDen^{B1}(r) + $$ + $$ + m^{B1Tip} = TipMass(1) + $$ + 这些公式定义了叶片单位长度质量 $\mu^{B1}(r)$ 和叶尖质量 $m^{B1Tip}$，它们可能包含调整系数 ($AdjBlMs^{B1}$) 或直接来自输入参数 ($BMassDen^{B1}(r)$, $TipMass(1)$)。 + +* **公式 (4):** + $$ + F_{r}^{*}\big|_{B1} = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, {}^{E}\boldsymbol{v}_{r}^{S1}(r) \cdot \left\{ \ldots \right\} dr - m^{B1Tip} {}^{E}\boldsymbol{v}_{r}^{S1}(BldFlexL) \cdot \left\{ \ldots \right\} + $$ + 此公式展开了加速度 ${}^E\boldsymbol{a}^{S1}(r)$。在多体动力学中，加速度通常可以分解为与广义加速度 $\ddot{q}_i$ 相关的项和与广义速度 $\dot{q}_i$ 相关的项（通常是科里奥利力和离心力项）。 + * 第一行大括号内是与 $\ddot{q}_i$ 相关的项，用于形成质量矩阵。 + * 第二行大括号内是与 $\dot{q}_i$ 相关的项，通常是科里奥利力和离心力，它们来自虚速度对时间的导数。 + +* **公式 (5):** + $$ + \left[C(q,t)\right]_{B1} (\text{Row,Col}) = \int_{0}^{BldFlexL} \mu^{B1}(r) \, {}^{E}\boldsymbol{v}_{\text{Row}}^{S1}(r) \cdot {}^{E}\boldsymbol{v}_{\text{Col}}^{S1}(r) \, dr + m^{B1Tip} {}^{E}\boldsymbol{v}_{\text{Row}}^{S1}(BldFlexL) \cdot {}^{E}\boldsymbol{v}_{\text{Col}}^{S1}(BldFlexL) + $$ + 这是系统的**广义质量矩阵** $[C(q,t)]_{B1}$ 的元素。它通过将公式 (4) 中与 $\ddot{q}$ 相关的项提取出来，并与虚速度 ${}^{E}\boldsymbol{v}_{\text{Row}}^{S1}(r)$ 进行点积积分得到。矩阵的每个元素代表了不同广义坐标之间的惯性耦合。 + +* **公式 (6):** + $$ + \left\{-f(\dot{q},q,t)\right\}_{B1} (\text{Row}) = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, {}^{E}\boldsymbol{v}_{\text{Row}}^{S1}(r) \cdot \left\{ \ldots \right\} dr - m^{B1Tip} {}^{E}\boldsymbol{v}_{\text{Row}}^{S1}(BldFlexL) \cdot \left\{ \ldots \right\} + $$ + 这是**广义科里奥利力和离心力**项 $\{-f(\dot{q},q,t)\}_{B1}$。它对应于公式 (4) 中与 $\dot{q}$ 相关的项（即 $\frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{S1}(r) \right) \dot{q}_{i}$），这些项在动力学方程中表现为非线性速度相关力。 + +### 2. 广义弹性力 ($F_{r}\big|_{ElasticB1}$) + +广义弹性力来源于叶片变形时储存的弹性势能。 + +* **公式 (7):** + $$ + F_{r}\big|_{ElasticB1} = -\frac{\partial V'^{B1}}{\partial q_r} + $$ + 这是弹性力的基本定义，即弹性势能 $V'^{B1}$ 对广义坐标 $q_r$ 的负偏导数。 + +* **公式 (8):** + $$ + \left. F_{r} \right|_{ElasticB1} = + \begin{cases} + - {k'}_{11}^{B1F} q_{B1F1} - {k'}_{12}^{B1F} q_{B1F2} + & \text{for } r = B1F1 \\ + - {k'}_{11}^{B1E} q_{B1E1} + & \text{for } r = B1E1 \\ + - {k'}_{21}^{B1F} q_{B1F1} - {k'}_{22}^{B1F} q_{B1F2} + & \text{for } r = B1F2 \\ + 0 + & \text{otherwise} + \end{cases} + $$ + 此公式给出了具体的广义弹性力表达式，它与广义位移 $q$ 呈线性关系，并通过广义刚度 $k'$ 连接。这表明叶片被建模为具有多个模态的弹性体，其中 $B1F1, B1F2$ 代表挥舞方向的模态， $B1E1$ 代表摆振方向的模态。 + +* **公式 (9) 和 (10):** + $$ + {k'}_{ij}^{B1F} = \sqrt{FlStTunr^{B1}(i) \, FlStTunr^{B1}(j)} \int_{0}^{BldFlexL} EI^{B1F}(r) \frac{d^2 \phi_i^{B1F}(r)}{dr^2} \frac{d^2 \phi_j^{B1F}(r)}{dr^2} dr + $$ + $$ + EI^{B1F}(r) = AdjFlSt^{B1} \cdot FlpStff^{B1}(r) + $$ + 这些公式定义了挥舞方向的广义刚度 ${k'}_{ij}^{B1F}$。它通过对叶片弯曲刚度 $EI^{B1F}(r)$ 和模态函数 $\phi$ 的二阶导数（曲率）进行积分得到。这来源于梁的弯曲应变能公式 $V = \frac{1}{2} \int EI (\frac{d^2w}{dx^2})^2 dx$。 + +* **公式 (11) 和 (12):** + $$ + {k'}_{11}^{B1E} = \int_{0}^{BldFlexL} EI^{B1E}(r) \left[ \frac{d^2 \phi_{1}^{B1E}(r)}{dr^2} \right]^2 dr + $$ + $$ + EI^{B1E}(r) = AdjEdSt^{B1} \cdot EdgStff^{B1}(r) + $$ + 这些公式定义了摆振方向的广义刚度 ${k'}_{11}^{B1E}$，与挥舞方向的推导类似。 + +### 3. 广义阻尼力 ($F_r\big|_{DampB1}$) + +广义阻尼力用于耗散系统能量，通常与广义速度成正比。这里采用了瑞利阻尼技术。 + +* **公式 (13):** + $$ + \left. F_r \right|_{DampB1} = + \begin{cases} + -\frac{\zeta_1^{B1F} {k'}_{11}^{B1F}}{\pi f^{\prime B 1 F}_1} \dot{q}_{B1F1} - \frac{\zeta_2^{B1F} {k'}_{12}^{B1F}}{\pi f^{\prime B 1 F}_2} \dot{q}_{B1F2} & \text{for } r = B1F1 \\ + -\frac{\zeta_1^{B1E} {k'}_{11}^{B1E}}{\pi f^{\prime B 1 E}_1} \dot{q}_{B1E1} & \text{for } r = B1E1 \\ + -\frac{\zeta_1^{B1F} {k'}_{21}^{B1F}}{\pi f^{\prime B 1 F}_1} \dot{q}_{B1F1} - \frac{\zeta_2^{B1F} {k'}_{22}^{B1F}}{\pi f^{\prime B 1 F}_2} \dot{q}_{B1F2} & \text{for } r = B1F2 \\ + 0 & \text{otherwise} + \end{cases} + $$ + 此公式给出了广义阻尼力，它与广义速度 $\dot{q}$ 呈线性关系。阻尼系数由结构阻尼比 $\zeta$、广义刚度 $k'$ 和固有频率 $f'$ 确定。瑞利阻尼假设阻尼矩阵与质量矩阵和/或刚度矩阵成比例。 + +* **公式 (14) 和 (15):** + $$ + f_i^{\prime B1F} = \frac{1}{2\pi} \sqrt{\frac{k_{ii}^{\prime B1F}}{m_{ii}^{\prime B1F}}} + $$ + $$ + f_i^{\prime B1E} = \frac{1}{2\pi} \sqrt{\frac{k_{ii}^{\prime B1E}}{m_{ii}^{\prime B1E}}} + $$ + 这些是单自由度系统固有频率的经典公式，其中 $k_{ii}^{\prime}$ 是广义刚度， $m_{ii}^{\prime}$ 是广义质量。 + +* **公式 (16) 和 (17):** + $$ + m_{ij}^{\prime B1F} = \int_{0}^{BldFlexL} \mu^{B1}(r) \phi_i^{B1F}(r) \phi_j^{B1F}(r) \, dr + $$ + $$ + m_{11}^{\prime B1E} = \int_{0}^{BldFlexL} \mu^{B1}(r) \left[ \phi_1^{B1E}(r) \right]^2 dr + $$ + 这些公式定义了广义质量 $m_{ij}^{\prime}$，它通过对叶片单位长度质量 $\mu^{B1}(r)$ 和模态函数 $\phi$ 的乘积进行积分得到。这来源于动能公式 $T = \frac{1}{2} \int \mu \dot{w}^2 dx$。 + +* **公式 (18), (19), (20):** + $$ + [C(q,t)]|_{ElasticB1} = 0 \quad \text{and} \quad [C(q,t)]|_{DampB1} = 0 + $$ + $$ + \left. \left\{ -f(\dot{q},q,t) \right\} \right|_{ElasticB1} = \begin{pmatrix} \vdots \\ \vdots \\ \vdots \\ -k_{11}^{\prime B1F} q_{B1F1} - k_{12}^{\prime B1F} q_{B1F2} \\ -k_{11}^{\prime B1E} q_{B1E1} \\ -k_{21}^{\prime B1F} q_{B1F1} - k_{22}^{\prime B1F} q_{B1F2} \\ \vdots \\ \vdots \\ \vdots \end{pmatrix} + $$ + $$ + \left. \left\{ -f(\dot{q},q,t) \right\} \right|_{DampB1} = \begin{pmatrix} \vdots \\ \vdots \\ \vdots \\ -\frac{\zeta_1^{B1F} k_{11}^{\prime B1F}}{\pi f_1^{\prime B1F}} \dot{q}_{B1F1} - \frac{\zeta_2^{B1F} k_{12}^{\prime B1F}}{\pi f_2^{\prime B1F}} \dot{q}_{B1F2} \\ -\frac{\zeta_1^{B1E} k_{11}^{\prime B1E}}{\pi f_1^{\prime B1E}} \dot{q}_{B1E1} \\ -\frac{\zeta_1^{B1F} k_{21}^{\prime B1F}}{\pi f_1^{\prime B1F}} \dot{q}_{B1F1} - \frac{\zeta_2^{B1F} k_{22}^{\prime B1F}}{\pi f_2^{\prime B1F}} \dot{q}_{B1F2} \\ \vdots \\ \vdots \\ \vdots \end{pmatrix} + $$ + 这些公式总结了弹性力和阻尼力在动力学方程中的作用。它们不贡献于广义质量矩阵 $[C(q,t)]$，而是作为广义力向量 $\{-f(\dot{q},q,t)\}$ 的一部分出现。 + +### 4. 广义重力 ($F_{r}\big|_{GravB1}$) + +广义重力来源于叶片在重力场中的势能。 + +* **公式 (21):** + $$ + F_{r}\big|_{GravB1} = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, \boldsymbol{g}^{E}\boldsymbol{v}_{r}^{S1}(r) \cdot \boldsymbol{{z}}_2 \, dr - m^{B1Tip} \boldsymbol{g}^{E}\boldsymbol{v}_{r}^{S1}(BldFlexL) \cdot \boldsymbol{{z}}_2 + $$ + 此公式表示广义重力，它由叶片分布式质量和叶尖质量的重力贡献组成。其中 $\boldsymbol{g}$ 是重力加速度，$\boldsymbol{z}_2$ 是重力方向的单位向量。点积表示重力在虚位移方向上所做的虚功。 + +* **公式 (22) 和 (23):力加速度，$\boldsymbol{z}_2$ 是重力方向的单位向量。点积表示重力在虚位移方向上所做的虚功。 + +* **公式 (22) 和 (23): +[Timestamp: 2026/01/07 10:32:38] + +**user**: 继续生成 +[Timestamp: 2026/01/07 10:52:38] + +**ai**: 好的，我将继续对笔记 [[Blade]] 中剩余的公式进行详细解释。 + +### 4. 广义重力 ($F_{r}\big|_{GravB1}$) (续) + +* **公式 (22):** + $$ + [C(q,t)]|_{GravB1} = 0 + $$ + 这个公式表明广义重力不贡献于广义质量矩阵 $[C(q,t)]$。重力是一个保守力，其作用只与位置有关，不直接与加速度相关。 + +* **公式 (23):** + $$ + \left\{ -f(\dot{q}, q, t) \right\}|_{GravB1} (\text{Row}) = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, g \, {}^E\boldsymbol{v}_{Row}^{S1}(r) \cdot \boldsymbol{z}_2 \, dr - m^{B1Tip} \, g \, {}^E\boldsymbol{v}_{Row}^{S1}(BldFlexL) \cdot \boldsymbol{z}_2 + $$ + 此公式给出了广义重力项 $\{-f(\dot{q}, q, t)\}|_{GravB1}$。它表示重力在广义坐标方向上的投影。其中 $g$ 是重力加速度的大小，$\boldsymbol{z}_2$ 是重力方向的单位向量。这个项是广义力向量的一部分，在动力学方程中作为外力项出现。 + +### 5. 广义气动力 ($\boldsymbol{F}_r \big|_{AeroB1}$) + +广义气动力是由于叶片与空气相互作用产生的力，包括升力、阻力以及气动俯仰力矩。 + +* **公式 (24):** + $$ + \boldsymbol{F}_r \big|_{AeroB1} = \int_{0}^{BldFlexL} \left[ {}^E\boldsymbol{v}_r^{S1}(r) \cdot \boldsymbol{F}_{AeroB1}^{S1}(r) + {}^E\boldsymbol{\omega}_r^{M1}(r) \cdot \boldsymbol{M}_{AeroB1}^{M1}(r) \right] dr + {}^E\boldsymbol{v}_r^{S1}(BldFlexL) \cdot \boldsymbol{F}_{TipDragB1}^{S1}(BldFlexL) + $$ + 这个公式表示叶片1的广义气动力，它由三部分组成： + 1. **分布式气动力项:** 沿叶片柔性段长度 $BldFlexL$ 的积分项。 + * ${}^E\boldsymbol{v}_r^{S1}(r) \cdot \boldsymbol{F}_{AeroB1}^{S1}(r)$: 表示单位展向长度上的气动力 $\boldsymbol{F}_{AeroB1}^{S1}(r)$ 在广义坐标 $q_r$ 对应的虚速度 ${}^E\boldsymbol{v}_r^{S1}(r)$ 方向上所做的虚功。 + * ${}^E\boldsymbol{\omega}_r^{M1}(r) \cdot \boldsymbol{M}_{AeroB1}^{M1}(r)$: 表示单位展向长度上的气动俯仰力矩 $\boldsymbol{M}_{AeroB1}^{M1}(r)$ 在广义坐标 $q_r$ 对应的虚角速度 ${}^E\boldsymbol{\omega}_r^{M1}(r)$ 方向上所做的虚功。 + 2. **叶尖阻力项:** ${}^E\boldsymbol{v}_r^{S1}(BldFlexL) \cdot \boldsymbol{F}_{TipDragB1}^{S1}(BldFlexL)$ 表示叶尖阻力 $\boldsymbol{F}_{TipDragB1}^{S1}(BldFlexL)$ 在叶尖虚速度方向上所做的虚功。 + + **物理意义:** 广义气动力是通过将作用在叶片上的实际气动力和气动扭矩（包括分布式载荷和叶尖载荷）投影到各个广义坐标的虚位移（或虚速度）方向上得到的。这符合虚功原理，即广义力等于实际力在虚位移上所做的虚功。 + + **其中描述:** + * $F_{AeroB1}^{S1}(r)$ 和 $M_{AeroB1}^{M1}(r)$ 是作用于叶片1上点S1处的单位展向长度上的气动力和气动俯仰力矩。 + * $M_{AeroB1}^{M1}(r)$ 可以包含直接气动俯仰力矩（如由翼型弯度引起的Cm）和由气动升力/阻力作用在偏离叶片元素质心的位置所产生的力矩。 + +* **公式 (25):** + $$ + [C(q,t)]|_{AeroB1} = 0 + $$ + 这个公式表明广义气动力不贡献于广义质量矩阵 $[C(q,t)]$。气动力通常被视为外力，其大小和方向取决于叶片的运动状态（速度、攻角等），但不直接与加速度线性相关以形成质量矩阵的项。 + +* **公式 (26):** + $$ + \left\{ -f(\dot{q}, q, t) \right\}|_{AeroB1} (\text{Row}) = \int_{0}^{BldFlexL} \left[ {}^E\boldsymbol{v}_{Row}^{S1}(r) \cdot \boldsymbol{F}_{AeroB1}^{S1}(r) + {}^E\boldsymbol{\omega}_{Row}^{M1}(r) \cdot \boldsymbol{M}_{AeroB1}^{M1}(r) \right] dr + {}^E\boldsymbol{v}_{Row}^{S1}(BldFlexL) \cdot \boldsymbol{F}_{TipDragB1}^{S1}(BldFlexL) + $$ + 此公式给出了广义气动力项 $\{-f(\dot{q}, q, t)\}|_{AeroB1}$。它与公式 (24) 的形式相同，表明广义气动力作为广义力向量的一部分，在动力学方程中作为外力项出现。气动力通常是速度和位移的复杂非线性函数。 + +### 总结 + +这份笔记 [[Blade]] 详细列出了风力机叶片动力学模型中各种广义力的最终表达式。这些表达式是基于多体动力学理论（如虚功原理或拉格朗日方程）推导而来的。 + +**推导过程的通用步骤通常包括：** + +1. **定义系统和广义坐标:** 确定叶片的柔性模态（挥舞、摆振）以及其他自由度（如桨距、偏航等），并选择一组广义坐标 $q_i$ 来描述系统的运动。 +2. **建立运动学关系:** 确定叶片上任意一点的位移、速度和加速度与广义坐标及其导数之间的关系。这通常涉及坐标系的转换和模态函数的应用。 +3. **计算动能和势能:** + * **动能 ($T$):** 积分叶片分布式质量和叶尖质量的动能。 + * **弹性势能 ($V_{elastic}$):** 积分叶片变形（弯曲、扭转）产生的应变能。 + * **重力势能 ($V_{gravity}$):** 积分叶片分布式质量和叶尖质量在重力场中的势能。 +4. **应用拉格朗日方程或虚功原理:** + * **拉格朗日方程:** $\frac{d}{dt}\left(\frac{\partial T}{\partial \dot{q}_r}\right) - \frac{\partial T}{\partial q_r} + \frac{\partial V}{\partial q_r} = Q_r$ + * 其中 $Q_r$ 是非保守广义力（如阻尼力、气动力）。 + * 通过对动能和势能求导，可以得到广义惯性力、广义弹性力和广义重力。 + * **虚功原理:** $\delta W = \sum_i (\boldsymbol{F}_i - m_i \boldsymbol{a}_i) \cdot \delta \boldsymbol{r}_i + \sum_j (\boldsymbol{M}_j - \boldsymbol{H}_j) \cdot \delta \boldsymbol{\theta}_j = 0$ + * 将所有主动力、惯性力、阻尼力等在虚位移上所做的虚功求和，并令其为零。通过提取广义坐标前的系数，可以得到广义力。 + +笔记中给出的公式是这些推导过程的最终结果，它们将叶片的物理特性（质量分布、刚度分布、阻尼比、模态形状）与广义坐标及其导数联系起来，形成了用于风力机叶片动力学分析的运动方程。 +[Timestamp: 2026/01/07 10:52:51] + +**user**: 基于凯恩方程原理，对上述公式做重新解释 +[Timestamp: 2026/01/07 10:55:38] + +**ai**: 这份笔记 [[Blade]] 中展示的公式与凯恩方程（Kane's Equations）原理高度一致。凯恩方程是一种系统地推导复杂机械系统运动方程的方法，尤其适用于多体系统和柔性体动力学。其核心思想是将系统的运动方程表示为广义主动力 ($F_r$) 和广义惯性力 ($F_r^*$) 之和为零的形式： + +$F_r + F_r^* = 0 \quad (r = 1, \ldots, n)$ + +其中，$n$ 是广义坐标的数量。 + +下面我将基于凯恩方程的原理，对笔记中的各项公式进行重新解释： + +### 凯恩方程基础概念回顾 + +1. **广义坐标 ($q_r$) 和广义速度 ($u_r$):** 描述系统构型的独立变量及其时间导数。在笔记中，广义坐标 $q_r$ 及其时间导数 $\dot{q}_r$ 被直接用作广义速度。 +2. **偏速度 (Partial Velocity) 和偏角速度 (Partial Angular Velocity):** 对于系统中的任意一点或刚体，其绝对速度 $\boldsymbol{v}$ 和绝对角速度 $\boldsymbol{\omega}$ 可以表示为广义速度的线性组合： + $\boldsymbol{v} = \sum_r \boldsymbol{v}_r \dot{q}_r + \boldsymbol{v}_t$ + $\boldsymbol{\omega} = \sum_r \boldsymbol{\omega}_r \dot{q}_r + \boldsymbol{\omega}_t$ + 其中，$\boldsymbol{v}_r$ 和 $\boldsymbol{\omega}_r$ 分别是该点或刚体相对于广义速度 $\dot{q}_r$ 的偏速度和偏角速度。笔记中的 ${}^E\boldsymbol{v}_{r}^{S1}(r)$ 和 ${}^E\boldsymbol{\omega}_{r}^{M1}(r)$ 正是这些偏速度和偏角速度。 +3. **广义主动力 ($F_r$):** 所有作用在系统上的**主动力**（如弹性力、阻尼力、重力、气动力等）在对应偏速度和偏角速度方向上的投影之和。 + $F_r = \sum_{k} (\boldsymbol{v}_r^{(k)} \cdot \boldsymbol{R}^{(k)} + \boldsymbol{\omega}_r^{(k)} \cdot \boldsymbol{M}^{(k)})$ + 其中 $\boldsymbol{R}^{(k)}$ 和 $\boldsymbol{M}^{(k)}$ 是作用在 $k$ 点的力和力矩。 +4. **广义惯性力 ($F_r^*$):** 所有作用在系统上的**惯性力**（由质量和加速度引起）在对应偏速度和偏角速度方向上的投影之和。 + $F_r^* = \sum_{k} (\boldsymbol{v}_r^{(k)} \cdot \boldsymbol{R}^{*(k)} + \boldsymbol{\omega}_r^{(k)} \cdot \boldsymbol{M}^{*(k)})$ + 其中 $\boldsymbol{R}^{*(k)} = -m^{(k)} \boldsymbol{a}^{(k)}$ 是惯性力，$\boldsymbol{M}^{*(k)}$ 是惯性力矩。 + +### 笔记公式的凯恩方程解释 + +#### 1. 广义惯性力 ($F_{r}^{*}\big|_{B1}$) + +凯恩方程中，广义惯性力是系统各部分惯性力在偏速度方向上的投影之和。 + +* **公式 (1):** + $$ + F_{r}^{*}\big|_{B1} = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, {}^E \, \boldsymbol{v}_{r}^{S1}(r) \cdot {}^E \, \boldsymbol{a}^{S1}(r) \, dr - m^{BldTip} \, {}^E \, \boldsymbol{v}_{r}^{S1}(BldFlexL) \cdot {}^E \, \boldsymbol{a}^{S1}(BldFlexL) + $$ + 这个公式直接体现了凯恩方程中广义惯性力的定义。 + * ${}^E \, \boldsymbol{v}_{r}^{S1}(r)$ 是叶片上径向位置 $r$ 处点 $S1$ 相对于广义坐标 $q_r$ 的**偏速度**。 + * $\mu^{B1}(r) \, dr$ 是微分质量元素。 + * $\mu^{B1}(r) \, {}^E \, \boldsymbol{a}^{S1}(r) \, dr$ 是作用在微分质量元素上的惯性力（质量乘以加速度）。 + * 点积 $\boldsymbol{v}_r \cdot (m \boldsymbol{a})$ 表示惯性力在偏速度方向上的投影。负号表示惯性力与加速度方向相反。 + * 积分项处理分布式质量的贡献，第二项处理叶尖集中质量的贡献。 + +* **公式 (4):** 此公式将绝对加速度 ${}^E\boldsymbol{a}^{S1}(r)$ 分解为与广义加速度 $\ddot{q}_i$ 相关的项和与广义速度 $\dot{q}_i$ 相关的项。 + * 与 $\ddot{q}_i$ 相关的项用于构建广义质量矩阵。 + * 与 $\dot{q}_i$ 相关的项（涉及偏速度对时间的导数 $\frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{S1}(r) \right)$）代表了**科里奥利加速度和离心加速度**等速度相关效应。这些项在凯恩方程中构成了广义惯性力中与速度相关的部分。 + +* **公式 (5) - 广义质量矩阵 $[C(q,t)]_{B1}$:** + 在凯恩方程中，当广义惯性力 $F_r^*$ 被写成 $\sum_s M_{rs} \ddot{q}_s + G_r$ 的形式时，$[C(q,t)]_{B1}$ 正是其中的 $M_{rs}$ 元素，即**广义质量矩阵**。它是由偏速度与偏速度的点积（加权质量）形成的，反映了系统在加速时的惯性耦合。 + +* **公式 (6) - 速度相关惯性力 $\{-f(\dot{q},q,t)\}_{B1}$:** + 此项对应于凯恩方程中广义惯性力 $F_r^*$ 的 $G_r$ 部分，包含了所有**速度相关的惯性效应**，如科里奥利力和离心力。这些力是非线性的，并且在运动方程中作为广义力向量的一部分出现。 + +#### 2. 广义弹性力 ($F_{r}\big|_{ElasticB1}$) + +* **公式 (7):** $F_{r}\big|_{ElasticB1} = -\frac{\partial V'^{B1}}{\partial q_r}$。 + 在凯恩方程中，保守力（如弹性力）的广义力可以通过其势能 $V'$ 对广义坐标的负偏导数来计算。这是一种将保守力纳入广义主动力 $F_r$ 的标准方法。 +* **公式 (8):** 提供了具体的线性弹性广义力表达式。这些力直接贡献于广义主动力 $F_r$。 +* **公式 (9)-(12):** 定义了广义刚度 $k'$，这些是弹性力表达式中的系数，来源于叶片的材料特性和模态形状。 + +#### 3. 广义阻尼力 ($F_r\big|_{DampB1}$) + +* **公式 (13):** 提供了广义阻尼力的表达式。 + 在凯恩方程中，阻尼力是非保守力，它们直接作为广义主动力 $F_r$ 的一部分。此处的瑞利阻尼假设使得阻尼力与广义速度成正比。 +* **公式 (14)-(17):** 定义了固有频率 $f'$ 和广义质量 $m'$，这些参数用于计算阻尼系数，进而确定广义阻尼力。 + +#### 4. 广义重力 ($F_{r}\big|_{GravB1}$) + +* **公式 (21):** + $$ + F_{r}\big|_{GravB1} = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, \boldsymbol{g}^{E}\boldsymbol{v}_{r}^{S1}(r) \cdot \boldsymbol{{z}}_2 \, dr - m^{B1Tip} \boldsymbol{g}^{E}\boldsymbol{v}_{r}^{S1}(BldFlexL) \cdot \boldsymbol{{z}}_2 + $$ + 此公式计算了广义重力。 + * $\mu^{B1}(r) \, \boldsymbol{g}$ 是微分质量元素上的重力。 + * ${}^E\boldsymbol{v}_{r}^{S1}(r) \cdot \boldsymbol{{z}}_2$ 表示重力在偏速度方向上的投影。 + * 这直接符合凯恩方程中广义主动力 $F_r$ 的定义，即主动力在偏速度方向上的投影。 + +#### 5. 广义气动力 ($\boldsymbol{F}_r \big|_{AeroB1}$) + +* **公式 (24):** + $$ + \boldsymbol{F}_r \big|_{AeroB1} = \int_{0}^{BldFlexL} \left[ {}^E\boldsymbol{v}_r^{S1}(r) \cdot \boldsymbol{F}_{AeroB1}^{S1}(r) + {}^E\boldsymbol{\omega}_r^{M1}(r) \cdot \boldsymbol{M}_{AeroB1}^{M1}(r) \right] dr + {}^E\boldsymbol{v}_r^{S1}(BldFlexL) \cdot \boldsymbol{F}_{TipDragB1}^{S1}(BldFlexL) + $$ + 这个公式是凯恩方程中广义主动力计算的典型示例。 + * $\boldsymbol{F}_{AeroB1}^{S1}(r)$ 和 $\boldsymbol{M}_{AeroB1}^{M1}(r)$ 是实际的气动力和气动俯仰力矩。 + * ${}^E\boldsymbol{v}_r^{S1}(r)$ 和 ${}^E\boldsymbol{\omega}_r^{M1}(r)$ 分别是作用点处的**偏速度**和**偏角速度**。 + * 通过将实际气动力和气动扭矩投影到对应的偏速度和偏角速度方向上，直接计算出广义气动力，这正是凯恩方程中广义主动力 $F_r$ 的定义。 + +### 总结 + +笔记 [[Blade]] 中的所有公式都完美契合凯恩方程的框架。它将复杂的叶片动力学问题分解为广义惯性力 ($F_r^*$) 和广义主动力 ($F_r$) 的计算，其中： + +* **广义惯性力**被进一步分解为广义质量矩阵 $[C(q,t)]$ 和速度相关的惯性力 $\{-f(\dot{q},q,t)\}_{B1}$。 +* **广义主动力**则包含了弹性力、阻尼力、重力、气动力等各项贡献。 + +最终，将这些项组合起来，就可以得到叶片的运动方程，其形式为： + +$\mathbf{M}(q,t)\ddot{\mathbf{q}} + \mathbf{g}(\dot{\mathbf{q}},\mathbf{q},t) = \mathbf{Q}(\dot{\mathbf{q}},\mathbf{q},t)$ + +其中，$\mathbf{M}(q,t)$ 是广义质量矩阵，$\mathbf{g}(\dot{\mathbf{q}},\mathbf{q},t)$ 包含了所有速度相关的惯性力，而 $\mathbf{Q}(\dot{\mathbf{q}},\mathbf{q},t)$ 则是所有广义主动力（弹性、阻尼、重力、气动力）的总和。这种结构正是凯恩方程在柔性多体系统动力学建模中的典型应用。 +[Timestamp: 2026/01/07 10:56:10] + +**user**: [[Blade]][[FillAugmat]]结合理论和代码，解释叶片部分是如何计算的 +[Timestamp: 2026/01/07 11:08:54] + +**ai**: 笔记 [[Blade]] 详细阐述了风力机叶片广义力（惯性力、弹性力、阻尼力、重力、气动力）的理论公式。而 `FillAugMat` 子程序（来自 [[FillAugmat]] 笔记）则是将这些理论公式转化为实际计算代码，用于构建系统的运动方程。 + +`FillAugMat` 的核心任务是填充一个**增广矩阵 (Augmented Matrix)** `AugMat`。这个矩阵通常用于表示系统的动力学方程： + +$\mathbf{M}(q,t)\ddot{\mathbf{q}} = \mathbf{F}(\dot{q},q,t)$ + +其中： +* $\mathbf{M}(q,t)$ 是广义质量矩阵，它包含了与广义加速度 $\ddot{\mathbf{q}}$ 相关的项。 +* $\mathbf{F}(\dot{q},q,t)$ 是广义力向量，它包含了所有主动力（弹性力、阻尼力、重力、气动力）以及速度相关的惯性力（科里奥利力、离心力）。 + +在 `FillAugMat` 中，`AugMat` 的前 `N` 列（假设有 `N` 个自由度）构成了广义质量矩阵 $\mathbf{M}$，而最后一列（通过 `p%NAug` 索引）则存储了广义力向量 $\mathbf{F}$。 + +下面将结合 [[Blade]] 中的理论公式和 `FillAugMat` 中的代码，解释叶片部分的计算过程。 + +### 叶片部分的计算 (`DO K = 1,p%NumBl` 循环内) + +`FillAugMat` 子程序首先会遍历所有叶片（`DO K = 1,p%NumBl`），对每个叶片计算其对系统运动方程的贡献。 + +#### 1. 广义惯性力 (Generalized Inertia Forces) + +凯恩方程中的广义惯性力 $F_r^*$ 包含两部分：与广义加速度 $\ddot{q}$ 相关的项（构成广义质量矩阵）和与广义速度 $\dot{q}$ 相关的项（构成广义力向量的一部分）。 + +* **理论依据:** [[Blade]] 中的公式 (1), (4), (5), (6)。 + * 公式 (5) 定义了广义质量矩阵 $[C(q,t)]_{B1}$ 的元素，它与 $\ddot{q}$ 相关。 + * 公式 (6) 定义了速度相关的惯性力 $\{-f(\dot{q},q,t)\}_{B1}$，它与 $\dot{q}$ 相关。 + +* **代码实现:** + * **叶尖质量的贡献:** + ```fortran + 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,:), & + RtHSdat%PLinVelES(K, p%TipNode, p%DOFs%PSBE(K,L),0,:) ) + ``` + 这行代码计算了叶尖质量对**广义质量矩阵**的贡献。它对应于 [[Blade]] 公式 (5) 中的叶尖质量项 $m^{B1Tip} {}^{E}\boldsymbol{v}_{\text{Row}}^{S1}(BldFlexL) \cdot {}^{E}\boldsymbol{v}_{\text{Col}}^{S1}(BldFlexL)$。 + * `p%TipMass(K)` 是叶尖质量 $m^{B1Tip}$。 + * `RtHSdat%PLinVelES(...)` 是偏速度 ${}^{E}\boldsymbol{v}_{\text{Row}}^{S1}$。 + * **分布式叶片元素的贡献:** + ```fortran + AugMat(p%DOFs%PSBE(K,I),p%DOFs%PSBE(K,L)) = AugMat(...) + p%BElmntMass(J,K)*& + DOT_PRODUCT( RtHSdat%PLinVelES(K,J,p%DOFs%PSBE(K,I),0,:), & + RtHSdat%PLinVelES(K,J,p%DOFs%PSBE(K,L),0,:) ) + ``` + 这行代码累加了每个叶片元素质量对**广义质量矩阵**的贡献。它对应于 [[Blade]] 公式 (5) 中的积分项 $\int_{0}^{BldFlexL} \mu^{B1}(r) \, {}^{E}\boldsymbol{v}_{\text{Row}}^{S1}(r) \cdot {}^{E}\boldsymbol{v}_{\text{Col}}^{S1}(r) \, dr$。 + * `p%BElmntMass(J,K)` 是叶片元素 $J$ 的质量（对应 $\mu^{B1}(r) \, dr$）。 + * `RtHSdat%PLinVelES(...)` 是偏速度 ${}^{E}\boldsymbol{v}_{\text{Row}}^{S1}(r)$。 + * **速度相关的惯性力 (科里奥利力、离心力):** + 在 `TmpVec1` 的计算中，包含了 `p%TipMass(K)*RtHSdat%LinAccESt(:,K,p%TipNode)` 和 `p%BElmntMass(J,K)*RtHSdat%LinAccESt(:,K,J)`。这里的 `RtHSdat%LinAccESt` 提供了叶片上点的绝对加速度，其中包含了速度相关的项。当这些项与偏速度进行点积并累加到 `AugMat(..., p%NAug)` 时，就形成了广义力向量中速度相关的惯性力部分，对应于 [[Blade]] 公式 (6)。 + +#### 2. 广义弹性力 (Generalized Elastic Forces) + +* **理论依据:** [[Blade]] 中的公式 (7), (8), (19)。 + * 弹性力是保守力，其广义力通过势能对广义坐标的负偏导数得到。 + * 公式 (8) 和 (19) 给出了具体的线性弹性力表达式，与广义位移 $q$ 成正比。 + +* **代码实现:** + ```fortran + 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)) & + ... + ENDIF + ``` + 这些代码块直接将弹性力项添加到广义力向量 `AugMat(..., p%NAug)` 中。 + * `p%KBF(K,i,j)` 和 `p%KBE(K,i,j)` 是叶片 $K$ 的广义刚度系数，对应于 [[Blade]] 公式 (8) 中的 ${k'}_{ij}^{B1F}$ 和 ${k'}_{11}^{B1E}$。 + * `x%QT(DOF)` 是广义坐标 $q$。 + * 负号表示弹性力是恢复力。 + +#### 3. 广义阻尼力 (Generalized Damping Forces) + +* **理论依据:** [[Blade]] 中的公式 (13), (20)。 + * 阻尼力是非保守力，与广义速度 $\dot{q}$ 成正比。 + +* **代码实现:** + ```fortran + IF ( p%DOF_Flag(DOF_BF(K,1)) ) THEN + AugMat( DOF_BF(K,1),p%NAug) = AugMat(DOF_BF(K,1),p%NAug) & + ... + - p%CBF(K,1,1)*x%QDT(DOF_BF(K,1)) & + - p%CBF(K,1,2)*x%QDT(DOF_BF(K,2)) + ENDIF + ``` + 这些代码块直接将阻尼力项添加到广义力向量 `AugMat(..., p%NAug)` 中。 + * `p%CBF(K,i,j)` 和 `p%CBE(K,i,j)` 是叶片 $K$ 的广义阻尼系数。 + * `x%QDT(DOF)` 是广义速度 $\dot{q}$。 + * 负号表示阻尼力与速度方向相反。 + +#### 4. 广义重力 (Generalized Gravity Forces) + +* **理论依据:** [[Blade]] 中的公式 (21), (23)。 + * 重力是保守力，其广义力通过将重力投影到偏速度方向上得到。 + +* **代码实现:** + * 在 `TmpVec1` 的计算中，包含了重力项： + ```fortran + TmpVec1 = RtHSdat%FSTipDrag(:,K) - p%TipMass(K)*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccESt(:,K,p%TipNode) ) + ... + TmpVec1 = RtHSdat%FSAero(:,K,J)*p%DRNodes(J) - p%BElmntMass(J,K)*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccESt(:,K,J) ) + ``` + * `p%Gravity*CoordSys%z2` 代表重力加速度向量 $\boldsymbol{g}\boldsymbol{z}_2$。 + * `p%TipMass(K)*( p%Gravity*CoordSys%z2 )` 和 `p%BElmntMass(J,K)*( p%Gravity*CoordSys%z2 )` 分别是叶尖和叶片元素的重力。 + * 这些重力项随后通过点积 `DOT_PRODUCT( RtHSdat%PLinVelES(...), TmpVec1 )` 被投影到偏速度方向上，并累加到 `AugMat(..., p%NAug)`，对应于 [[Blade]] 公式 (23)。 + +#### 5. 广义气动力 (Generalized Aerodynamic Forces) + +* **理论依据:** [[Blade]] 中的公式 (24), (26)。 + * 气动力和气动扭矩是非保守力，其广义力通过将它们投影到偏速度和偏角速度方向上得到。 + +* **代码实现:** + * 在 `TmpVec1` 和 `TmpVec3` 的计算中，包含了气动力和气动扭矩项： + ```fortran + TmpVec1 = RtHSdat%FSTipDrag(:,K) - ... ! Tip drag + ... + TmpVec1 = RtHSdat%FSAero(:,K,J)*p%DRNodes(J) - ... ! Distributed aerodynamic force + TmpVec3 = RtHSdat%MMAero(:,K,J)*p%DRNodes(J) ! Distributed aerodynamic moment + ``` + * `RtHSdat%FSTipDrag(:,K)` 是叶尖阻力 $\boldsymbol{F}_{TipDragB1}^{S1}(BldFlexL)$。 + * `RtHSdat%FSAero(:,K,J)*p%DRNodes(J)` 是叶片元素 $J$ 上的气动力 $\boldsymbol{F}_{AeroB1}^{S1}(r) \, dr$。 + * `RtHSdat%MMAero(:,K,J)*p%DRNodes(J)` 是叶片元素 $J$ 上的气动扭矩 $\boldsymbol{M}_{AeroB1}^{M1}(r) \, dr$。 + * 这些力/扭矩随后通过点积 `DOT_PRODUCT( RtHSdat%PLinVelES(...), TmpVec1 )` 和 `DOT_PRODUCT( RtHSdat%PAngVelEM(...), TmpVec3 )` 被投影到偏速度和偏角速度方向上，并累加到 `AugMat(..., p%NAug)`，对应于 [[Blade]] 公式 (26)。 + +### 总结 + +`FillAugMat` 子程序通过迭代叶片的各个部分（叶尖、分布式元素），并利用预先计算好的偏速度、偏角速度、绝对加速度以及各种力/力矩数据 (`RtHSdat` 结构体)，系统地构建了风力机叶片部分的广义质量矩阵和广义力向量。它精确地将 [[Blade]] 中定义的理论公式（基于凯恩方程原理）转化为可计算的数值形式，从而为整个风力机系统的动力学仿真奠定基础。 + +值得注意的是，代码中还包含了对质量矩阵的对称性处理（先填充下三角，再镜像到上三角），以及其他系统组件（如塔架、传动链、偏航系统）的计算，但上述解释主要聚焦于叶片部分的实现。 +[Timestamp: 2026/01/07 11:09:28] + +**user**: [[Blade]][[CalculateForcesMoments]]结合理论和代码，解释叶片部分是如何计算的 +[Timestamp: 2026/01/07 14:04:15] + +**ai**: `CalculateForcesMoments` 子程序在风力机动力学仿真中扮演着关键角色，它负责计算系统中各个部件（包括叶片）所受到的物理力和力矩。这些计算结果随后会被 `FillAugMat` 子程序（如前所述）用于构建系统的广义质量矩阵和广义力向量。 + +本解释将聚焦于 `CalculateForcesMoments` 中与叶片相关的计算部分，并将其与 [[Blade]] 笔记中的理论公式进行对照。 + +### `CalculateForcesMoments` 子程序中的叶片部分计算 + +`CalculateForcesMoments` 子程序首先遍历所有叶片（`DO K = 1,p%NumBl`），然后对每个叶片上的各个节点/元素进行计算（`DO J = 1,p%BldNodes`）。 + +#### 1. 气动力和气动扭矩的计算 (`FSAero`, `MMAero`, `FSTipDrag`) + +这部分代码主要从外部气动模块（通过 `u` 输入）获取气动载荷，并将其转换为内部使用的格式。 + +* **理论依据:** [[Blade]] 中的公式 (24) 描述了广义气动力，其中包含了单位展向长度上的气动力 $\boldsymbol{F}_{AeroB1}^{S1}(r)$ 和气动俯仰力矩 $\boldsymbol{M}_{AeroB1}^{M1}(r)$，以及叶尖阻力 $\boldsymbol{F}_{TipDragB1}^{S1}(BldFlexL)$。 + +* **代码实现:** + * **叶尖阻力 (`RtHSdat%FSTipDrag`):** + ```fortran + RtHSdat%FSTipDrag = 0.0_ReKi ! Calculate the tip drag forces if necessary + !bjj: add this back when we've figured out how to handle the tip brakes: + !RtHSdat%FSTipDrag(:,K) = m%CoordSys%m2(K,p%BldNodes,:)*SIGN( 0.5*p%AirDens*(RtHSdat%LinVelESm2(K)**2)*u%TBDrCon(K), -1.*RtHSdat%LinVelESm2(K) ) + ``` + 目前代码中将叶尖阻力初始化为零，并注释掉了实际的计算逻辑。这表明在当前版本的实现中，叶尖阻力可能被简化或通过其他方式处理。如果启用，它将直接对应 [[Blade]] 公式 (24) 中的叶尖阻力项 $\boldsymbol{F}_{TipDragB1}^{S1}(BldFlexL)$。 + * **气动中心位置 (`RtHSdat%rSAerCen`):** + ```fortran + RtHSdat%rSAerCen(:,J,K) = p%rSAerCenn1(K,J)*CoordSys%n1(K,J,:) + p%rSAerCenn2(K,J)*CoordSys%n2(K,J,:) + ``` + 这行代码计算了叶片元素 $J$ 上点 $S$ 到其气动中心的位置向量 `rSAerCen`。这是计算气动俯仰力矩的关键，因为力矩通常是力与力臂的叉积。 + * **分布式气动力 (`RtHSdat%FSAero`) 和气动扭矩 (`RtHSdat%MMAero`):** + ```fortran + if (p%UseAD14) then + RtHSdat%FSAero(:,K,J) = ( u%BladePtLoads(K)%Force(1,NodeNum) * CoordSys%te1(K,J,:) & + + u%BladePtLoads(K)%Force(2,NodeNum) * CoordSys%te2(K,J,:) ) / p%DRNodes(J) + + RtHSdat%MMAero(:,K,J) = CROSS_PRODUCT( RtHSdat%rSAerCen(:,J,K), RtHSdat%FSAero(:,K,J) )& + + u%BladePtLoads(K)%Moment(3,NodeNum)/p%DRNodes(J) * CoordSys%te3(K,J,:) + else + ! ... (不同的坐标系转换和分量赋值) + end if + ``` + 这部分代码从输入 `u%BladePtLoads` 中获取每个叶片元素的气动力和气动扭矩。 + * `u%BladePtLoads(K)%Force` 和 `u%BladePtLoads(K)%Moment` 是由外部气动模块（如 AeroDyn）计算并提供的物理力/力矩。 + * `p%DRNodes(J)` 是叶片元素 $J$ 的长度，用于将总力/力矩转换为单位展向长度上的力/力矩（对应 [[Blade]] 公式 (24) 中的 $F_{AeroB1}^{S1}(r)$ 和 $M_{AeroB1}^{M1}(r)$）。 + * `CROSS_PRODUCT( RtHSdat%rSAerCen(:,J,K), RtHSdat%FSAero(:,K,J) )` 计算了气动力在气动中心产生的力矩，这与 [[Blade]] 公式 (24) 中对 $M_{AeroB1}^{M1}(r)$ 的描述一致，即可以包含由气动偏移引起的力矩。 + +#### 2. 叶片根部偏力与偏力矩的计算 (`PFrcS0B`, `PMomH0B`) + +这部分计算的是叶片根部（点 S(0)）的**偏力 (Partial Force)** 和**偏力矩 (Partial Moment)**。这些“偏”量是惯性力对广义加速度的系数，用于构建广义质量矩阵。 + +* **理论依据:** [[Blade]] 中的公式 (1) 和 (4) 描述了广义惯性力，其中加速度项被分解为与 $\ddot{q}_i$ 相关的部分（用于质量矩阵）和与 $\dot{q}_i$ 相关的部分（用于广义力向量）。这里的偏力/偏力矩与广义质量矩阵的构建直接相关。 + +* **代码实现:** + * **叶尖质量的贡献:** + ```fortran + TmpVec1 = -p%TipMass(K)*RtHSdat%PLinVelES(K,p%TipNode,p%DOFs%PSE(K,I),0,:) + RtHSdat%PFrcS0B(:,K,p%DOFs%PSE(K,I)) = TmpVec1 + RtHSdat%PMomH0B(:,K,p%DOFs%PSE(K,I)) = CROSS_PRODUCT( RtHSdat%rS0S(:,K,p%TipNode), TmpVec1 ) + ``` + * `p%TipMass(K)` 是叶尖质量。 + * `RtHSdat%PLinVelES(...)` 是叶尖点 S(BldFlexL) 相对于广义坐标 $q_I$ 的**偏速度** ${}^{E}\boldsymbol{v}_{I}^{S1}(BldFlexL)$。 + * `TmpVec1` 实际上是叶尖质量乘以偏速度，这构成了广义质量矩阵的组成部分。负号表示惯性力的方向。 + * `CROSS_PRODUCT( RtHSdat%rS0S(:,K,p%TipNode), TmpVec1 )` 将这个惯性力转换为关于叶片根部（点 S(0)）的惯性力矩。`RtHSdat%rS0S` 是从叶片根部到叶尖的位置向量。 + * **分布式叶片元素的贡献:** + ```fortran + TmpVec1 = -p%BElmntMass(J,K)*RtHSdat%PLinVelES(K,J,p%DOFs%PSE(K,I),0,:) + RtHSdat%PFrcS0B(:,K,p%DOFs%PSE(K,I)) = RtHSdat%PFrcS0B(:,K,p%DOFs%PSE(K,I)) + TmpVec1 + RtHSdat%PMomH0B(:,K,p%DOFs%PSE(K,I)) = RtHSdat%PMomH0B(:,K,p%DOFs%PSE(K,I)) + & + CROSS_PRODUCT( RtHSdat%rS0S(:,K,J), TmpVec1 ) + ``` + 这部分代码以类似的方式累加了每个叶片元素质量对叶片根部偏力/偏力矩的贡献。`p%BElmntMass(J,K)` 是叶片元素 $J$ 的质量，`RtHSdat%PLinVelES(...)` 是叶片元素 $J$ 处点 $S1$ 相对于广义坐标 $q_I$ 的偏速度 ${}^{E}\boldsymbol{v}_{I}^{S1}(r)$。 + +#### 3. 叶片根部总力与总力矩的计算 (`FrcS0Bt`, `MomH0Bt`) + +这部分计算的是叶片根部（点 S(0)）的**总力 (Total Force)** 和**总力矩 (Total Moment)**。这些“总”量包含了所有主动力（气动力、重力）以及速度相关的惯性力（科里奥利力、离心力），它们共同构成了广义力向量。 + +* **理论依据:** [[Blade]] 中的公式 (6) 描述了速度相关的惯性力，公式 (23) 描述了广义重力，公式 (26) 描述了广义气动力。这些都是广义力向量的组成部分。 + +* **代码实现:** + * **叶尖贡献:** + ```fortran + TmpVec1 = RtHSdat%FSTipDrag(:,K) - p%TipMass(K)*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccESt(:,K,p%TipNode) ) + RtHSdat%FrcS0Bt(:,K) = TmpVec1 + RtHSdat%MomH0Bt(:,K) = CROSS_PRODUCT( RtHSdat%rS0S(:,K,p%TipNode), TmpVec1 ) + ``` + * `RtHSdat%FSTipDrag(:,K)`: 叶尖气动阻力。 + * `p%TipMass(K)*p%Gravity*CoordSys%z2`: 叶尖质量的重力。这对应于 [[Blade]] 公式 (21) 中的叶尖重力项。 + * `p%TipMass(K)*RtHSdat%LinAccESt(:,K,p%TipNode)`: 叶尖质量乘以其绝对加速度。`RtHSdat%LinAccESt` 包含了科里奥利加速度和离心加速度等速度相关的项。当这些项被投影到偏速度方向时，就构成了 [[Blade]] 公式 (6) 中的速度相关惯性力。 + * `CROSS_PRODUCT(...)` 将这些力转换为关于叶片根部的力矩。 + * **分布式叶片元素贡献:** + ```fortran + TmpVec1 = RtHSdat%FSAero(:,K,J)*p%DRNodes(J) - p%BElmntMass(J,K)*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccESt(:,K,J) ) + TmpVec2 = CROSS_PRODUCT( RtHSdat%rS0S(:,K,J), TmpVec1 ) + TmpVec3 = RtHSdat%MMAero(:,K,J)*p%DRNodes(J) + + RtHSdat%FrcS0Bt(:,K) = RtHSdat%FrcS0Bt(:,K) + TmpVec1 + RtHSdat%MomH0Bt(:,K) = RtHSdat%MomH0Bt(:,K) + TmpVec2 + TmpVec3 + ``` + * `RtHSdat%FSAero(:,K,J)*p%DRNodes(J)`: 叶片元素 $J$ 上的气动力。这对应于 [[Blade]] 公式 (24) 中的分布式气动力项。 + * `p%BElmntMass(J,K)*p%Gravity*CoordSys%z2`: 叶片元素 $J$ 的重力。这对应于 [[Blade]] 公式 (21) 中的分布式重力项。 + * `p%BElmntMass(J,K)*RtHSdat%LinAccESt(:,K,J)`: 叶片元素 $J$ 质量乘以其绝对加速度。同样，这包含了速度相关的惯性力。 + * `TmpVec3 = RtHSdat%MMAero(:,K,J)*p%DRNodes(J)`: 叶片元素 $J$ 上的气动扭矩。这对应于 [[Blade]] 公式 (24) 中的分布式气动扭矩项。 + * `CROSS_PRODUCT(...)` 将力转换为关于叶片根部的力矩。 + +### 与广义力的联系 + +`CalculateForcesMoments` 子程序计算并存储在 `RtHSdat` 结构体中的这些物理力（`FSAero`, `FSTipDrag`, `FrcS0Bt`, `PFrcS0B` 等）和力矩（`MMAero`, `MomH0Bt`, `PMomH0B` 等）是构建系统动力学方程的**中间结果**。 + +* **`PFrcS0B` 和 `PMomH0B`** (Partial Forces/Moments): 这些是与广义加速度 $\ddot{q}$ 相关的系数，它们在 `FillAugMat` 中与偏速度/偏角速度进行点积，用于填充**广义质量矩阵** $[C(q,t)]$。 +* **`FrcS0Bt` 和 `MomH0Bt`** (Total Forces/Moments): 这些是包含了所有主动力（气动力、重力、弹性力、阻尼力）以及速度相关的惯性力（科里奥利力、离心力）的总和。它们在 `FillAugMat` 中与偏速度/偏角速度进行点积，用于填充**广义力向量** $\{-f(\dot{q},q,t)\}$。 + +简而言之，`CalculateForcesMoments` 负责将外部输入（如气动载荷）和系统内部状态（如位置、速度、加速度）转化为作用在叶片上的具体物理力和力矩，并进一步分解为“偏”量和“总”量，为后续在 `FillAugMat` 中构建完整的广义动力学方程做准备。它将 [[Blade]] 笔记中描述的各种力（惯性、重力、气动）的物理表达式，通过具体的数值计算和向量运算实现出来。 +[Timestamp: 2026/01/07 14:08:00] + +**user**: [[Blade]][[CalculateAngularPosVelPAcc]]结合理论和代码，解释叶片部分是如何计算的 +[Timestamp: 2026/01/07 14:23:01] + +**ai**: `CalculateAngularPosVelPAcc` 子程序是风力机动力学模型中一个非常基础且关键的模块，它负责计算系统中各个刚体和柔性体元素的**角位置、角速度、偏角速度及其一阶导数（用于角加速度）**。这些运动学量是构建广义惯性力（特别是科里奥利力和离心力项）以及将物理力/力矩投影到广义坐标系中的基础。 + +本解释将聚焦于 `CalculateAngularPosVelPAcc` 中与叶片相关的计算部分，并将其与 [[Blade]] 笔记中的理论公式进行对照。 + +### 凯恩方程中的运动学基础 + +在凯恩方程中，任何点的速度 $\boldsymbol{v}$ 和任何刚体的角速度 $\boldsymbol{\omega}$ 都可以表示为广义速度 $\dot{q}_i$ 的线性组合： + +$\boldsymbol{v} = \sum_i \boldsymbol{v}_i \dot{q}_i + \boldsymbol{v}_t$ +$\boldsymbol{\omega} = \sum_i \boldsymbol{\omega}_i \dot{q}_i + \boldsymbol{\omega}_t$ + +其中： +* $\boldsymbol{v}_i$ 和 $\boldsymbol{\omega}_i$ 分别是相对于广义速度 $\dot{q}_i$ 的**偏速度 (Partial Velocity)** 和**偏角速度 (Partial Angular Velocity)**。 +* $\boldsymbol{v}_t$ 和 $\boldsymbol{\omega}_t$ 是与广义速度无关的项（通常是由于参考系运动或时间显式依赖）。 + +加速度 $\boldsymbol{a}$ 和角加速度 $\boldsymbol{\alpha}$ 则涉及偏速度/偏角速度对时间的导数： + +$\boldsymbol{a} = \sum_i \boldsymbol{v}_i \ddot{q}_i + \sum_i \frac{d\boldsymbol{v}_i}{dt} \dot{q}_i + \frac{d\boldsymbol{v}_t}{dt}$ +$\boldsymbol{\alpha} = \sum_i \boldsymbol{\omega}_i \ddot{q}_i + \sum_i \frac{d\boldsymbol{\omega}_i}{dt} \dot{q}_i + \frac{d\boldsymbol{\omega}_t}{dt}$ + +其中，$\frac{d\boldsymbol{v}_i}{dt}$ 和 $\frac{d\boldsymbol{\omega}_i}{dt}$ 是偏速度/偏角速度的**一阶导数**，它们在惯性力中产生了科里奥利力和离心力项。 + +`CalculateAngularPosVelPAcc` 子程序正是计算这些偏速度、偏角速度及其导数，以及总角速度和总角加速度中与 $\ddot{q}$ 无关的部分。 + +### `CalculateAngularPosVelPAcc` 子程序中的叶片部分计算 + +该子程序首先计算了从平台到轮毂的各个刚体的角运动学量，然后进入叶片循环 (`DO K = 1,p%NumBl`) 来计算叶片的柔性运动。 + +#### 1. 轮毂的角运动学量 (`RtHSdat%AngVelEH`, `RtHSdat%PAngVelEH`, `RtHSdat%AngAccEHt`) + +在计算叶片运动之前，需要先确定叶片根部所连接的轮毂（Hub，体 H）在惯性系 E 中的角运动学量。 + +* **`RtHSdat%PAngVelEH(:,0,:)` (轮毂的偏角速度):** + ```fortran + RtHSdat%PAngVelEH( :,0,:)= RtHSdat%PAngVelEL(:,0,:) + RtHSdat%AngVelEH = RtHSdat%AngVelEL + IF ( p%NumBl == 2 ) THEN ! 2-blader + RtHSdat%PAngVelEH(DOF_Teet,0,:)= CoordSys%f2 + RtHSdat%AngVelEH = RtHSdat%AngVelEH + x%QDT(DOF_Teet)*RtHSdat%PAngVelEH(DOF_Teet,0,:) + ENDIF + ``` + 这部分代码计算了轮毂相对于惯性系 E 的偏角速度。它基于低速轴 (LSS，体 L) 的偏角速度 `RtHSdat%PAngVelEL`，并考虑了可能的**摆振 (Teeter)** 自由度（如果存在）。`CoordSys%f2` 是摆振轴的单位向量。 +* **`RtHSdat%AngVelEH` (轮毂的总角速度):** + 通过将轮毂的偏角速度与对应的广义速度 `x%QDT(DOF)` 相乘并求和，得到轮毂在惯性系 E 中的总角速度。 +* **`RtHSdat%PAngVelEH(:,1,:)` (轮毂偏角速度的一阶导数):** + ```fortran + RtHSdat%PAngVelEH( :,1,:) = RtHSdat%PAngVelEL(:,1,:) + IF ( p%NumBl == 2 ) THEN ! 2-blader + RtHSdat%PAngVelEH(DOF_Teet,1,:) = CROSS_PRODUCT( RtHSdat%AngVelEH, RtHSdat%PAngVelEH(DOF_Teet,0,:) ) + ... + ENDIF + ``` + 这部分计算了轮毂偏角速度的一阶导数。对于旋转的参考系，偏角速度的导数通常涉及当前角速度与偏角速度的叉积（例如，科里奥利效应）。 +* **`RtHSdat%AngAccEHt` (轮毂角加速度中与 $\ddot{q}$ 无关的部分):** + ```fortran + RtHSdat%AngAccEHt = AngAccELt + IF ( p%NumBl == 2 ) THEN ! 2-blader + ... + RtHSdat%AngAccEHt = RtHSdat%AngAccEHt + x%QDT(DOF_Teet)*RtHSdat%PAngVelEH(DOF_Teet,1,:) + ENDIF + ``` + 这部分计算了轮毂在惯性系 E 中的角加速度中与广义加速度 $\ddot{q}$ 无关的项。它包含了由于广义速度 $\dot{q}$ 引起的科里奥利和离心效应。 + +#### 2. 叶片柔性点的角运动学量 (`RtHSdat%PAngVelEM`, `RtHSdat%AngVelEM`, `RtHSdat%AngAccEKt`) + +接下来，子程序遍历每个叶片（`K`）上的每个节点/元素（`J`），计算这些柔性点的角运动学量。 + +* **理论依据:** [[Blade]] 中的公式 (4) 描述了叶片上点 $S1$ 的绝对加速度，其中包含了与偏速度一阶导数相关的项。这些偏速度一阶导数正是由角运动学量（特别是偏角速度的导数）推导而来。公式 (24) 中的广义气动力也需要叶片元素的偏角速度 $\boldsymbol{\omega}_r^{M1}(r)$。 + +* **代码实现:** + * **叶片元素偏角速度 (`RtHSdat%PAngVelEM(K,J,:,0,:)`):** + ```fortran + RtHSdat%PAngVelEM(K,J, :,0,:) = RtHSdat%PAngVelEH(:,0,:) + RtHSdat%PAngVelEM(K,J,DOF_BF(K,1),0,:) = - p%TwistedSF(K,2,1,J,1)*CoordSys%j1(K,:) & + + p%TwistedSF(K,1,1,J,1)*CoordSys%j2(K,:) + RtHSdat%PAngVelEM(K,J,DOF_BF(K,2),0,:) = - p%TwistedSF(K,2,2,J,1)*CoordSys%j1(K,:) & + + p%TwistedSF(K,1,2,J,1)*CoordSys%j2(K,:) + RtHSdat%PAngVelEM(K,J,DOF_BE(K,1),0,:) = - p%TwistedSF(K,2,3,J,1)*CoordSys%j1(K,:) & + + p%TwistedSF(K,1,3,J,1)*CoordSys%j2(K,:) + ``` + 这部分代码计算了叶片元素 $J$ 相对于惯性系 E 的偏角速度。它首先继承了轮毂的偏角速度 `RtHSdat%PAngVelEH`，然后添加了由叶片柔性模态（挥舞 `DOF_BF` 和摆振 `DOF_BE`）引起的额外偏角速度。 + * `p%TwistedSF` 是**扭转形状函数 (Twisted Shape Function)**，它描述了叶片在柔性变形时，其横截面的扭转和弯曲如何影响局部坐标系的方向。这些形状函数与 `CoordSys%j1(K,:)` 和 `CoordSys%j2(K,:)`（叶片局部坐标系的单位向量）结合，定义了柔性模态引起的偏角速度分量。 + * 这些偏角速度 `RtHSdat%PAngVelEM` 直接对应于 [[Blade]] 公式 (24) 中的 ${}^E\boldsymbol{\omega}_r^{M1}(r)$。 + * **叶片元素总角速度 (`RtHSdat%AngVelEM`):** + ```fortran + AngVelHM = x%QDT(DOF_BF(K,1))*RtHSdat%PAngVelEM(K,J,DOF_BF(K,1),0,:) & + + x%QDT(DOF_BF(K,2))*RtHSdat%PAngVelEM(K,J,DOF_BF(K,2),0,:) & + + x%QDT(DOF_BE(K,1))*RtHSdat%PAngVelEM(K,J,DOF_BE(K,1),0,:) + RtHSdat%AngVelEM(:,J,K ) = RtHSdat%AngVelEH + AngVelHM + ``` + 叶片元素的总角速度是轮毂的总角速度 `RtHSdat%AngVelEH` 加上由叶片柔性模态引起的相对角速度 `AngVelHM`。 + * **叶片元素偏角速度的一阶导数 (`RtHSdat%PAngVelEM(K,J,:,1,:)`):** + 代码中这部分被注释掉了： + ```fortran + ! RtHSdat%PAngVelEM(K,J, :,1,:) = RtHSdat%PAngVelEH(:,1,:) + ! RtHSdat%PAngVelEM(K,J,DOF_BF(K,1),1,:) = CROSS_PRODUCT( RtHSdat%AngVelEH, PAngVelEM(K,J,DOF_BF(K,1),0,:) ) + ! RtHSdat%PAngVelEM(K,J,DOF_BF(K,2),1,:) = CROSS_PRODUCT( RtHSdat%AngVelEH, PAngVelEM(K,J,DOF_BF(K,2),0,:) ) + ! RtHSdat%PAngVelEM(K,J,DOF_BE(K,1),1,:) = CROSS_PRODUCT( RtHSdat%AngVelEH, PAngVelEM(K,J,DOF_BE(K,1),0,:) ) + ``` + 注释说明这些计算目前未被使用。如果需要考虑叶片元素的惯性力矩（例如，当叶片元素具有显著的转动惯量时），则需要计算这些偏角速度的一阶导数。这些导数在凯恩方程中用于计算角加速度中与广义速度相关的项。 +* **`RtHSdat%AngAccEKt` (叶片元素角加速度中与 $\ddot{q}$ 无关的部分):** + ```fortran + RtHSdat%AngAccEKt(:,J ,K) = RtHSdat%AngAccEHt + x%QDT(DOF_BF(K,1))*RtHSdat%PAngVelEM(K,J,DOF_BF(K,1),1,:) & + + x%QDT(DOF_BF(K,2))*RtHSdat%PAngVelEM(K,J,DOF_BF(K,2),1,:) & + + x%QDT(DOF_BE(K,1))*RtHSdat%PAngVelEM(K,J,DOF_BE(K,1),1,:) + ``` + 这部分计算了叶片元素在惯性系 E 中的角加速度中与广义加速度 $\ddot{q}$ 无关的项。它结合了轮毂的相应项 `RtHSdat%AngAccEHt` 和由叶片柔性模态引起的、与广义速度 `x%QDT(DOF)` 及其偏角速度一阶导数 `RtHSdat%PAngVelEM(K,J,DOF,1,:)` 相关的项。这些项在 [[Blade]] 公式 (4) 中，当计算绝对加速度时，会通过叉积运算产生科里奥利和离心加速度。 + +### 总结 + +`CalculateAngularPosVelPAcc` 子程序是整个风力机动力学模型中运动学计算的基石。它通过层层递进的方式，从基础的平台运动开始，逐步计算出轮毂以及每个叶片柔性点的角位置、角速度、偏角速度及其一阶导数。 + +这些计算结果对于以下方面至关重要： + +1. **广义惯性力的计算:** [[Blade]] 公式 (4) 和 (6) 中的速度相关惯性力项（科里奥利力和离心力）需要偏速度和偏角速度的一阶导数。 +2. **广义气动力的计算:** [[Blade]] 公式 (24) 中的气动扭矩项需要叶片元素的偏角速度。 +3. **坐标系转换:** 准确的角位置和角速度是进行各种力（如气动力、重力）在不同坐标系之间转换的基础。 + +通过这些详细的运动学计算，`CalculateAngularPosVelPAcc` 为后续的力和力矩计算（如 `CalculateForcesMoments`）和最终的运动方程构建（如 `FillAugMat`）提供了所有必要的运动学输入。 +[Timestamp: 2026/01/07 14:23:28] + +**user**: [[Blade]][[CalculatePositions]]结合理论和代码，解释代码的作用 +[Timestamp: 2026/01/07 14:29:33] + +**ai**: `CalculatePositions` 子程序在风力机动力学模型中负责计算系统中各个关键点和分布式元素（包括叶片）的**位置向量**。这些位置信息是进行运动学分析、计算力臂、以及最终构建广义力方程的基础。 + +本解释将聚焦于 `CalculatePositions` 中与叶片相关的计算部分，并将其与 [[Blade]] 笔记中的理论公式进行对照。 + +### 凯恩方程中的位置学基础 + +在凯恩方程中，系统的构型由广义坐标 $q_i$ 描述。系统中任意一点的绝对位置向量 $\boldsymbol{r}$ 可以表示为广义坐标的函数： + +$\boldsymbol{r} = \boldsymbol{r}(q_1, q_2, \ldots, q_n, t)$ + +这些位置向量是计算偏速度、偏角速度以及后续的力和力矩（特别是力矩的力臂）的基础。 + +### `CalculatePositions` 子程序中的叶片部分计算 + +`CalculatePositions` 子程序首先计算了从平台到轮毂的各个刚体之间的位置向量，然后进入叶片循环 (`DO K = 1,p%NumBl`) 来计算叶片柔性点的位移。 + +#### 1. 叶片根部到叶片柔性点的相对位置 (`RtHSdat%rS0S`) + +这部分代码计算了从叶片根部（点 S(0)）到叶片上任意柔性点（点 S(r)）的相对位置向量。这个相对位置包含了叶片的柔性变形。 + +* **理论依据:** 叶片的柔性变形通常通过模态叠加法来描述。叶片上任意一点的位移可以表示为各阶模态形状函数与对应广义坐标的乘积之和。虽然 [[Blade]] 笔记中没有直接给出位置向量的公式，但广义弹性力 (公式 7, 8) 和广义质量 (公式 16, 17) 都依赖于模态形状函数 $\phi(r)$，这些形状函数正是描述叶片柔性变形的基础。 + +* **代码实现:** + * **叶尖位置 (`RtHSdat%rS0S(:,K,p%TipNode)`):** + ```fortran + RtHSdat%rS0S(:,K,p%TipNode) = ( p%TwistedSF(K,1,1,p%TipNode,0)*x%QT( DOF_BF(K,1) ) & + + p%TwistedSF(K,1,2,p%TipNode,0)*x%QT( DOF_BF(K,2) ) & + + p%TwistedSF(K,1,3,p%TipNode,0)*x%QT( DOF_BE(K,1) ) )*CoordSys%j1(K,:) & + + ( p%TwistedSF(K,2,1,p%TipNode,0)*x%QT( DOF_BF(K,1) ) & + + p%TwistedSF(K,2,2,p%TipNode,0)*x%QT( DOF_BF(K,2) ) & + + p%TwistedSF(K,2,3,p%TipNode,0)*x%QT( DOF_BE(K,1) ) )*CoordSys%j2(K,:) & + + ( p%BldFlexL - 0.5* & + ( p%AxRedBld(K,1,1,p%TipNode)*x%QT( DOF_BF(K,1) )*x%QT( DOF_BF(K,1) ) & + + p%AxRedBld(K,2,2,p%TipNode)*x%QT( DOF_BF(K,2) )*x%QT( DOF_BF(K,2) ) & + + p%AxRedBld(K,3,3,p%TipNode)*x%QT( DOF_BE(K,1) )*x%QT( DOF_BE(K,1) ) & + + 2.*p%AxRedBld(K,1,2,p%TipNode)*x%QT( DOF_BF(K,1) )*x%QT( DOF_BF(K,2) ) & + + 2.*p%AxRedBld(K,2,3,p%TipNode)*x%QT( DOF_BF(K,2) )*x%QT( DOF_BE(K,1) ) & + + 2.*p%AxRedBld(K,1,3,p%TipNode)*x%QT( DOF_BF(K,1) )*x%QT( DOF_BE(K,1) ) ) )*CoordSys%j3(K,:) + ``` + 这行代码计算了从叶片根部到叶尖的相对位置向量 `rS0S`。 + * `p%TwistedSF(K,i,j,Node,0)`: 这是**扭转形状函数 (Twisted Shape Function)**，它描述了叶片在柔性变形时，其在叶片局部坐标系 `j1`, `j2` 方向上的位移。`x%QT(DOF_BF(K,1))` 等是对应的广义坐标（模态幅值）。 + * `CoordSys%j1(K,:)`, `CoordSys%j2(K,:)`, `CoordSys%j3(K,:)`: 叶片 $K$ 的局部坐标系单位向量。 + * `p%BldFlexL`: 叶片柔性段的长度。 + * `p%AxRedBld(K,i,j,Node)`: **轴向缩短系数 (Axial Shortening Coefficient)**。由于叶片的弯曲变形会导致其轴向长度的微小缩短（二阶效应），这些系数用于修正沿叶片轴向 `j3` 方向的位移。这是一个非线性项，与广义坐标的平方项相关。 + * **物理意义:** 整个表达式通过模态叠加和轴向缩短效应，精确地描述了叶片在柔性变形下的三维位移。 + * **分布式叶片元素位置 (`RtHSdat%rS0S(:,K,J)`):** + 在 `DO J = 1,p%BldNodes` 循环中，以类似的方式计算了从叶片根部到每个叶片元素 $J$ 的相对位置向量。 + +#### 2. 轮毂旋转中心到叶片柔性点的相对位置 (`RtHSdat%rQS`) + +这部分计算了从轮毂旋转中心（点 Q）到叶片上任意柔性点（点 S(r)）的相对位置向量。 + +* **代码实现:** + ```fortran + RtHSdat%rQS (:,K,p%TipNode) = RtHSdat%rS0S(:,K,p%TipNode) + p%HubRad*CoordSys%j3(K,:) + ``` + 这行代码将叶片根部到叶尖的相对位置 `rS0S` 加上从轮毂旋转中心到叶片根部的径向距离 `p%HubRad*CoordSys%j3(K,:)`，从而得到从轮毂旋转中心到叶尖的相对位置 `rQS`。 + * `p%HubRad`: 轮毂半径，即从轮毂旋转中心到叶片根部的径向距离。 + +#### 3. 惯性系原点到叶片柔性点的绝对位置 (`RtHSdat%rS`) + +这部分计算了从惯性系原点（通常是地面）到叶片上任意柔性点（点 S(r)）的绝对位置向量。 + +* **代码实现:** + ```fortran + RtHSdat%rS (:,K,p%TipNode) = RtHSdat%rQS (:,K,p%TipNode) + RtHSdat%rQ + ``` + 这行代码将从轮毂旋转中心到叶尖的相对位置 `rQS` 加上从惯性系原点到轮毂旋转中心的绝对位置 `RtHSdat%rQ`，从而得到从惯性系原点到叶尖的绝对位置 `rS`。 + * `RtHSdat%rQ`: 这是在叶片循环之前计算的，表示从惯性系原点到轮毂旋转中心（点 Q）的绝对位置向量。 + +#### 4. 其他关键位置向量 + +除了叶片柔性点的位置，子程序还计算了其他与叶片相关的关键位置向量： + +* **`RtHSdat%rPS0(:,K)` (从摆振销到叶片根部的位置):** + ```fortran + RtHSdat%rPS0(:,K) = RtHSdat%rPQ + p%HubRad*CoordSys%j3(K,:) + ``` + 这表示从摆振销（Teeter Pin，点 P）到叶片根部（点 S(0)）的位置向量。`RtHSdat%rPQ` 是从摆振销到轮毂旋转中心的位置向量。这个向量在计算叶片根部力矩时非常重要。 + +### 总结 + +`CalculatePositions` 子程序通过以下步骤，为叶片部分的动力学分析提供了全面的位置信息： + +1. **定义基准位置:** 首先计算了风力机各个主要刚体（平台、塔架、机舱、轮毂）之间的相对和绝对位置向量。 +2. **处理叶片柔性:** 针对每个叶片，利用**扭转形状函数**和**轴向缩短系数**，结合广义坐标 `x%QT`，精确地计算了叶片上每个柔性点相对于叶片根部、轮毂旋转中心以及惯性系原点的三维位置。 +3. **提供力臂信息:** 这些位置向量是计算力矩的关键力臂。例如，在 `CalculateForcesMoments` 中计算叶片根部力矩时，就需要用到 `RtHSdat%rS0S` 和 `RtHSdat%rPS0` 等作为力臂。 +4. **支持运动学计算:** 准确的位置信息是计算速度和加速度（包括偏速度和偏角速度）的基础，这些运动学量在 `CalculateAngularPosVelPAcc` 和 `CalculateLinearVelPAcc` 中被进一步计算。 + +通过这些详细的位置计算，`CalculatePositions` 确保了风力机叶片在复杂三维运动和柔性变形下的几何构型被准确地捕捉，为后续的动力学方程构建提供了坚实的基础。 +[Timestamp: 2026/01/07 14:29:57] \ No newline at end of file diff --git a/书籍/力学书籍/OpenFast/FASTKinetics/auto/Blade.md b/书籍/力学书籍/OpenFast/FASTKinetics/auto/Blade.md index 409b47c..15ce885 100644 --- a/书籍/力学书籍/OpenFast/FASTKinetics/auto/Blade.md +++ b/书籍/力学书籍/OpenFast/FASTKinetics/auto/Blade.md @@ -1,51 +1,51 @@ The distributed properties of blade 1 bring about generalized inertia forces and generalized active forces associated with blade elasticity, blade damping, blade weight, and blade aerodynamics. - +叶片1的分布式特性导致与叶片弹性、叶片阻尼、叶片重量及叶片气动力相关的广义惯性力和广义主动力。 $$ -F_{r}^{*}\big|_{Bl} = -\int_{0}^{Bl \, dFlexL} \mu^{Bl}(r) \, {}^E \, \boldsymbol{v}_{r}^{SI}(r) \cdot E \, \boldsymbol{a}^{SI}(r) \, dr - m^{BldTip} \, {}^E \, \boldsymbol{v}_{r}^{SI}(Bl \, dFlexL) \cdot {}^E \, \boldsymbol{a}^{SI}(Bl \, dFlexL) \quad (r = 1, 2, \ldots, 22) +F_{r}^{*}\big|_{B1} = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, {}^E \, \boldsymbol{v}_{r}^{S1}(r) \cdot {}^E \, \boldsymbol{a}^{S1}(r) \, dr - m^{BldTip} \, {}^E \, \boldsymbol{v}_{r}^{S1}(BldFlexL) \cdot {}^E \, \boldsymbol{a}^{S1}(BldFlexL) \quad (r = 1, 2, \ldots, 22) \tag{1} $$ 其中 $$ -\mu^{B1}(r) = AdjBlMs^{B1} \cdot BMassDen^{B1}(r) +\mu^{B1}(r) = AdjBlMs^{B1} \cdot BMassDen^{B1}(r) \tag{2} $$ $$ -m^{B1Tip} = TipMass(1) +m^{B1Tip} = TipMass(1) \tag{3} $$ or, $$ -F_{r}^{*}\big|_{B1} = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, {}^E\boldsymbol{v}_{r}^{SI}(r) \cdot \left\{ +F_{r}^{*}\big|_{B1} = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, {}^E\boldsymbol{v}_{r}^{S1}(r) \cdot \left\{ \begin{array}{l} -\left( \sum_{i=1}^{14} {}^{E}\boldsymbol{v}_{i}^{SI}(r) \ddot{q}_{i} \right) + \left( \sum_{i=16}^{18} {}^{E}\boldsymbol{v}_{i}^{SI}(r) \ddot{q}_{i} \right) + {}^{E}\boldsymbol{v}_{Teet}^{SI}(r) \ddot{q}_{Teet} \\ -+ \left[ \sum_{i=4}^{14} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{SI}(r) \right) \dot{q}_{i} \right] + \left[ \sum_{i=16}^{18} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{SI}(r) \right) \dot{q}_{i} \right] + \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{Teet}^{SI}(r) \right) \dot{q}_{Teet} +\left( \sum_{i=1}^{14} {}^{E}\boldsymbol{v}_{i}^{S1}(r) \ddot{q}_{i} \right) + \left( \sum_{i=16}^{18} {}^{E}\boldsymbol{v}_{i}^{S1}(r) \ddot{q}_{i} \right) + {}^{E}\boldsymbol{v}_{Teet}^{S1}(r) \ddot{q}_{Teet} \\ ++ \left[ \sum_{i=4}^{14} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{S1}(r) \right) \dot{q}_{i} \right] + \left[ \sum_{i=16}^{18} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{S1}(r) \right) \dot{q}_{i} \right] + \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{Teet}^{S1}(r) \right) \dot{q}_{Teet} \end{array} \right\} dr \\ -- m^{B1Tip} {}^{E}\boldsymbol{v}_{r}^{SI}(BldFlexL) \cdot \left\{ +- m^{B1Tip} {}^{E}\boldsymbol{v}_{r}^{S1}(BldFlexL) \cdot \left\{ \begin{array}{l} -\left( \sum_{i=1}^{14} {}^{E}\boldsymbol{v}_{i}^{SI}(BldFlexL) \ddot{q}_{i} \right) + \left( \sum_{i=16}^{18} {}^{E}\boldsymbol{v}_{i}^{SI}(BldFlexL) \ddot{q}_{i} \right) + {}^{E}\boldsymbol{v}_{Teet}^{SI}(BldFlexL) \ddot{q}_{Teet} \\ -+ \left[ \sum_{i=4}^{14} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{SI}(BldFlexL) \right) \dot{q}_{i} \right] + \left[ \sum_{i=16}^{18} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{SI}(BldFlexL) \right) \dot{q}_{i} \right] \\ -+ \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{Teet}^{SI}(BldFlexL) \right) \dot{q}_{Teet} +\left( \sum_{i=1}^{14} {}^{E}\boldsymbol{v}_{i}^{S1}(BldFlexL) \ddot{q}_{i} \right) + \left( \sum_{i=16}^{18} {}^{E}\boldsymbol{v}_{i}^{S1}(BldFlexL) \ddot{q}_{i} \right) + {}^{E}\boldsymbol{v}_{Teet}^{S1}(BldFlexL) \ddot{q}_{Teet} \\ ++ \left[ \sum_{i=4}^{14} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{S1}(BldFlexL) \right) \dot{q}_{i} \right] + \left[ \sum_{i=16}^{18} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{S1}(BldFlexL) \right) \dot{q}_{i} \right] \\ ++ \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{Teet}^{S1}(BldFlexL) \right) \dot{q}_{Teet} \end{array} -\right\} \quad (r = 1, 2, \ldots, 14; 16, 17, 18; Teet) +\right\} \quad (r = 1, 2, \ldots, 14; 16, 17, 18; Teet) \tag{4} $$ 因此， $$ -\left[C(q,t)\right]_{B1} (\text{Row,Col}) = \int_{0}^{BldFlexL} \mu^{B1}(r) \, {}^{E}\boldsymbol{v}_{\text{Row}}^{S1}(r) \cdot {}^{E}\boldsymbol{v}_{\text{Col}}^{SI}(r) \, dr + m^{B1Tip} {}^{E}\boldsymbol{v}_{\text{Row}}^{SI}(BldFlexL) \cdot {}^{E}\boldsymbol{v}_{\text{Col}}^{SI}(BldFlexL) \quad (\text{Row, Col} = 1,2,\ldots,14;16,17,18;22) +\left[C(q,t)\right]_{B1} (\text{Row,Col}) = \int_{0}^{BldFlexL} \mu^{B1}(r) \, {}^{E}\boldsymbol{v}_{\text{Row}}^{S1}(r) \cdot {}^{E}\boldsymbol{v}_{\text{Col}}^{S1}(r) \, dr + m^{B1Tip} {}^{E}\boldsymbol{v}_{\text{Row}}^{S1}(BldFlexL) \cdot {}^{E}\boldsymbol{v}_{\text{Col}}^{S1}(BldFlexL) \quad (\text{Row, Col} = 1,2,\ldots,14;16,17,18;22) \tag{5} $$ $$ -\left\{-f(\dot{q},q,t)\right\}_{B1} (\text{Row}) = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, {}^{E}\boldsymbol{v}_{\text{Row}}^{SI}(r) \cdot \left\{ \left[ \sum_{i=4}^{14} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{SI}(r) \right) \dot{q}_{i} \right] + \left[ \sum_{i=16}^{18} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{SI}(r) \right) \dot{q}_{i} \right] + \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{Teet}^{SI}(r) \right) \dot{q}_{Teet} \right\} dr \\ -- m^{B1Tip} {}^{E}\boldsymbol{v}_{\text{Row}}^{SI}(BldFlexL) \cdot \left\{ \left[ \sum_{i=4}^{14} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{SI}(BldFlexL) \right) \dot{q}_{i} \right] + \left[ \sum_{i=16}^{18} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{SI}(BldFlexL) \right) \dot{q}_{i} \right] + \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{Teet}^{SI}(BldFlexL) \right) \dot{q}_{Teet} \right\} \quad (\text{Row} = 1,2,\ldots,14;16,17,18;22) +\left\{-f(\dot{q},q,t)\right\}_{B1} (\text{Row}) = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, {}^{E}\boldsymbol{v}_{\text{Row}}^{S1}(r) \cdot \left\{ \left[ \sum_{i=4}^{14} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{S1}(r) \right) \dot{q}_{i} \right] + \left[ \sum_{i=16}^{18} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{S1}(r) \right) \dot{q}_{i} \right] + \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{Teet}^{S1}(r) \right) \dot{q}_{Teet} \right\} dr \\ +- m^{B1Tip} {}^{E}\boldsymbol{v}_{\text{Row}}^{S1}(BldFlexL) \cdot \left\{ \left[ \sum_{i=4}^{14} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{S1}(BldFlexL) \right) \dot{q}_{i} \right] + \left[ \sum_{i=16}^{18} \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{i}^{S1}(BldFlexL) \right) \dot{q}_{i} \right] + \frac{d}{dt} \left( {}^{E}\boldsymbol{v}_{Teet}^{S1}(BldFlexL) \right) \dot{q}_{Teet} \right\} \quad (\text{Row} = 1,2,\ldots,14;16,17,18;22) \tag{6} $$ $$ -F_{r}\big|_{ElasticB1} = -\frac{\partial V'^{B1}}{\partial q_r} \quad (r = 1, 2, \ldots, 22) +F_{r}\big|_{ElasticB1} = -\frac{\partial V'^{B1}}{\partial q_r} \quad (r = 1, 2, \ldots, 22) \tag{7} $$ $$ \left. F_{r} \right|_{ElasticB1} = @@ -58,65 +58,65 @@ $$ & \text{for } r = B1F2 \\ 0 & \text{otherwise} -\end{cases} +\end{cases} \tag{8} $$ where $k_{\ i j}^{\,\prime B I F}$ and $\boldsymbol{k\,}_{I I}^{\prime B I E}$ are the generalized stiffnesses of blade 1 in the local flap and local edge directions respectively when centrifugal-stiffening effects are not included as follows: - +其中，$k_{ij}^{\,\prime BIF}$ 和 $\boldsymbol{k}_{II}^{\prime BIE}$ 分别为未考虑离心刚化效应时，叶片1在局部挥舞方向和局部摆振方向的广义刚度，具体如下： $$ -{k'}_{ij}^{B1F} = \sqrt{FlStTunr^{B1}(i) \, FlStTunr^{B1}(j)} \int_{0}^{BldFlexL} EI^{B1F}(r) \frac{d^2 \phi_i^{B1F}(r)}{dr^2} \frac{d^2 \phi_j^{B1F}(r)}{dr^2} dr \quad (i, j = 1, 2) +{k'}_{ij}^{B1F} = \sqrt{FlStTunr^{B1}(i) \, FlStTunr^{B1}(j)} \int_{0}^{BldFlexL} EI^{B1F}(r) \frac{d^2 \phi_i^{B1F}(r)}{dr^2} \frac{d^2 \phi_j^{B1F}(r)}{dr^2} dr \quad (i, j = 1, 2) \tag{9} $$ 其中， $$ -EI^{B1F}(r) = AdjFlSt^{B1} \cdot FlpStff^{B1}(r) +EI^{B1F}(r) = AdjFlSt^{B1} \cdot FlpStff^{B1}(r) \tag{10} $$ $$ -{k'}_{11}^{B1E} = \int_{0}^{BldFlexL} EI^{B1E}(r) \left[ \frac{d^2 \phi_{1}^{B1E}(r)}{dr^2} \right]^2 dr +{k'}_{11}^{B1E} = \int_{0}^{BldFlexL} EI^{B1E}(r) \left[ \frac{d^2 \phi_{1}^{B1E}(r)}{dr^2} \right]^2 dr \tag{11} $$ 其中， $$ -EI^{B1E}(r) = AdjEdSt^{B1} \cdot EdgStff^{B1}(r) +EI^{B1E}(r) = AdjEdSt^{B1} \cdot EdgStff^{B1}(r) \tag{12} $$ Similarly, when using the Rayleigh damping technique where the damping is assumed proportional to the stiffness, then - +同样，在使用瑞利阻尼技术（假设阻尼与刚度成正比）时，则 $$ \left. F_r \right|_{DampB1} = \begin{cases} --\frac{\zeta_1^{B1F} {k'}_{11}^{B1F}}{\pi f^{\prime,B 1 F}_1} \dot{q}_{B1F1} - \frac{\zeta_2^{B1F} {k'}_{12}^{B1F}}{\pi f^{\prime,B 1 F}_2} \dot{q}_{B1F2} & \text{for } r = B1F1 \\ --\frac{\zeta_1^{B1E} {k'}_{11}^{B1E}}{\pi f^{\prime,B 1 E}_1} \dot{q}_{B1E1} & \text{for } r = B1E1 \\ --\frac{\zeta_1^{B1F} {k'}_{21}^{B1F}}{\pi f^{\prime,B 1 F}_1} \dot{q}_{B1F1} - \frac{\zeta_2^{B1F} {k'}_{22}^{B1F}}{\pi f^{\prime,B 1 F}_2} \dot{q}_{B1F2} & \text{for } r = B1F2 \\ +-\frac{\zeta_1^{B1F} {k'}_{11}^{B1F}}{\pi f^{\prime B 1 F}_1} \dot{q}_{B1F1} - \frac{\zeta_2^{B1F} {k'}_{12}^{B1F}}{\pi f^{\prime B 1 F}_2} \dot{q}_{B1F2} & \text{for } r = B1F1 \\ +-\frac{\zeta_1^{B1E} {k'}_{11}^{B1E}}{\pi f^{\prime B 1 E}_1} \dot{q}_{B1E1} & \text{for } r = B1E1 \\ +-\frac{\zeta_1^{B1F} {k'}_{21}^{B1F}}{\pi f^{\prime B 1 F}_1} \dot{q}_{B1F1} - \frac{\zeta_2^{B1F} {k'}_{22}^{B1F}}{\pi f^{\prime B 1 F}_2} \dot{q}_{B1F2} & \text{for } r = B1F2 \\ 0 & \text{otherwise} -\end{cases} +\end{cases} \tag{13} $$ where $\zeta_{i}^{B I F}$ and $\zeta_{i}^{B I E}$ represent the structural damping ratio of blade 1 for the $i^{\mathrm{th}}$ mode in the local flap and edge directions, $B l d F l D m p^{B 1}(i)/100$ and $B l d E d D m p^{B 1}\left(i\right)/100$ respectively. Also, ${f^{\prime}}_{i}^{B I F}$ and $\boldsymbol{f}_{\ i}^{\prime B I E}$ represent the natural frequency of blade 1 for the $i^{\mathrm{th}}$ mode in the local flap and edge directions respectively without centrifugal-stiffening effects. That is, - +其中，$\zeta_{i}^{BIF}$ 和 $\zeta_{i}^{BIE}$ 分别代表叶片1在第 $i$ 阶模态下在**挥舞**方向和**摆振**方向的结构阻尼比，对应表达式为 $BldFlDmp^{B1}(i)/100$ 和 $BldEdDmp^{B1}(i)/100$。此外，${f^{\prime}}_{i}^{B I F}$ 和 $\boldsymbol{f}_{\ i}^{\prime B I E}$ 分别表示叶片1在第 $i$ 阶模态下，**不考虑离心刚度效应**时的**挥舞**方向和**摆振**方向的固有频率。也就是说， $$ -f_i^{\prime B1F} = \frac{1}{2\pi} \sqrt{\frac{k_{ii}^{\prime B1F}}{m_{ii}^{\prime B1F}}} +f_i^{\prime B1F} = \frac{1}{2\pi} \sqrt{\frac{k_{ii}^{\prime B1F}}{m_{ii}^{\prime B1F}}} \tag{14} $$ $$ -f_i^{\prime B1E} = \frac{1}{2\pi} \sqrt{\frac{k_{ii}^{\prime B1E}}{m_{ii}^{\prime B1E}}} +f_i^{\prime B1E} = \frac{1}{2\pi} \sqrt{\frac{k_{ii}^{\prime B1E}}{m_{ii}^{\prime B1E}}} \tag{15} $$ where $m_{\ i i}^{\prime B I F}$ and $m_{\ i i}^{\prime B I E}$ represent the generalized mass of blade 1 for the $i^{\mathrm{th}}$ mode in the local flap and edge directions respectively without centrifugal-stiffening and tip mass effects as follows: - +其中，$m_{ii}^{\prime BIF}$ 和 $m_{ii}^{\prime BIE}$ 分别表示在未考虑离心刚化和叶尖质量效应的情况下，第 $i$ 阶模态下叶片1在**挥舞方向**和**摆振方向**的广义质量，具体定义如下： $$ -m_{ij}^{\prime B1F} = \int_{0}^{BldFlexL} \mu^{B1}(r) \phi_i^{B1F}(r) \phi_j^{B1F}(r) \, dr \quad (i, j = 1, 2) +m_{ij}^{\prime B1F} = \int_{0}^{BldFlexL} \mu^{B1}(r) \phi_i^{B1F}(r) \phi_j^{B1F}(r) \, dr \quad (i, j = 1, 2) \tag{16} $$ $$ -m_{11}^{\prime B1E} = \int_{0}^{BldFlexL} \mu^{B1}(r) \left[ \phi_1^{B1E}(r) \right]^2 dr +m_{11}^{\prime B1E} = \int_{0}^{BldFlexL} \mu^{B1}(r) \left[ \phi_1^{B1E}(r) \right]^2 dr \tag{17} $$ Thus $$ -[C(q,t)]_{ElasticB1} = 0 \quad \text{and} \quad [C(q,t)]_{DampB1} = 0 +[C(q,t)]|_{ElasticB1} = 0 \quad \text{and} \quad [C(q,t)]|_{DampB1} = 0 \tag{18} $$ $$ @@ -131,7 +131,7 @@ $$ \vdots \\ \vdots \\ \vdots -\end{pmatrix} +\end{pmatrix} \tag{19} $$ $$ @@ -146,10 +146,30 @@ $$ \vdots \\ \vdots \\ \vdots -\end{pmatrix} +\end{pmatrix} \tag{20} $$$$ -F_{r}\big|_{GravB1} = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, \boldsymbol{g}^{E}\boldsymbol{v}_{r}^{SI}(r) \cdot \boldsymbol{{z}}_2 \, dr - m^{B1Tip} \boldsymbol{g}^{E}\boldsymbol{v}_{r}^{S1}(BldFlexL) \cdot \boldsymbol{{z}}_2 \quad (r = 3,4,\ldots,14;16,17,18;Teet) +F_{r}\big|_{GravB1} = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, \boldsymbol{g}^{E}\boldsymbol{v}_{r}^{S1}(r) \cdot \boldsymbol{{z}}_2 \, dr - m^{B1Tip} \boldsymbol{g}^{E}\boldsymbol{v}_{r}^{S1}(BldFlexL) \cdot \boldsymbol{{z}}_2 \quad (r = 3,4,\ldots,14;16,17,18;Teet)\tag{21} +$$ + +Thus, +$$ +[C(q,t)]|_{GravB1} = 0\tag{22} +$$ +$$ +\left\{ -f(\dot{q}, q, t) \right\}|_{GravB1} (\text{Row}) = -\int_{0}^{BldFlexL} \mu^{B1}(r) \, g \, {}^E\boldsymbol{v}_{Row}^{S1}(r) \cdot \boldsymbol{z}_2 \, dr - m^{B1Tip} \, g \, {}^E\boldsymbol{v}_{Row}^{S1}(BldFlexL) \cdot \boldsymbol{z}_2 \quad (\text{Row} = 3, 4, \ldots, 14; 16, 17, 18; 22)\tag{23} +$$ + + +$$ +\boldsymbol{F}_r \big|_{AeroB1} = \int_{0}^{BldFlexL} \left[ {}^E\boldsymbol{v}_r^{S1}(r) \cdot \boldsymbol{F}_{AeroB1}^{S1}(r) + {}^E\boldsymbol{\omega}_r^{M1}(r) \cdot \boldsymbol{M}_{AeroB1}^{M1}(r) \right] dr + {}^E\boldsymbol{v}_r^{S1}(BldFlexL) \cdot \boldsymbol{F}_{TipDragB1}^{S1}(BldFlexL)(r = 1, 2, \ldots, 14; 16, 17, 18; Teet)\tag{24} $$ where $F_{A e r o B 1}^{S 1}(r)$ and $M_{A e r o B 1}^{M 1}\left(r\right)$ are aerodynamic forces and pitching moments applied to point S1 on blade 1 respectively expressed per unit span. Note that $M_{A e r o B 1}^{M 1}\left(r\right)$ can include effects of both direct aerodynamic pitching moments (i.e., $\mathrm{Cm})$ ) and aerodynamic pitching moments caused by an aerodynamic offset (i.e., moments due to aerodynamic lift and drag forces acting at a distance away from the center of mass of the blade element along the aerodynamic chord). +其中，$F_{\text{AeroB1}}^{S1}(r)$ 和 $M_{\text{AeroB1}}^{M1}(r)$ 分别表示作用于叶片1上点S1处的单位展向长度上的气动力和气动俯仰力矩。需要注意的是，$M_{\text{AeroB1}}^{M1}(r)$ 可包含两部分气动俯仰力矩的影响：一是直接气动俯仰力矩（即 $\mathrm{Cm}$），二是由气动力（升力和阻力）在叶片元素气动弦线上相对于质心产生的偏移距离所引起的气动俯仰力矩。 +$$ +[C(q,t)]|_{AeroB1} = 0\tag{25} +$$ +$$ +\left\{ -f(\dot{q}, q, t) \right\}|_{AeroB1} (\text{Row}) = \int_{0}^{BldFlexL} \left[ {}^E\boldsymbol{v}_{Row}^{S1}(r) \cdot \boldsymbol{F}_{AeroB1}^{S1}(r) + {}^E\boldsymbol{\omega}_{Row}^{M1}(r) \cdot \boldsymbol{M}_{AeroB1}^{M1}(r) \right] dr + {}^E\boldsymbol{v}_{Row}^{S1}(BldFlexL) \cdot \boldsymbol{F}_{TipDragB1}^{S1}(BldFlexL) \quad (\text{Row} = 1, 2, \ldots, 14; 16, 17, 18; 22)\tag{26} +$$ \ No newline at end of file diff --git a/书籍/力学书籍/OpenFast/FASTKinetics/auto/CalculateAngularPosVelPAcc.md b/书籍/力学书籍/OpenFast/FASTKinetics/auto/CalculateAngularPosVelPAcc.md new file mode 100644 index 0000000..66797e5 --- /dev/null +++ b/书籍/力学书籍/OpenFast/FASTKinetics/auto/CalculateAngularPosVelPAcc.md @@ -0,0 +1,222 @@ + + +```fortran +SUBROUTINE CalculateAngularPosVelPAcc( p, x, CoordSys, RtHSdat ) +!.................................................................................................................................. + + ! 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_RtHndSide), INTENT(INOUT) :: RtHSdat !< data from the RtHndSid module (contains positions to be set) + + !Local variables + + REAL(ReKi) :: AngVelHM (3) ! Angular velocity of eleMent J of blade K (body M) in the hub (body H). +! REAL(ReKi) :: AngVelEN (3) ! Angular velocity of the nacelle (body N) in the inertia frame (body E for earth). + REAL(ReKi) :: AngAccELt (3) ! Portion of the angular acceleration of the low-speed shaft (body L) in the inertia frame (body E for earth) associated with everything but the QD2T()'s. + INTEGER(IntKi) :: J ! Counter for elements + INTEGER(IntKi) :: K ! Counter for blades + + !------------------------------------------------------------------------------------------------- + ! Angular and partial angular velocities + !------------------------------------------------------------------------------------------------- + + ! Define the angular and partial angular velocities of all of the rigid bodies in the inertia frame: + ! NOTE: PAngVelEN(I,D,:) = the Dth-derivative of the partial angular velocity of DOF I for body N in body E. + + RtHSdat%PAngVelEX( :,0,:) = 0.0 + RtHSdat%PAngVelEX(DOF_R ,0,:) = CoordSys%z1 + RtHSdat%PAngVelEX(DOF_P ,0,:) = -CoordSys%z3 + RtHSdat%PAngVelEX(DOF_Y ,0,:) = CoordSys%z2 + RtHSdat%AngVelEX = x%QDT(DOF_R )*RtHSdat%PAngVelEX(DOF_R ,0,:) & + + x%QDT(DOF_P )*RtHSdat%PAngVelEX(DOF_P ,0,:) & + + x%QDT(DOF_Y )*RtHSdat%PAngVelEX(DOF_Y ,0,:) + RtHSdat%AngPosEX = x%QT (DOF_R )*RtHSdat%PAngVelEX(DOF_R ,0,:) & + + x%QT (DOF_P )*RtHSdat%PAngVelEX(DOF_P ,0,:) & + + x%QT (DOF_Y )*RtHSdat%PAngVelEX(DOF_Y ,0,:) + + RtHSdat%PAngVelEB( :,0,:) = RtHSdat%PAngVelEX(:,0,:) + RtHSdat%PAngVelEB(DOF_TFA1,0,:) = -p%TwrFASF(1,p%TTopNode,1)*CoordSys%a3 + RtHSdat%PAngVelEB(DOF_TSS1,0,:) = p%TwrSSSF(1,p%TTopNode,1)*CoordSys%a1 + RtHSdat%PAngVelEB(DOF_TFA2,0,:) = -p%TwrFASF(2,p%TTopNode,1)*CoordSys%a3 + RtHSdat%PAngVelEB(DOF_TSS2,0,:) = p%TwrSSSF(2,p%TTopNode,1)*CoordSys%a1 + RtHSdat%AngVelEB = RtHSdat%AngVelEX + x%QDT(DOF_TFA1)*RtHSdat%PAngVelEB(DOF_TFA1,0,:) & + + x%QDT(DOF_TSS1)*RtHSdat%PAngVelEB(DOF_TSS1,0,:) & + + x%QDT(DOF_TFA2)*RtHSdat%PAngVelEB(DOF_TFA2,0,:) & + + x%QDT(DOF_TSS2)*RtHSdat%PAngVelEB(DOF_TSS2,0,:) + RtHSdat%AngPosXB = x%QT (DOF_TFA1)*RtHSdat%PAngVelEB(DOF_TFA1,0,:) & + + x%QT (DOF_TSS1)*RtHSdat%PAngVelEB(DOF_TSS1,0,:) & + + x%QT (DOF_TFA2)*RtHSdat%PAngVelEB(DOF_TFA2,0,:) & + + x%QT (DOF_TSS2)*RtHSdat%PAngVelEB(DOF_TSS2,0,:) + + RtHSdat%PAngVelEN( :,0,:)= RtHSdat%PAngVelEB(:,0,:) + RtHSdat%PAngVelEN(DOF_Yaw ,0,:)= CoordSys%d2 + RtHSdat%AngVelEN = RtHSdat%AngVelEB + x%QDT(DOF_Yaw )*RtHSdat%PAngVelEN(DOF_Yaw ,0,:) + + RtHSdat%PAngVelER( :,0,:)= RtHSdat%PAngVelEN(:,0,:) + RtHSdat%PAngVelER(DOF_RFrl,0,:)= CoordSys%rfa + RtHSdat%AngVelER = RtHSdat%AngVelEN + x%QDT(DOF_RFrl)*RtHSdat%PAngVelER(DOF_RFrl,0,:) + + RtHSdat%PAngVelEL( :,0,:)= RtHSdat%PAngVelER(:,0,:) + RtHSdat%PAngVelEL(DOF_GeAz,0,:)= CoordSys%c1 + RtHSdat%PAngVelEL(DOF_DrTr,0,:)= CoordSys%c1 + RtHSdat%AngVelEL = RtHSdat%AngVelER + x%QDT(DOF_GeAz)*RtHSdat%PAngVelEL(DOF_GeAz,0,:) & + + x%QDT(DOF_DrTr)*RtHSdat%PAngVelEL(DOF_DrTr,0,:) + + RtHSdat%PAngVelEH( :,0,:)= RtHSdat%PAngVelEL(:,0,:) + RtHSdat%AngVelEH = RtHSdat%AngVelEL +IF ( p%NumBl == 2 ) THEN ! 2-blader + RtHSdat%PAngVelEH(DOF_Teet,0,:)= CoordSys%f2 + RtHSdat%AngVelEH = RtHSdat%AngVelEH + x%QDT(DOF_Teet)*RtHSdat%PAngVelEH(DOF_Teet,0,:) +ENDIF + + RtHSdat%PAngVelEG( :,0,:) = RtHSdat%PAngVelER(:,0,:) + RtHSdat%PAngVelEG(DOF_GeAz,0,:) = p%GBRatio*CoordSys%c1 + RtHSdat%AngVelEG = RtHSdat%AngVelER + x%QDT(DOF_GeAz)*RtHSdat%PAngVelEG(DOF_GeAz,0,:) + + RtHSdat%PAngVelEA( :,0,:) = RtHSdat%PAngVelEN(:,0,:) + RtHSdat%PAngVelEA(DOF_TFrl,0,:) = CoordSys%tfa + RtHSdat%AngVelEA = RtHSdat%AngVelEN + x%QDT(DOF_TFrl)*RtHSdat%PAngVelEA(DOF_TFrl,0,:) + + + + ! Define the 1st derivatives of the partial angular velocities of all + ! of the rigid bodies in the inertia frame and the portion of the angular + ! acceleration of the rigid bodies in the inertia frame associated with + ! everything but the QD2T()'s: + + RtHSdat%PAngVelEX( :,1,:) = 0.0 + RtHSdat%AngAccEXt = 0.0 + + RtHSdat%PAngVelEB( :,1,:) = RtHSdat%PAngVelEX(:,1,:) + RtHSdat%PAngVelEB(DOF_TFA1,1,:) = CROSS_PRODUCT( RtHSdat%AngVelEX, RtHSdat%PAngVelEB(DOF_TFA1,0,:) ) + RtHSdat%PAngVelEB(DOF_TSS1,1,:) = CROSS_PRODUCT( RtHSdat%AngVelEX, RtHSdat%PAngVelEB(DOF_TSS1,0,:) ) + RtHSdat%PAngVelEB(DOF_TFA2,1,:) = CROSS_PRODUCT( RtHSdat%AngVelEX, RtHSdat%PAngVelEB(DOF_TFA2,0,:) ) + RtHSdat%PAngVelEB(DOF_TSS2,1,:) = CROSS_PRODUCT( RtHSdat%AngVelEX, RtHSdat%PAngVelEB(DOF_TSS2,0,:) ) + RtHSdat%AngAccEBt = RtHSdat%AngAccEXt + x%QDT(DOF_TFA1)*RtHSdat%PAngVelEB(DOF_TFA1,1,:) & + + x%QDT(DOF_TSS1)*RtHSdat%PAngVelEB(DOF_TSS1,1,:) & + + x%QDT(DOF_TFA2)*RtHSdat%PAngVelEB(DOF_TFA2,1,:) & + + x%QDT(DOF_TSS2)*RtHSdat%PAngVelEB(DOF_TSS2,1,:) + + RtHSdat%PAngVelEN( :,1,:) = RtHSdat%PAngVelEB(:,1,:) + RtHSdat%PAngVelEN(DOF_Yaw ,1,:) = CROSS_PRODUCT( RtHSdat%AngVelEB, RtHSdat%PAngVelEN(DOF_Yaw ,0,:) ) + RtHSdat%AngAccENt = RtHSdat%AngAccEBt + x%QDT(DOF_Yaw )*RtHSdat%PAngVelEN(DOF_Yaw ,1,:) + + RtHSdat%PAngVelER( :,1,:) = RtHSdat%PAngVelEN(:,1,:) + RtHSdat%PAngVelER(DOF_RFrl,1,:) = CROSS_PRODUCT( RtHSdat%AngVelEN, RtHSdat%PAngVelER(DOF_RFrl,0,:) ) + RtHSdat%AngAccERt = RtHSdat%AngAccENt + x%QDT(DOF_RFrl)*RtHSdat%PAngVelER(DOF_RFrl,1,:) + + RtHSdat%PAngVelEL( :,1,:) = RtHSdat%PAngVelER(:,1,:) + RtHSdat%PAngVelEL(DOF_GeAz,1,:) = CROSS_PRODUCT( RtHSdat%AngVelER, RtHSdat%PAngVelEL(DOF_GeAz,0,:) ) + RtHSdat%PAngVelEL(DOF_DrTr,1,:) = CROSS_PRODUCT( RtHSdat%AngVelER, RtHSdat%PAngVelEL(DOF_DrTr,0,:) ) + AngAccELt = RtHSdat%AngAccERt + x%QDT(DOF_GeAz)*RtHSdat%PAngVelEL(DOF_GeAz,1,:) & + + x%QDT(DOF_DrTr)*RtHSdat%PAngVelEL(DOF_DrTr,1,:) + + RtHSdat%PAngVelEH( :,1,:) = RtHSdat%PAngVelEL(:,1,:) + RtHSdat%AngAccEHt = AngAccELt +IF ( p%NumBl == 2 ) THEN ! 2-blader + RtHSdat%PAngVelEH(DOF_Teet,1,:) = CROSS_PRODUCT( RtHSdat%AngVelEH, RtHSdat%PAngVelEH(DOF_Teet,0,:) ) + RtHSdat%AngAccEHt = RtHSdat%AngAccEHt + x%QDT(DOF_Teet)*RtHSdat%PAngVelEH(DOF_Teet,1,:) +ENDIF + + RtHSdat%PAngVelEG( :,1,:) = RtHSdat%PAngVelER(:,1,:) + RtHSdat%PAngVelEG(DOF_GeAz,1,:) = CROSS_PRODUCT( RtHSdat%AngVelER, RtHSdat%PAngVelEG(DOF_GeAz,0,:) ) + RtHSdat%AngAccEGt = RtHSdat%AngAccERt + x%QDT(DOF_GeAz)*RtHSdat%PAngVelEG(DOF_GeAz,1,:) + + RtHSdat%PAngVelEA( :,1,:) = RtHSdat%PAngVelEN(:,1,:) + RtHSdat%PAngVelEA(DOF_TFrl,1,:) = CROSS_PRODUCT( RtHSdat%AngVelEN, RtHSdat%PAngVelEA(DOF_TFrl,0,:) ) + RtHSdat%AngAccEAt = RtHSdat%AngAccENt + x%QDT(DOF_TFrl)*RtHSdat%PAngVelEA(DOF_TFrl,1,:) + + + + DO K = 1,p%NumBl ! Loop through all blades + + DO J = 0,p%TipNode ! Loop through the blade nodes / elements + ! Define the partial angular velocities of the current node (body M(RNodes(J))) in the inertia frame: + ! NOTE: PAngVelEM(K,J,I,D,:) = the Dth-derivative of the partial angular velocity + ! of DOF I for body M of blade K, element J in body E. + + RtHSdat%PAngVelEM(K,J, :,0,:) = RtHSdat%PAngVelEH(:,0,:) + RtHSdat%PAngVelEM(K,J,DOF_BF(K,1),0,:) = - p%TwistedSF(K,2,1,J,1)*CoordSys%j1(K,:) & + + p%TwistedSF(K,1,1,J,1)*CoordSys%j2(K,:) + RtHSdat%PAngVelEM(K,J,DOF_BF(K,2),0,:) = - p%TwistedSF(K,2,2,J,1)*CoordSys%j1(K,:) & + + p%TwistedSF(K,1,2,J,1)*CoordSys%j2(K,:) + RtHSdat%PAngVelEM(K,J,DOF_BE(K,1),0,:) = - p%TwistedSF(K,2,3,J,1)*CoordSys%j1(K,:) & + + p%TwistedSF(K,1,3,J,1)*CoordSys%j2(K,:) + AngVelHM = x%QDT(DOF_BF(K,1))*RtHSdat%PAngVelEM(K,J,DOF_BF(K,1),0,:) & + + x%QDT(DOF_BF(K,2))*RtHSdat%PAngVelEM(K,J,DOF_BF(K,2),0,:) & + + x%QDT(DOF_BE(K,1))*RtHSdat%PAngVelEM(K,J,DOF_BE(K,1),0,:) + RtHSdat%AngVelEM(:,J,K ) = RtHSdat%AngVelEH + AngVelHM + RtHSdat%AngPosHM(:,K,J ) = x%QT (DOF_BF(K,1))*RtHSdat%PAngVelEM(K,J,DOF_BF(K,1),0,:) & + + x%QT (DOF_BF(K,2))*RtHSdat%PAngVelEM(K,J,DOF_BF(K,2),0,:) & + + x%QT (DOF_BE(K,1))*RtHSdat%PAngVelEM(K,J,DOF_BE(K,1),0,:) + RtHSdat%AngAccEKt(:,J ,K) = RtHSdat%AngAccEHt + x%QDT(DOF_BF(K,1))*RtHSdat%PAngVelEM(K,J,DOF_BF(K,1),1,:) & + + x%QDT(DOF_BF(K,2))*RtHSdat%PAngVelEM(K,J,DOF_BF(K,2),1,:) & + + x%QDT(DOF_BE(K,1))*RtHSdat%PAngVelEM(K,J,DOF_BE(K,1),1,:) + + ! Define the 1st derivatives of the partial angular velocities of the current node (body M(RNodes(J))) in the inertia frame: + + ! NOTE: These are currently unused by the code, therefore, they need not + ! be calculated. Thus, they are currently commented out. If it + ! turns out that they are ever needed (i.e., if inertias of the + ! blade elements are ever added, etc...) simply uncomment out these computations: + ! RtHSdat%PAngVelEM(K,J, :,1,:) = RtHSdat%PAngVelEH(:,1,:) + ! RtHSdat%PAngVelEM(K,J,DOF_BF(K,1),1,:) = CROSS_PRODUCT( RtHSdat%AngVelEH, PAngVelEM(K,J,DOF_BF(K,1),0,:) ) + ! RtHSdat%PAngVelEM(K,J,DOF_BF(K,2),1,:) = CROSS_PRODUCT( RtHSdat%AngVelEH, PAngVelEM(K,J,DOF_BF(K,2),0,:) ) + ! RtHSdat%PAngVelEM(K,J,DOF_BE(K,1),1,:) = CROSS_PRODUCT( RtHSdat%AngVelEH, PAngVelEM(K,J,DOF_BE(K,1),0,:) ) + + + END DO !J = 1,p%BldNodes ! Loop through the blade nodes / elements + + END DO !K = 1,p%NumBl + + + !............... + ! tower values: + !............... + + DO J = 0,p%TwrNodes ! Loop through the tower nodes / elements + + ! Define the partial angular velocities (and their 1st derivatives) of the + ! current node (body F(HNodes(J)) in the inertia frame. + ! Also define the overall angular velocity of the current node in the inertia frame. + ! Also, define the portion of the angular acceleration of the current node + ! in the inertia frame associated with everything but the QD2T()'s: + + ! NOTE: PAngVelEF(J,I,D,:) = the Dth-derivative of the partial angular velocity + ! of DOF I for body F of element J in body E. + + RtHSdat%PAngVelEF (J, :,0,:) = RtHSdat%PAngVelEX(:,0,:) + RtHSdat%PAngVelEF (J,DOF_TFA1,0,:) = -p%TwrFASF(1,J,1)*CoordSys%a3 + RtHSdat%PAngVelEF (J,DOF_TSS1,0,:) = p%TwrSSSF(1,J,1)*CoordSys%a1 + RtHSdat%PAngVelEF (J,DOF_TFA2,0,:) = -p%TwrFASF(2,J,1)*CoordSys%a3 + RtHSdat%PAngVelEF (J,DOF_TSS2,0,:) = p%TwrSSSF(2,J,1)*CoordSys%a1 + + RtHSdat%PAngVelEF (J, :,1,:) = RtHSdat%PAngVelEX(:,1,:) + RtHSdat%PAngVelEF (J,DOF_TFA1,1,:) = CROSS_PRODUCT( RtHSdat%AngVelEX , RtHSdat%PAngVelEF(J,DOF_TFA1,0,:) ) + RtHSdat%PAngVelEF (J,DOF_TSS1,1,:) = CROSS_PRODUCT( RtHSdat%AngVelEX , RtHSdat%PAngVelEF(J,DOF_TSS1,0,:) ) + RtHSdat%PAngVelEF (J,DOF_TFA2,1,:) = CROSS_PRODUCT( RtHSdat%AngVelEX , RtHSdat%PAngVelEF(J,DOF_TFA2,0,:) ) + RtHSdat%PAngVelEF (J,DOF_TSS2,1,:) = CROSS_PRODUCT( RtHSdat%AngVelEX , RtHSdat%PAngVelEF(J,DOF_TSS2,0,:) ) + + + RtHSdat%AngVelEF (:,J) = RtHSdat%AngVelEX + x%QDT(DOF_TFA1)*RtHSdat%PAngVelEF(J,DOF_TFA1,0,:) & + + x%QDT(DOF_TSS1)*RtHSdat%PAngVelEF(J,DOF_TSS1,0,:) & + + x%QDT(DOF_TFA2)*RtHSdat%PAngVelEF(J,DOF_TFA2,0,:) & + + x%QDT(DOF_TSS2)*RtHSdat%PAngVelEF(J,DOF_TSS2,0,:) + + RtHSdat%AngPosXF (:,J) = x%QT (DOF_TFA1)*RtHSdat%PAngVelEF(J,DOF_TFA1,0,:) & + + x%QT (DOF_TSS1)*RtHSdat%PAngVelEF(J,DOF_TSS1,0,:) & + + x%QT (DOF_TFA2)*RtHSdat%PAngVelEF(J,DOF_TFA2,0,:) & + + x%QT (DOF_TSS2)*RtHSdat%PAngVelEF(J,DOF_TSS2,0,:) + RtHSdat%AngPosEF (:,J) = RtHSdat%AngPosEX + RtHSdat%AngPosXF(:,J) + RtHSdat%AngAccEFt(:,J) = RtHSdat%AngAccEXt + x%QDT(DOF_TFA1)*RtHSdat%PAngVelEF(J,DOF_TFA1,1,:) & + + x%QDT(DOF_TSS1)*RtHSdat%PAngVelEF(J,DOF_TSS1,1,:) & + + x%QDT(DOF_TFA2)*RtHSdat%PAngVelEF(J,DOF_TFA2,1,:) & + + x%QDT(DOF_TSS2)*RtHSdat%PAngVelEF(J,DOF_TSS2,1,:) + + END DO ! J + + +END SUBROUTINE CalculateAngularPosVelPAcc +``` \ No newline at end of file diff --git a/书籍/力学书籍/OpenFast/FASTKinetics/auto/CalculateForcesMoments.md b/书籍/力学书籍/OpenFast/FASTKinetics/auto/CalculateForcesMoments.md new file mode 100644 index 0000000..4aa1426 --- /dev/null +++ b/书籍/力学书籍/OpenFast/FASTKinetics/auto/CalculateForcesMoments.md @@ -0,0 +1,634 @@ + +```fortran +SUBROUTINE CalculateForcesMoments( p, x, CoordSys, u, RtHSdat ) +!.................................................................................................................................. + + ! 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) & nacelle forces/moments + TYPE(ED_RtHndSide), INTENT(INOUT) :: RtHSdat !< data from the RtHndSid module (contains positions to be set) + + ! 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) :: TmpVec2 (3) ! A temporary vector used in various computations. + REAL(ReKi) :: TmpVec3 (3) ! A temporary vector used in various computations. + REAL(ReKi) :: TmpVec4 (3) ! A temporary vector used in various computations. + REAL(ReKi) :: TmpVec5 (3) ! A temporary vector used in various computations. + REAL(ReKi) :: Force(3) ! External force (e.g. from AeroDyn) + REAL(ReKi) :: Moment(3) ! External moment (e.g. from AeroDyn) + 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) :: NodeNum ! Node number for blade element (on a single mesh) + +!..................................... +! Compute forces and moments from properties related to Aero inputs +!FSTipDrag +!FSAero and MMAero +!..................................... + DO K = 1,p%NumBl ! Loop through all blades + + ! Calculate the tip drag forces if necessary: + !bjj: add this back when we've figured out how to handle the tip brakes: + !RtHSdat%FSTipDrag(:,K) = m%CoordSys%m2(K,p%BldNodes,:)*SIGN( 0.5*p%AirDens*(RtHSdat%LinVelESm2(K)**2)*u%TBDrCon(K), -1.*RtHSdat%LinVelESm2(K) ) + RtHSdat%FSTipDrag = 0.0_ReKi ! Calculate the tip drag forces if necessary + + + ! Calculate the normal and tangential aerodynamic forces and the aerodynamic + ! pitching moment at the current element per unit span by calling AeroDyn, + ! if necessary: + + DO J = 1,p%BldNodes ! Loop through the blade nodes / elements + + + ! Calculate the aerodynamic pitching moment arm (i.e., the position vector + ! from point S on the blade to the aerodynamic center of the element): + + RtHSdat%rSAerCen(:,J,K) = p%rSAerCenn1(K,J)*CoordSys%n1(K,J,:) + p%rSAerCenn2(K,J)*CoordSys%n2(K,J,:) + +! rPAerCen = m%RtHS%rPQ + m%RtHS%rQS(:,K,J) + m%RtHS%rSAerCen(:,J,K) ! Position vector from teeter pin (point P) to blade analysis node aerodynamic center. +! rAerCen = m%RtHS%rS (:,K,J) + m%RtHS%rSAerCen(:,J,K) ! Position vector from inertial frame origin to blade analysis node aerodynamic center. + + + ! fill FSAero() and MMAero() with the forces resulting from inputs u%BladeLn2Mesh(K)%Force(1:2,:) and u%BladeLn2Mesh(K)%Moment(3,:): + ! [except, we're ignoring the additional nodes we added on the mesh end points] + + NodeNum = J ! we're ignoring the root and tip + + if (p%UseAD14) then + RtHSdat%FSAero(:,K,J) = ( u%BladePtLoads(K)%Force(1,NodeNum) * CoordSys%te1(K,J,:) & + + u%BladePtLoads(K)%Force(2,NodeNum) * CoordSys%te2(K,J,:) ) / p%DRNodes(J) + + RtHSdat%MMAero(:,K,J) = CROSS_PRODUCT( RtHSdat%rSAerCen(:,J,K), RtHSdat%FSAero(:,K,J) )& + + u%BladePtLoads(K)%Moment(3,NodeNum)/p%DRNodes(J) * CoordSys%te3(K,J,:) + else + RtHSdat%FSAero(1,K,J) = u%BladePtLoads(K)%Force(1,NodeNum) / p%DRNodes(J) + RtHSdat%FSAero(2,K,J) = u%BladePtLoads(K)%Force(3,NodeNum) / p%DRNodes(J) + RtHSdat%FSAero(3,K,J) = -u%BladePtLoads(K)%Force(2,NodeNum) / p%DRNodes(J) + + RtHSdat%MMAero(1,K,J) = u%BladePtLoads(K)%Moment(1,NodeNum) / p%DRNodes(J) + RtHSdat%MMAero(2,K,J) = u%BladePtLoads(K)%Moment(3,NodeNum) / p%DRNodes(J) + RtHSdat%MMAero(3,K,J) = -u%BladePtLoads(K)%Moment(2,NodeNum) / p%DRNodes(J) + end if + + + END DO !J + END DO ! K + + +!..................................... +! PFrcS0B and PMomH0B +!..................................... +! RtHSdat%PFrcS0B = 0.0 +! RtHSdat%PMomH0B = 0.0 +DO K = 1,p%NumBl ! Loop through all blades + + ! Initialize the partial forces and moments (including those associated + ! with the QD2T()'s and those that are not) at the blade root (point S(0)) + ! using the tip brake effects: + + RtHSdat%PFrcS0B(:,K,:) = 0.0 ! Initialize these partial + RtHSdat%PMomH0B(:,K,:) = 0.0 ! forces and moments to zero + DO I = 1,p%DOFs%NPSE(K) ! Loop through all active (enabled) DOFs that contribute to the QD2T-related linear accelerations of blade K + + TmpVec1 = -p%TipMass(K)*RtHSdat%PLinVelES(K,p%TipNode,p%DOFs%PSE(K,I),0,:) ! The portion of PFrcS0B associated with the tip brake + + RtHSdat%PFrcS0B(:,K,p%DOFs%PSE(K,I)) = TmpVec1 + RtHSdat%PMomH0B(:,K,p%DOFs%PSE(K,I)) = CROSS_PRODUCT( RtHSdat%rS0S(:,K,p%TipNode), TmpVec1 ) ! The portion of PMomH0B associated with the tip brake + + ENDDO ! I - All active (enabled) DOFs that contribute to the QD2T-related linear accelerations of blade K + + + DO J = 1,p%BldNodes ! Loop through the blade nodes / elements + + ! Integrate to find the partial forces and moments (including those associated + ! with the QD2T()'s and those that are not) at the blade root (point S(0)): + + DO I = 1,p%DOFs%NPSE(K) ! Loop through all active (enabled) DOFs that contribute to the QD2T-related linear accelerations of blade K + + TmpVec1 = -p%BElmntMass(J,K)*RtHSdat%PLinVelES(K,J,p%DOFs%PSE(K,I),0,:) ! The portion of PFrcS0B associated with blade element J + + RtHSdat%PFrcS0B(:,K,p%DOFs%PSE(K,I)) = RtHSdat%PFrcS0B(:,K,p%DOFs%PSE(K,I)) + TmpVec1 + RtHSdat%PMomH0B(:,K,p%DOFs%PSE(K,I)) = RtHSdat%PMomH0B(:,K,p%DOFs%PSE(K,I)) + & + CROSS_PRODUCT( RtHSdat%rS0S(:,K,J), TmpVec1 ) ! The portion of PMomH0B associated with blade element J + + ENDDO ! I - All active (enabled) DOFs that contribute to the QD2T-related linear accelerations of blade K + END DO + + + END DO + + +!..................................... +! FrcS0Bt and MomH0Bt +!..................................... + DO K = 1,p%NumBl ! Loop through all blades + + 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 + RtHSdat%FrcS0Bt(:,K) = TmpVec1 + RtHSdat%MomH0Bt(:,K) = CROSS_PRODUCT( RtHSdat%rS0S(:,K,p%TipNode), TmpVec1 ) ! The portion of MomH0Bt associated with the tip brake + + DO J = 1,p%BldNodes ! Loop through the blade nodes / elements + + 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 + TmpVec2 = CROSS_PRODUCT( RtHSdat%rS0S(:,K,J), TmpVec1 ) ! The portion of MomH0Bt associated with blade element J + TmpVec3 = RtHSdat%MMAero(:,K,J)*p%DRNodes(J) ! The total external moment applied to blade element J + + RtHSdat%FrcS0Bt(:,K) = RtHSdat%FrcS0Bt(:,K) + TmpVec1 + RtHSdat%MomH0Bt(:,K) = RtHSdat%MomH0Bt(:,K) + TmpVec2 + TmpVec3 + + END DO !J + + END DO !K + + +!..................................... +! PFrcPRot AND PMomLPRot: +! ( requires PFrcS0B and PMomH0B) +!..................................... + ! Initialize the partial forces and moments (including those associated + ! with the QD2T()'s and those that are not) at the teeter pin (point P) using the hub mass effects: + + RtHSdat%PFrcPRot = 0.0 ! Initialize these partial + RtHSdat%PMomLPRot = 0.0 ! forces and moments to zero + DO I = 1,p%DOFs%NPCE ! Loop through all active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the hub center of mass (point C) + + TmpVec1 = -p%HubMass*RtHSdat%PLinVelEC(p%DOFs%PCE(I),0,:) ! The portion of PFrcPRot associated with the HubMass + TmpVec2 = CROSS_PRODUCT( RtHSdat%rPC, TmpVec1 ) ! The portion of PMomLPRot associated with the HubMass + + RtHSdat%PFrcPRot (:,p%DOFs%PCE(I)) = TmpVec1 + RtHSdat%PMomLPRot(:,p%DOFs%PCE(I)) = TmpVec2 - p%Hubg1Iner*CoordSys%g1*DOT_PRODUCT( CoordSys%g1, RtHSdat%PAngVelEH(p%DOFs%PCE(I),0,:) ) & + - p%Hubg2Iner*CoordSys%g2*DOT_PRODUCT( CoordSys%g2, RtHSdat%PAngVelEH(p%DOFs%PCE(I),0,:) ) + + ENDDO ! I - All active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the hub center of mass (point C) + + + DO K = 1,p%NumBl ! Loop through all blades + + ! Calculate the position vector from the teeter pin to the blade root: + + !rPS0(:,K) = RtHSdat%rPQ + p%HubRad*CoordSys%j3(K,:) ! Position vector from teeter pin (point P) to blade root (point S(0)). + + ! Add the blade effects to the partial forces and moments (including those associated with the QD2T()'s and those that are + ! not) at the teeter pin (point P): + + DO I = 1,p%DOFs%NPSE(K) ! Loop through all active (enabled) DOFs that contribute to the QD2T-related linear accelerations of blade K + + TmpVec = CROSS_PRODUCT( RtHSdat%rPS0(:,K), RtHSdat%PFrcS0B(:,K,p%DOFs%PSE(K,I)) ) ! The portion of PMomLPRot associated with PFrcS0B. + + RtHSdat%PFrcPRot (:,p%DOFs%PSE(K,I)) = RtHSdat%PFrcPRot (:,p%DOFs%PSE(K,I)) + RtHSdat%PFrcS0B(:,K,p%DOFs%PSE(K,I)) + RtHSdat%PMomLPRot(:,p%DOFs%PSE(K,I)) = RtHSdat%PMomLPRot(:,p%DOFs%PSE(K,I)) + RtHSdat%PMomH0B(:,K,p%DOFs%PSE(K,I))+TmpVec + + ENDDO ! I - All active (enabled) DOFs that contribute to the QD2T-related linear accelerations of blade K + + + END DO ! K + !PRINT '(F30.20)', RtHSdat%PMomLPRot(:,11) +!..................................... +! FrcPRott and MomLPRott: +! (requires FrcS0Bt and MomH0Bt) +!..................................... + + TmpVec1 = -p%HubMass*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccECt ) ! The portion of FrcPRott associated with the HubMass + TmpVec2 = CROSS_PRODUCT( RtHSdat%rPC, TmpVec1 ) ! The portion of MomLPRott associated with the HubMass + TmpVec = p%Hubg1Iner*CoordSys%g1*DOT_PRODUCT( CoordSys%g1, RtHSdat%AngVelEH ) & ! = ( Hub inertia dyadic ) dot ( angular velocity of hub in the inertia frame ) + + p%Hubg2Iner*CoordSys%g2*DOT_PRODUCT( CoordSys%g2, RtHSdat%AngVelEH ) + TmpVec3 = CROSS_PRODUCT( -RtHSdat%AngVelEH, TmpVec ) ! = ( -angular velocity of hub in the inertia frame ) cross ( TmpVec ) + + RtHSdat%FrcPRott(1) = TmpVec1(1) + u%HubPtLoad%Force(1,1) + RtHSdat%FrcPRott(2) = TmpVec1(2) + u%HubPtLoad%Force(3,1) + RtHSdat%FrcPRott(3) = TmpVec1(3) - u%HubPtLoad%Force(2,1) + + RtHSdat%MomLPRott = TmpVec2 + TmpVec3 - p%Hubg1Iner*CoordSys%g1*DOT_PRODUCT( CoordSys%g1, RtHSdat%AngAccEHt ) & + - p%Hubg2Iner*CoordSys%g2*DOT_PRODUCT( CoordSys%g2, RtHSdat%AngAccEHt ) + + RtHSdat%MomLPRott(1) = RtHSdat%MomLPRott(1) + u%HubPtLoad%Moment(1,1) + RtHSdat%MomLPRott(2) = RtHSdat%MomLPRott(2) + u%HubPtLoad%Moment(3,1) + RtHSdat%MomLPRott(3) = RtHSdat%MomLPRott(3) - u%HubPtLoad%Moment(2,1) + + DO K = 1,p%NumBl ! Loop through all blades + + ! Calculate the position vector from the teeter pin to the blade root: + + !rPS0 = RtHSdat%rPQ + p%HubRad*m%CoordSys%j3(K,:) ! Position vector from teeter pin (point P) to blade root (point S(0)). + + TmpVec = CROSS_PRODUCT( RtHSdat%rPS0(:,K), RtHSdat%FrcS0Bt(:,K) ) ! The portion of MomLPRott associated with FrcS0Bt. + + RtHSdat%FrcPRott = RtHSdat%FrcPRott + RtHSdat%FrcS0Bt(:,K) + RtHSdat%MomLPRott = RtHSdat%MomLPRott + RtHSdat%MomH0Bt(:,K) + TmpVec + + END DO ! K + + + ! Define the partial forces and moments (including those associated with + ! the QD2T()'s and those that are not) at the specified point on the + ! rotor-furl axis (point V) / nacelle (body N) using the structure that + ! furls with the rotor, generator, and rotor effects. + +!..................................... +! PMomNGnRt and PFrcVGnRt +! (requires PMomLPRot and PFrcPRot) +!..................................... + RtHSdat%PFrcVGnRt = RtHSdat%PFrcPRot ! Initialize these partial forces and + RtHSdat%PMomNGnRt = RtHSdat%PMomLPRot ! moments using the rotor effects + DO I = 1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs + + TmpVec = CROSS_PRODUCT( RtHSdat%rVP, RtHSdat%PFrcPRot(:,p%DOFs%SrtPS(I)) ) ! The portion of PMomNGnRt associated with the PFrcPRot + + RtHSdat%PMomNGnRt(:,p%DOFs%SrtPS(I)) = RtHSdat%PMomNGnRt(:,p%DOFs%SrtPS(I)) + TmpVec + + ENDDO ! I - All active (enabled) DOFs + DO I = 1,p%DOFs%NPDE ! Loop through all active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the center of mass of the structure that furls with the rotor (not including rotor) (point D) + + TmpVec1 = -p%RFrlMass*RtHSdat%PLinVelED(p%DOFs%PDE(I) ,0,:) ! The portion of PFrcVGnRt associated with the RFrlMass + TmpVec2 = CROSS_PRODUCT( RtHSdat%rVD, TmpVec1 ) ! The portion of PMomNGnRt associated with the RFrlMass + + RtHSdat%PFrcVGnRt(:,p%DOFs%PDE(I) ) = RtHSdat%PFrcVGnRt(:,p%DOFs%PDE(I) ) + TmpVec1 + + RtHSdat%PMomNGnRt(:,p%DOFs%PDE(I) ) = RtHSdat%PMomNGnRt(:,p%DOFs%PDE(I) ) + TmpVec2 & + - p%RrfaIner*CoordSys%rfa*DOT_PRODUCT( CoordSys%rfa, RtHSdat%PAngVelER(p%DOFs%PDE(I) ,0,:) ) & + - p%GenIner*CoordSys%c1 *DOT_PRODUCT( CoordSys%c1 , RtHSdat%PAngVelEG(p%DOFs%PDE(I) ,0,:) ) + + ENDDO ! I - All active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the center of mass of the structure that furls with the rotor (not including rotor) (point D) + IF ( p%DOF_Flag(DOF_GeAz) ) THEN + + RtHSdat%PMomNGnRt(:,DOF_GeAz) = RtHSdat%PMomNGnRt(:,DOF_GeAz) & ! The previous loop (DO I = 1,NPDE) misses the DOF_GeAz-contribution to: ( Generator inertia dyadic ) dot ( partial angular velocity of the generator in the inertia frame ) + - p%GenIner*CoordSys%c1 *DOT_PRODUCT( CoordSys%c1, RtHSdat%PAngVelEG(DOF_GeAz,0,:) ) ! Thus, add this contribution if necessary. + + ENDIF + +!..................................... +! FrcVGnRtt and MomNGnRtt +! (requires FrcPRott and MomLPRott) +!..................................... + TmpVec1 = -p%RFrlMass*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccEDt ) ! The portion of FrcVGnRtt associated with the RFrlMass + TmpVec2 = CROSS_PRODUCT( RtHSdat%rVD , TmpVec1 ) ! The portion of MomNGnRtt associated with the RFrlMass + TmpVec3 = CROSS_PRODUCT( RtHSdat%rVP , RtHSdat%FrcPRott ) ! The portion of MomNGnRtt associated with the FrcPRott + TmpVec = p%RrfaIner*CoordSys%rfa*DOT_PRODUCT( CoordSys%rfa, RtHSdat%AngVelER ) ! = ( R inertia dyadic ) dot ( angular velocity of structure that furls with the rotor in the inertia frame ) + TmpVec4 = CROSS_PRODUCT( -RtHSdat%AngVelER, TmpVec ) ! = ( -angular velocity of structure that furls with the rotor in the inertia frame ) cross ( TmpVec ) + TmpVec = p%GenIner*CoordSys%c1* DOT_PRODUCT( CoordSys%c1 , RtHSdat%AngVelEG ) ! = ( Generator inertia dyadic ) dot ( angular velocity of generator in the inertia frame ) + TmpVec5 = CROSS_PRODUCT( -RtHSdat%AngVelEG, TmpVec ) ! = ( -angular velocity of generator in the inertia frame ) cross ( TmpVec ) + + RtHSdat%FrcVGnRtt = RtHSdat%FrcPRott + TmpVec1 + RtHSdat%MomNGnRtt = RtHSdat%MomLPRott + TmpVec2 + TmpVec3 + TmpVec4 + TmpVec5 & + - p%RrfaIner*CoordSys%rfa*DOT_PRODUCT( CoordSys%rfa, RtHSdat%AngAccERt ) & + - p%GenIner*CoordSys%c1 *DOT_PRODUCT( CoordSys%c1 , RtHSdat%AngAccEGt ) + + +!..................................... +! PFrcWTail and PMomNTail +!..................................... + ! Define the partial forces and moments (including those associated with the QD2T()'s and + ! those that are not) at the specified point on the tail-furl axis (point W) / nacelle (body N) using the tail effects. + + RtHSdat%PFrcWTail = 0.0 ! Initialize these partial + RtHSdat%PMomNTail = 0.0 ! forces and moments to zero + DO I = 1,p%DOFs%NPIE ! Loop through all active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the tail boom center of mass (point I) + + TmpVec1 = -p%BoomMass*RtHSdat%PLinVelEI(p%DOFs%PIE(I),0,:) ! The portion of PFrcWTail associated with the BoomMass + TmpVec2 = -p%TFinMass*RtHSdat%PLinVelEJ(p%DOFs%PIE(I),0,:) ! The portion of PFrcWTail associated with the TFinMass + TmpVec3 = CROSS_PRODUCT( RtHSdat%rWI, TmpVec1 ) ! The portion of PMomNTail associated with the BoomMass + TmpVec4 = CROSS_PRODUCT( RtHSdat%rWJ, TmpVec2 ) ! The portion of PMomNTail associated with the TFinMass + + RtHSdat%PFrcWTail(:,p%DOFs%PIE(I)) = TmpVec1 + TmpVec2 + RtHSdat%PMomNTail(:,p%DOFs%PIE(I)) = TmpVec3 + TmpVec4 - p%AtfaIner*CoordSys%tfa* & + DOT_PRODUCT( CoordSys%tfa, RtHSdat%PAngVelEA(p%DOFs%PIE(I),0,:) ) + + ENDDO ! I - All active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the tail boom center of mass (point I) + +!..................................... +! FrcWTailt and MomNTailt - Forces on the tailfin +!..................................... + ! Aerodynamic loads on TailFin CM (point K), with change of coordinate system + Force(1:3) = (/ u%TFinCMLoads%Force (1,1), u%TFinCMLoads%Force (3,1), -u%TFinCMLoads%Force (2,1) /) + Moment(1:3) = (/ u%TFinCMLoads%Moment(1,1), u%TFinCMLoads%Moment(3,1), -u%TFinCMLoads%Moment(2,1) /) + + TmpVec1 = -p%BoomMass*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccEIt ) ! The portion of FrcWTailt associated with the BoomMass + TmpVec2 = -p%TFinMass*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccEJt ) ! The portion of FrcWTailt associated with the TFinMass + TmpVec3 = CROSS_PRODUCT( RtHSdat%rWI , TmpVec1 ) ! The portion of MomNTailt associated with the BoomMass + TmpVec4 = CROSS_PRODUCT( RtHSdat%rWJ , TmpVec2 ) ! The portion of MomNTailt associated with the TFinMass + TmpVec = p%AtfaIner*CoordSys%tfa*DOT_PRODUCT( CoordSys%tfa, RtHSdat%AngVelEA ) ! = ( A inertia dyadic ) dot ( angular velocity of the tail in the inertia frame ) + TmpVec5 = CROSS_PRODUCT( -RtHSdat%AngVelEA, TmpVec ) ! = ( -angular velocity of the tail in the inertia frame ) cross ( TmpVec ) + + RtHSdat%FrcWTailt = Force + TmpVec1 + TmpVec2 + RtHSdat%MomNTailt = Moment + TmpVec3 + TmpVec4 + TmpVec5 & + + CROSS_PRODUCT( RtHSdat%rWJ , Force ) & ! The portion of MomNTailt associated with Force with lever arm WK + - p%AtfaIner*CoordSys%tfa*DOT_PRODUCT( CoordSys%tfa, RtHSdat%AngAccEAt ) + +!..................................... +! PFrcONcRt and PMomBNcRt +! (requires PFrcVGnRt, PMomNGnRt, PFrcWTail, PMomNTail, ) +!..................................... + + ! Define the partial forces and moments (including those associated with + ! the QD2T()'s and those that are not) at the yaw bearing (point O) / + ! base plate (body B) using the nacelle, generator, rotor, and tail effects. + + RtHSdat%PFrcONcRt = RtHSdat%PFrcVGnRt + RtHSdat%PFrcWTail ! Initialize these partial forces and moments using + RtHSdat%PMomBNcRt = RtHSdat%PMomNGnRt + RtHSdat%PMomNTail ! the rotor, rotor-furl, generator, and tail effects + DO I = 1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs + + TmpVec = CROSS_PRODUCT( RtHSdat%rOV, RtHSdat%PFrcVGnRt(:,p%DOFs%SrtPS(I)) ) ! The portion of PMomBNcRt associated with the PFrcVGnRt + + RtHSdat%PMomBNcRt(:,p%DOFs%SrtPS(I)) = RtHSdat%PMomBNcRt(:,p%DOFs%SrtPS(I)) + TmpVec + + ENDDO ! I - All active (enabled) DOFs + DO I = 1,p%DOFs%NPIE ! Loop through all active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the tail boom center of mass (point I) + + TmpVec = CROSS_PRODUCT( RtHSdat%rOW, RtHSdat%PFrcWTail(:,p%DOFs%PIE(I) ) ) ! The portion of PMomBNcRt associated with the PFrcWTail + + RtHSdat%PMomBNcRt(:,p%DOFs%PIE(I) ) = RtHSdat%PMomBNcRt(:,p%DOFs%PIE(I) ) + TmpVec + + ENDDO ! I - All active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the tail boom center of mass (point I) + DO I = 1,p%DOFs%NPUE ! Loop through all active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the nacelle center of mass (point U) + + TmpVec1 = -p%NacMass*RtHSdat%PLinVelEU(p%DOFs%PUE(I),0,:) ! The portion of PFrcONcRt associated with the NacMass + TmpVec2 = CROSS_PRODUCT( RtHSdat%rOU, TmpVec1 ) ! The portion of PMomBNcRt associated with the NacMass + + RtHSdat%PFrcONcRt(:,p%DOFs%PUE(I) ) = RtHSdat%PFrcONcRt(:,p%DOFs%PUE(I) ) + TmpVec1 + RtHSdat%PMomBNcRt(:,p%DOFs%PUE(I) ) = RtHSdat%PMomBNcRt(:,p%DOFs%PUE(I) ) + TmpVec2 - p%Nacd2Iner*CoordSys%d2* & + DOT_PRODUCT( CoordSys%d2, RtHSdat%PAngVelEN(p%DOFs%PUE(I),0,:) ) + + ENDDO ! I - All active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the nacelle center of mass (point U) + + +!..................................... +! FrcONcRtt and MomBNcRtt +! (requires FrcVGnRtt, MomNGnRtt, FrcWTailt, MomNTailt) +!..................................... + + TmpVec1 = -p%NacMass*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccEUt ) ! The portion of FrcONcRtt associated with the NacMass + TmpVec2 = CROSS_PRODUCT( RtHSdat%rOU, TmpVec1 ) ! The portion of MomBNcRtt associated with the NacMass + TmpVec3 = CROSS_PRODUCT( RtHSdat%rOV, RtHSdat%FrcVGnRtt ) ! The portion of MomBNcRtt associated with the FrcVGnRtt + TmpVec4 = CROSS_PRODUCT( RtHSdat%rOW, RtHSdat%FrcWTailt ) ! The portion of MomBNcRtt associated with the FrcWTailt + TmpVec = p%Nacd2Iner*CoordSys%d2*DOT_PRODUCT( CoordSys%d2, RtHSdat%AngVelEN ) ! = ( Nacelle inertia dyadic ) dot ( angular velocity of nacelle in the inertia frame ) + + RtHSdat%FrcONcRtt = RtHSdat%FrcVGnRtt + RtHSdat%FrcWTailt + TmpVec1 + (/ u%NacelleLoads%Force(1,1), u%NacelleLoads%Force(3,1), -u%NacelleLoads%Force(2,1) /) + + RtHSdat%MomBNcRtt = RtHSdat%MomNGnRtt + RtHSdat%MomNTailt + TmpVec2 + TmpVec3 + TmpVec4 & + + CROSS_PRODUCT( -RtHSdat%AngVelEN, TmpVec ) & ! = ( -angular velocity of nacelle in the inertia frame ) cross ( TmpVec ) & + - p%Nacd2Iner*CoordSys%d2*DOT_PRODUCT( CoordSys%d2, RtHSdat%AngAccENt ) & + + (/ u%NacelleLoads%Moment(1,1), u%NacelleLoads%Moment(3,1), -u%NacelleLoads%Moment(2,1) /) + +!..................................... +! PFTHydro and PMFHydro +! (requires TwrAddedMass) +!..................................... + + ! Compute the partial hydrodynamic forces and moments per unit length + ! (including those associated with the QD2T()'s and those that are not) at the current tower element (point T) / (body F): + + ! NOTE: These forces are named PFTHydro, PMFHydro, FTHydrot, and MFHydrot. However, the names should not imply that the + ! forces are a result of hydrodynamic contributions only. These tower forces contain contributions from any external + ! load acting on the tower other than loads transmitted from aerodynamics. For example, these tower forces contain + ! contributions from foundation stiffness and damping [not floating] or mooring line restoring and damping, + ! as well as hydrostatic and hydrodynamic contributions [offshore]. + + DO J=1,p%TwrNodes + + RtHSdat%PFTHydro(:,J,:) = 0.0 + RtHSdat%PMFHydro(:,J,:) = 0.0 + DO I = 1,p%DOFs%NPTE ! Loop through all active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the tower + + RtHSdat%PFTHydro(:,J,p%DOFs%PTE(I)) = & + CoordSys%z1*( - u%TwrAddedMass(DOF_Sg,DOF_Sg,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,1) & + + u%TwrAddedMass(DOF_Sg,DOF_Sw,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,3) & + - u%TwrAddedMass(DOF_Sg,DOF_Hv,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,2) & + - u%TwrAddedMass(DOF_Sg,DOF_R ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,1) & + + u%TwrAddedMass(DOF_Sg,DOF_P ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,3) & + - u%TwrAddedMass(DOF_Sg,DOF_Y ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,2) ) & + - CoordSys%z3*( - u%TwrAddedMass(DOF_Sw,DOF_Sg,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,1) & + + u%TwrAddedMass(DOF_Sw,DOF_Sw,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,3) & + - u%TwrAddedMass(DOF_Sw,DOF_Hv,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,2) & + - u%TwrAddedMass(DOF_Sw,DOF_R ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,1) & + + u%TwrAddedMass(DOF_Sw,DOF_P ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,3) & + - u%TwrAddedMass(DOF_Sw,DOF_Y ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,2) ) & + + CoordSys%z2*( - u%TwrAddedMass(DOF_Hv,DOF_Sg,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,1) & + + u%TwrAddedMass(DOF_Hv,DOF_Sw,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,3) & + - u%TwrAddedMass(DOF_Hv,DOF_Hv,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,2) & + - u%TwrAddedMass(DOF_Hv,DOF_R ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,1) & + + u%TwrAddedMass(DOF_Hv,DOF_P ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,3) & + - u%TwrAddedMass(DOF_Hv,DOF_Y ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,2) ) + RtHSdat%PMFHydro(:,J,p%DOFs%PTE(I)) = & + CoordSys%z1*( - u%TwrAddedMass(DOF_R ,DOF_Sg,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,1) & + + u%TwrAddedMass(DOF_R ,DOF_Sw,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,3) & + - u%TwrAddedMass(DOF_R ,DOF_Hv,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,2) & + - u%TwrAddedMass(DOF_R ,DOF_R ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,1) & + + u%TwrAddedMass(DOF_R ,DOF_P ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,3) & + - u%TwrAddedMass(DOF_R ,DOF_Y ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,2) ) & + - CoordSys%z3*( - u%TwrAddedMass(DOF_P ,DOF_Sg,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,1) & + + u%TwrAddedMass(DOF_P ,DOF_Sw,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,3) & + - u%TwrAddedMass(DOF_P ,DOF_Hv,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,2) & + - u%TwrAddedMass(DOF_P ,DOF_R ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,1) & + + u%TwrAddedMass(DOF_P ,DOF_P ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,3) & + - u%TwrAddedMass(DOF_P ,DOF_Y ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,2) ) & + + CoordSys%z2*( - u%TwrAddedMass(DOF_Y ,DOF_Sg,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,1) & + + u%TwrAddedMass(DOF_Y ,DOF_Sw,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,3) & + - u%TwrAddedMass(DOF_Y ,DOF_Hv,J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,2) & + - u%TwrAddedMass(DOF_Y ,DOF_R ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,1) & + + u%TwrAddedMass(DOF_Y ,DOF_P ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,3) & + - u%TwrAddedMass(DOF_Y ,DOF_Y ,J)*RtHSdat%PAngVelEF(J,p%DOFs%PTE(I),0,2) ) + + END DO ! I - All active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the tower + END DO !J + +!..................................... +! FTHydrot and MFHydrot +! (requires TwrAddedMass) +!..................................... + + DO J=1,p%TwrNodes + RtHSdat%FTHydrot(:,J) = CoordSys%z1*( u%TowerPtLoads%Force(DOF_Sg,J)/abs(p%DHNodes(J)) & + - u%TwrAddedMass(DOF_Sg,DOF_Sg,J)*RtHSdat%LinAccETt(1,J) & + + u%TwrAddedMass(DOF_Sg,DOF_Sw,J)*RtHSdat%LinAccETt(3,J) & + - u%TwrAddedMass(DOF_Sg,DOF_Hv,J)*RtHSdat%LinAccETt(2,J) & + - u%TwrAddedMass(DOF_Sg,DOF_R ,J)*RtHSdat%AngAccEFt(1,J) & + + u%TwrAddedMass(DOF_Sg,DOF_P ,J)*RtHSdat%AngAccEFt(3,J) & + - u%TwrAddedMass(DOF_Sg,DOF_Y ,J)*RtHSdat%AngAccEFt(2,J) ) & + - CoordSys%z3*( u%TowerPtLoads%Force(DOF_Sw,J)/abs(p%DHNodes(J)) & + - u%TwrAddedMass(DOF_Sw,DOF_Sg,J)*RtHSdat%LinAccETt(1,J) & + + u%TwrAddedMass(DOF_Sw,DOF_Sw,J)*RtHSdat%LinAccETt(3,J) & + - u%TwrAddedMass(DOF_Sw,DOF_Hv,J)*RtHSdat%LinAccETt(2,J) & + - u%TwrAddedMass(DOF_Sw,DOF_R ,J)*RtHSdat%AngAccEFt(1,J) & + + u%TwrAddedMass(DOF_Sw,DOF_P ,J)*RtHSdat%AngAccEFt(3,J) & + - u%TwrAddedMass(DOF_Sw,DOF_Y ,J)*RtHSdat%AngAccEFt(2,J) ) & + + CoordSys%z2*( u%TowerPtLoads%Force(DOF_Hv,J)/abs(p%DHNodes(J)) & + - u%TwrAddedMass(DOF_Hv,DOF_Sg,J)*RtHSdat%LinAccETt(1,J) & + + u%TwrAddedMass(DOF_Hv,DOF_Sw,J)*RtHSdat%LinAccETt(3,J) & + - u%TwrAddedMass(DOF_Hv,DOF_Hv,J)*RtHSdat%LinAccETt(2,J) & + - u%TwrAddedMass(DOF_Hv,DOF_R ,J)*RtHSdat%AngAccEFt(1,J) & + + u%TwrAddedMass(DOF_Hv,DOF_P ,J)*RtHSdat%AngAccEFt(3,J) & + - u%TwrAddedMass(DOF_Hv,DOF_Y ,J)*RtHSdat%AngAccEFt(2,J) ) + RtHSdat%MFHydrot(:,J) = CoordSys%z1*( u%TowerPtLoads%Moment(DOF_R-3,J)/abs(p%DHNodes(J)) & + - u%TwrAddedMass(DOF_R ,DOF_Sg,J)*RtHSdat%LinAccETt(1,J) & + + u%TwrAddedMass(DOF_R ,DOF_Sw,J)*RtHSdat%LinAccETt(3,J) & + - u%TwrAddedMass(DOF_R ,DOF_Hv,J)*RtHSdat%LinAccETt(2,J) & + - u%TwrAddedMass(DOF_R ,DOF_R ,J)*RtHSdat%AngAccEFt(1,J) & + + u%TwrAddedMass(DOF_R ,DOF_P ,J)*RtHSdat%AngAccEFt(3,J) & + - u%TwrAddedMass(DOF_R ,DOF_Y ,J)*RtHSdat%AngAccEFt(2,J) ) & + - CoordSys%z3*( u%TowerPtLoads%Moment(DOF_P-3 ,J)/abs(p%DHNodes(J)) & + - u%TwrAddedMass(DOF_P ,DOF_Sg,J)*RtHSdat%LinAccETt(1,J) & + + u%TwrAddedMass(DOF_P ,DOF_Sw,J)*RtHSdat%LinAccETt(3,J) & + - u%TwrAddedMass(DOF_P ,DOF_Hv,J)*RtHSdat%LinAccETt(2,J) & + - u%TwrAddedMass(DOF_P ,DOF_R ,J)*RtHSdat%AngAccEFt(1,J) & + + u%TwrAddedMass(DOF_P ,DOF_P ,J)*RtHSdat%AngAccEFt(3,J) & + - u%TwrAddedMass(DOF_P ,DOF_Y ,J)*RtHSdat%AngAccEFt(2,J) ) & + + CoordSys%z2*( u%TowerPtLoads%Moment(DOF_Y-3 ,J)/abs(p%DHNodes(J)) & + - u%TwrAddedMass(DOF_Y ,DOF_Sg,J)*RtHSdat%LinAccETt(1,J) & + + u%TwrAddedMass(DOF_Y ,DOF_Sw,J)*RtHSdat%LinAccETt(3,J) & + - u%TwrAddedMass(DOF_Y ,DOF_Hv,J)*RtHSdat%LinAccETt(2,J) & + - u%TwrAddedMass(DOF_Y ,DOF_R ,J)*RtHSdat%AngAccEFt(1,J) & + + u%TwrAddedMass(DOF_Y ,DOF_P ,J)*RtHSdat%AngAccEFt(3,J) & + - u%TwrAddedMass(DOF_Y ,DOF_Y ,J)*RtHSdat%AngAccEFt(2,J) ) + + END DO !J + +!..................................... +! PFrcT0Trb and PMomX0Trb +! (requires PFrcONcRt, PMomBNcRt, PFrcT0Trb, PMomX0Trb, PFTHydro, PMFHydro) +!..................................... + + ! Initialize the partial forces and moments (including those associated + ! with the QD2T()'s and those that are not) at the tower base (point T(0)) using everything but the tower: + + RtHSdat%PFrcT0Trb = RtHSdat%PFrcONcRt ! Initialize these partial forces and moments + RtHSdat%PMomX0Trb = RtHSdat%PMomBNcRt ! using all of the effects above the yaw bearing + DO I = 1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs + + TmpVec = CROSS_PRODUCT( RtHSdat%rT0O, RtHSdat%PFrcONcRt(:,p%DOFs%SrtPS(I)) ) ! The portion of PMomX0Trb associated with the PFrcONcRt + + RtHSdat%PMomX0Trb(:,p%DOFs%SrtPS(I)) = RtHSdat%PMomX0Trb(:,p%DOFs%SrtPS(I)) + TmpVec + + ENDDO ! I - All active (enabled) DOFs + DO I = 1,p%DOFs%NPTE ! Loop through all active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the yaw bearing center of mass (point O) + + TmpVec1 = -p%YawBrMass*RtHSdat%PLinVelEO(p%DOFs%PTE(I),0,:) ! The portion of PFrcT0Trb associated with the YawBrMass + TmpVec2 = CROSS_PRODUCT( RtHSdat%rT0O, TmpVec1 ) ! The portion of PMomX0Trb associated with the YawBrMass + + RtHSdat%PFrcT0Trb(:,p%DOFs%PTE(I) ) = RtHSdat%PFrcT0Trb(:,p%DOFs%PTE(I) ) + TmpVec1 + RtHSdat%PMomX0Trb(:,p%DOFs%PTE(I) ) = RtHSdat%PMomX0Trb(:,p%DOFs%PTE(I) ) + TmpVec2 + + ENDDO ! I - All active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the yaw bearing center of mass (point O) + + TmpVec1 = -p%YawBrMass*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccEOt ) ! The portion of FrcT0Trbt associated with the YawBrMass + TmpVec2 = CROSS_PRODUCT( RtHSdat%rT0O, TmpVec1 ) ! The portion of MomX0Trbt associated with the YawBrMass + TmpVec3 = CROSS_PRODUCT( RtHSdat%rT0O, RtHSdat%FrcONcRtt ) ! The portion of MomX0Trbt associated with the FrcONcRtt + + RtHSdat%FrcT0Trbt = RtHSdat%FrcONcRtt + TmpVec1 + RtHSdat%MomX0Trbt = RtHSdat%MomBNcRtt + TmpVec2 + TmpVec3 + + ! Integrate to find the total partial forces and moments (including those + ! associated with the QD2T()'s and those that are not) at the tower base (point T(0)): + DO J=1,p%TwrNodes + + DO I = 1,p%DOFs%NPTE ! Loop through all active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the tower + + TmpVec1 = RtHSdat%PFTHydro(:,J,p%DOFs%PTE(I))*abs(p%DHNodes(J)) - p%TElmntMass(J)*RtHSdat%PLinVelET(J,p%DOFs%PTE(I),0,:) ! The portion of PFrcT0Trb associated with tower element J + TmpVec2 = CROSS_PRODUCT( RtHSdat%rT0T(:,J), TmpVec1 ) ! The portion of PMomX0Trb associated with tower element J + TmpVec3 = RtHSdat%PMFHydro(:,J,p%DOFs%PTE(I))*abs(p%DHNodes(J)) ! The added moment applied at tower element J + + RtHSdat%PFrcT0Trb(:,p%DOFs%PTE(I)) = RtHSdat%PFrcT0Trb(:,p%DOFs%PTE(I)) + TmpVec1 + RtHSdat%PMomX0Trb(:,p%DOFs%PTE(I)) = RtHSdat%PMomX0Trb(:,p%DOFs%PTE(I)) + TmpVec2 + TmpVec3 + + ENDDO ! I - All active (enabled) 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 + TmpVec2 = CROSS_PRODUCT( RtHSdat%rT0T(:,J), TmpVec1 ) ! The portion of MomX0Trbt associated with tower element J + TmpVec3 = ( RtHSdat%MFHydrot(:,J) )*abs(p%DHNodes(J)) ! The external moment applied to tower element J + + RtHSdat%FrcT0Trbt = RtHSdat%FrcT0Trbt + TmpVec1 + + RtHSdat%MomX0Trbt = RtHSdat%MomX0Trbt + TmpVec2 + TmpVec3 + + END DO !J + +!..................................... +! PFZHydro and PMXHydro +! ( requires PtfmAddedMass ) +!..................................... + + !.................................................................................................................................. + ! Compute the partial platform forces and moments (including those associated with the QD2T()'s and those that are not) at the + ! platform reference (point Z) / (body X). + ! + ! NOTE: These forces are named PFZHydro, PMXHydro, FZHydrot, and MXHydrot. However, the names should not imply that the forces + ! are a result of hydrodynamic contributions only. These platform forces contain contributions from any external load acting + ! on the platform other than loads transmitted from the wind turbine. For example, these platform forces contain contributions + ! from foundation stiffness and damping [not floating] or mooring line restoring and damping [floating], as well as hydrostatic + ! and hydrodynamic contributions [offshore]. + !bjj: m%RtHS%PFZHydro, %PMXHydro, %FZHydrot, and %MXHydrot are not used in the output routine anymore + ! (because of their dependence on inputs, u) + + RtHSdat%PFZHydro = 0.0 + RtHSdat%PMXHydro = 0.0 + DO I = 1,p%DOFs%NPYE ! Loop through all active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the platform center of mass (point Y) + + RtHSdat%PFZHydro(p%DOFs%PYE(I),:) = - u%PtfmAddedMass(DOF_Sg,p%DOFs%PYE(I))*RtHSdat%PLinVelEZ(DOF_Sg,0,:) & + - u%PtfmAddedMass(DOF_Sw,p%DOFs%PYE(I))*RtHSdat%PLinVelEZ(DOF_Sw,0,:) & + - u%PtfmAddedMass(DOF_Hv,p%DOFs%PYE(I))*RtHSdat%PLinVelEZ(DOF_Hv,0,:) + RtHSdat%PMXHydro(p%DOFs%PYE(I),:) = - u%PtfmAddedMass(DOF_R ,p%DOFs%PYE(I))*RtHSdat%PAngVelEX(DOF_R ,0,:) & + - u%PtfmAddedMass(DOF_P ,p%DOFs%PYE(I))*RtHSdat%PAngVelEX(DOF_P ,0,:) & + - u%PtfmAddedMass(DOF_Y ,p%DOFs%PYE(I))*RtHSdat%PAngVelEX(DOF_Y ,0,:) + + ENDDO ! I - All active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the platform center of mass (point Y) + + RtHSdat%FZHydrot = u%PlatformPtMesh%Force(DOF_Sg,1)*RtHSdat%PLinVelEZ(DOF_Sg,0,:) & + + u%PlatformPtMesh%Force(DOF_Sw,1)*RtHSdat%PLinVelEZ(DOF_Sw,0,:) & + + u%PlatformPtMesh%Force(DOF_Hv,1)*RtHSdat%PLinVelEZ(DOF_Hv,0,:) + RtHSdat%MXHydrot = u%PlatformPtMesh%Moment(DOF_R-3,1)*RtHSdat%PAngVelEX(DOF_R ,0,:) & + + u%PlatformPtMesh%Moment(DOF_P-3,1)*RtHSdat%PAngVelEX(DOF_P ,0,:) & + + u%PlatformPtMesh%Moment(DOF_Y-3,1)*RtHSdat%PAngVelEX(DOF_Y ,0,:) + +!..................................... +! PFrcZAll and PMomXAll +! (requires PFrcT0Trb, PMomX0Trb, PFZHydro, PMXHydro ) +!..................................... + + ! Define the partial forces and moments (including those associated with the QD2T()'s and those that are not) at the + ! platform reference (point Z) / (body X) using the turbine and platform effects: + + RtHSdat%PFrcZAll = RtHSdat%PFrcT0Trb ! Initialize these partial forces and moments + RtHSdat%PMomXAll = RtHSdat%PMomX0Trb ! using the effects from the wind turbine + DO I = 1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs + + TmpVec = CROSS_PRODUCT( RtHSdat%rZT0, RtHSdat%PFrcT0Trb(:,p%DOFs%SrtPS(I)) ) ! The portion of PMomXAll associated with the PFrcT0Trb + + RtHSdat%PMomXAll(:,p%DOFs%SrtPS(I)) = RtHSdat%PMomXAll(:,p%DOFs%SrtPS(I)) + TmpVec + + ENDDO ! I - All active (enabled) DOFs + DO I = 1,p%DOFs%NPYE ! Loop through all active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the platform center of mass (point Y) + + TmpVec1 = -p%PtfmMass*RtHSdat%PLinVelEY(p%DOFs%PYE(I),0,:) ! The portion of PFrcZAll associated with the PtfmMass + TmpVec2 = CROSS_PRODUCT( RtHSdat%rZY , TmpVec1 ) ! The portion of PMomXAll associated with the PtfmMass + + RtHSdat%PFrcZAll(:,p%DOFs%PYE(I)) = RtHSdat%PFrcZAll(:,p%DOFs%PYE(I) ) + RtHSdat%PFZHydro(p%DOFs%PYE(I),:) + TmpVec1 + RtHSdat%PMomXAll(:,p%DOFs%PYE(I)) = RtHSdat%PMomXAll(:,p%DOFs%PYE(I) ) + RtHSdat%PMXHydro(p%DOFs%PYE(I),:) + TmpVec2 & + - p%PtfmRIner*CoordSys%a1*DOT_PRODUCT( CoordSys%a1, RtHSdat%PAngVelEX(p%DOFs%PYE(I),0,:) ) & + - p%PtfmYIner*CoordSys%a2*DOT_PRODUCT( CoordSys%a2, RtHSdat%PAngVelEX(p%DOFs%PYE(I),0,:) ) & + - p%PtfmPIner*CoordSys%a3*DOT_PRODUCT( CoordSys%a3, RtHSdat%PAngVelEX(p%DOFs%PYE(I),0,:) ) + + ENDDO ! I - All active (enabled) DOFs that contribute to the QD2T-related linear accelerations of the platform center of mass (point Y) + +!..................................... +! FrcZAllt and MomXAllt +! (requires FrcT0Trbt, MomX0Trbt) +!..................................... + + TmpVec1 = -p%PtfmMass*( p%Gravity*CoordSys%z2 + RtHSdat%LinAccEYt ) ! The portion of FrcZAllt associated with the PtfmMass + TmpVec2 = CROSS_PRODUCT( RtHSdat%rZY , TmpVec1 ) ! The portion of MomXAllt associated with the PtfmMass + TmpVec3 = CROSS_PRODUCT( RtHSdat%rZT0 , RtHSdat%FrcT0Trbt ) ! The portion of MomXAllt associated with the FrcT0Trbt + TmpVec = p%PtfmRIner*CoordSys%a1*DOT_PRODUCT( CoordSys%a1, RtHSdat%AngVelEX ) & ! = ( Platform inertia dyadic ) dot ( angular velocity of platform in the inertia frame ) + + p%PtfmYIner*CoordSys%a2*DOT_PRODUCT( CoordSys%a2, RtHSdat%AngVelEX ) & + + p%PtfmPIner*CoordSys%a3*DOT_PRODUCT( CoordSys%a3, RtHSdat%AngVelEX ) + TmpVec4 = CROSS_PRODUCT( -RtHSdat%AngVelEX, TmpVec ) ! = ( -angular velocity of platform in the inertia frame ) cross ( TmpVec ) + + RtHSdat%FrcZAllt = RtHSdat%FrcT0Trbt + RtHSdat%FZHydrot + TmpVec1 + RtHSdat%MomXAllt = RtHSdat%MomX0Trbt + RtHSdat%MXHydrot + TmpVec2 + TmpVec3 + TmpVec4 + + +END SUBROUTINE CalculateForcesMoments +``` \ No newline at end of file diff --git a/书籍/力学书籍/OpenFast/FASTKinetics/auto/CalculatePositions.md b/书籍/力学书籍/OpenFast/FASTKinetics/auto/CalculatePositions.md new file mode 100644 index 0000000..0a550b5 --- /dev/null +++ b/书籍/力学书籍/OpenFast/FASTKinetics/auto/CalculatePositions.md @@ -0,0 +1,143 @@ + +```fortran +SUBROUTINE CalculatePositions( p, x, CoordSys, RtHSdat ) +!.................................................................................................................................. + + ! 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_RtHndSide), INTENT(INOUT) :: RtHSdat !< data from the RtHndSid module (contains positions to be set) + + !Local variables + !REAL(R8Ki) :: rQ (3) ! Position vector from inertial frame origin to apex of rotation (point Q). + + INTEGER(IntKi) :: J ! Counter for elements + INTEGER(IntKi) :: K ! Counter for blades + + !------------------------------------------------------------------------------------------------- + ! Positions + !------------------------------------------------------------------------------------------------- + + ! Define the position vectors between the various points on the wind turbine + ! that are not dependent on the distributed tower or blade parameters: + + RtHSdat%rZ = x%QT(DOF_Sg)* CoordSys%z1 + x%QT(DOF_Hv)* CoordSys%z2 - x%QT(DOF_Sw)* CoordSys%z3 ! Position vector from inertia frame origin to platform reference (point Z). + RtHSdat%rZY = p%rZYzt* CoordSys%a2 + p%PtfmCMxt*CoordSys%a1 - p%PtfmCMyt*CoordSys%a3 ! Position vector from platform reference (point Z) to platform mass center (point Y). + RtHSdat%rZT0 = p%rZT0zt* CoordSys%a2 ! Position vector from platform reference (point Z) to tower base (point T(0)) + RtHSdat%rZO = ( x%QT(DOF_TFA1) + x%QT(DOF_TFA2) )*CoordSys%a1 & ! Position vector from platform reference (point Z) to tower-top / base plate (point O). + + ( p%RefTwrHt - 0.5*( p%AxRedTFA(1,1,p%TTopNode)*x%QT(DOF_TFA1)*x%QT(DOF_TFA1) & + + p%AxRedTFA(2,2,p%TTopNode)*x%QT(DOF_TFA2)*x%QT(DOF_TFA2) & + + 2.0*p%AxRedTFA(1,2,p%TTopNode)*x%QT(DOF_TFA1)*x%QT(DOF_TFA2) & + + p%AxRedTSS(1,1,p%TTopNode)*x%QT(DOF_TSS1)*x%QT(DOF_TSS1) & + + p%AxRedTSS(2,2,p%TTopNode)*x%QT(DOF_TSS2)*x%QT(DOF_TSS2) & + + 2.0*p%AxRedTSS(1,2,p%TTopNode)*x%QT(DOF_TSS1)*x%QT(DOF_TSS2) ) )*CoordSys%a2 & + + ( x%QT(DOF_TSS1) + x%QT(DOF_TSS2) )*CoordSys%a3 + RtHSdat%rOU = p%NacCMxn*CoordSys%d1 + p%NacCMzn *CoordSys%d2 - p%NacCMyn *CoordSys%d3 ! Position vector from tower-top / base plate (point O) to nacelle center of mass (point U). + RtHSdat%rOV = p%RFrlPnt_n(1)*CoordSys%d1 + p%RFrlPnt_n(3)*CoordSys%d2 - p%RFrlPnt_n(2)*CoordSys%d3 ! Position vector from tower-top / base plate (point O) to specified point on rotor-furl axis (point V). + RtHSdat%rVIMU = p%rVIMUxn*CoordSys%rf1 + p%rVIMUzn *CoordSys%rf2 - p%rVIMUyn *CoordSys%rf3 ! Position vector from specified point on rotor-furl axis (point V) to nacelle IMU (point IMU). + RtHSdat%rVD = p%rVDxn*CoordSys%rf1 + p%rVDzn *CoordSys%rf2 - p%rVDyn *CoordSys%rf3 ! Position vector from specified point on rotor-furl axis (point V) to center of mass of structure that furls with the rotor (not including rotor) (point D). + RtHSdat%rVP = p%rVPxn*CoordSys%rf1 + p%rVPzn *CoordSys%rf2 - p%rVPyn *CoordSys%rf3 + p%OverHang*CoordSys%c1 ! Position vector from specified point on rotor-furl axis (point V) to teeter pin (point P). + RtHSdat%rPQ = -p%UndSling*CoordSys%g1 ! Position vector from teeter pin (point P) to apex of rotation (point Q). + RtHSdat%rQC = p%HubCM*CoordSys%g1 ! Position vector from apex of rotation (point Q) to hub center of mass (point C). + RtHSdat%rOW = p%TFrlPnt_n(1)*CoordSys%d1 + p%TFrlPnt_n(3) *CoordSys%d2 - p%TFrlPnt_n(2)*CoordSys%d3 ! Position vector from tower-top / base plate (point O) to specified point on tail-furl axis (point W). + RtHSdat%rWI = p%rWIxn*CoordSys%tf1 + p%rWIzn*CoordSys%tf2 - p%rWIyn*CoordSys%tf3 ! Position vector from specified point on tail-furl axis (point W) to tail boom center of mass (point I). + RtHSdat%rWJ = p%rWJxn*CoordSys%tf1 + p%rWJzn*CoordSys%tf2 - p%rWJyn*CoordSys%tf3 ! Position vector from specified point on tail-furl axis (point W) to tail fin center of mass (point J). + RtHSdat%rPC = RtHSdat%rPQ + RtHSdat%rQC ! Position vector from teeter pin (point P) to hub center of mass (point C). + RtHSdat%rT0O = RtHSdat%rZO - RtHSdat%rZT0 ! Position vector from the tower base (point T(0)) to tower-top / base plate (point O). + RtHSdat%rO = RtHSdat%rZ + RtHSdat%rZO ! Position vector from inertial frame origin to tower-top / base plate (point O). + RtHSdat%rV = RtHSdat%rO + RtHSdat%rOV ! Position vector from inertial frame origin to specified point on rotor-furl axis (point V) + !RtHSdat%rP = RtHSdat%rO + RtHSdat%rOV + RtHSdat%rVP ! Position vector from inertial frame origin to teeter pin (point P). + RtHSdat%rP = RtHSdat%rV + RtHSdat%rVP ! Position vector from inertial frame origin to teeter pin (point P). + RtHSdat%rQ = RtHSdat%rP + RtHSdat%rPQ ! Position vector from inertial frame origin to apex of rotation (point Q). + RtHSdat%rJ = RtHSdat%rO + RtHSdat%rOW + RtHSdat%rWJ ! Position vector from inertial frame origin to tail fin center of mass (point J). + + + DO K = 1,p%NumBl ! Loop through all blades + + ! Calculate the position vector of the tip: + RtHSdat%rS0S(:,K,p%TipNode) = ( p%TwistedSF(K,1,1,p%TipNode,0)*x%QT( DOF_BF(K,1) ) & ! Position vector from the blade root (point S(0)) to the blade tip (point S(p%BldFlexL)). + + p%TwistedSF(K,1,2,p%TipNode,0)*x%QT( DOF_BF(K,2) ) & + + p%TwistedSF(K,1,3,p%TipNode,0)*x%QT( DOF_BE(K,1) ) )*CoordSys%j1(K,:) & + + ( p%TwistedSF(K,2,1,p%TipNode,0)*x%QT( DOF_BF(K,1) ) & + + p%TwistedSF(K,2,2,p%TipNode,0)*x%QT( DOF_BF(K,2) ) & + + p%TwistedSF(K,2,3,p%TipNode,0)*x%QT( DOF_BE(K,1) ) )*CoordSys%j2(K,:) & + + ( p%BldFlexL - 0.5* & + ( p%AxRedBld(K,1,1,p%TipNode)*x%QT( DOF_BF(K,1) )*x%QT( DOF_BF(K,1) ) & + + p%AxRedBld(K,2,2,p%TipNode)*x%QT( DOF_BF(K,2) )*x%QT( DOF_BF(K,2) ) & + + p%AxRedBld(K,3,3,p%TipNode)*x%QT( DOF_BE(K,1) )*x%QT( DOF_BE(K,1) ) & + + 2.*p%AxRedBld(K,1,2,p%TipNode)*x%QT( DOF_BF(K,1) )*x%QT( DOF_BF(K,2) ) & + + 2.*p%AxRedBld(K,2,3,p%TipNode)*x%QT( DOF_BF(K,2) )*x%QT( DOF_BE(K,1) ) & + + 2.*p%AxRedBld(K,1,3,p%TipNode)*x%QT( DOF_BF(K,1) )*x%QT( DOF_BE(K,1) ) ) )*CoordSys%j3(K,:) + RtHSdat%rQS (:,K,p%TipNode) = RtHSdat%rS0S(:,K,p%TipNode) + p%HubRad*CoordSys%j3(K,:) ! Position vector from apex of rotation (point Q) to the blade tip (point S(p%BldFlexL)). + RtHSdat%rS (:,K,p%TipNode) = RtHSdat%rQS (:,K,p%TipNode) + RtHSdat%rQ ! Position vector from inertial frame origin to the blade tip (point S(p%BldFlexL)). + + ! position vectors for blade root node: + RtHSdat%rQS (:,K,0) = p%HubRad*CoordSys%j3(K,:) + RtHSdat%rS (:,K,0) = p%HubRad*CoordSys%j3(K,:) + RtHSdat%rQ + + + ! Calculate the position vector from the teeter pin to the blade root: + + RtHSdat%rPS0(:,K) = RtHSdat%rPQ + p%HubRad*CoordSys%j3(K,:) ! Position vector from teeter pin (point P) to blade root (point S(0)). + + + DO J = 1,p%BldNodes ! Loop through the blade nodes / elements + + + ! Calculate the position vector of the current node: + + RtHSdat%rS0S(:,K,J) = ( p%TwistedSF(K,1,1,J,0)*x%QT( DOF_BF(K,1) ) & ! Position vector from the blade root (point S(0)) to the current node (point S(RNodes(J)). + + p%TwistedSF(K,1,2,J,0)*x%QT( DOF_BF(K,2) ) & + + p%TwistedSF(K,1,3,J,0)*x%QT( DOF_BE(K,1) ) )*CoordSys%j1(K,:) & + + ( p%TwistedSF(K,2,1,J,0)*x%QT( DOF_BF(K,1) ) & + + p%TwistedSF(K,2,2,J,0)*x%QT( DOF_BF(K,2) ) & + + p%TwistedSF(K,2,3,J,0)*x%QT( DOF_BE(K,1) ) )*CoordSys%j2(K,:) & + + ( p%RNodes(J) - 0.5* & + ( p%AxRedBld(K,1,1,J)*x%QT( DOF_BF(K,1) )*x%QT( DOF_BF(K,1) ) & + + p%AxRedBld(K,2,2,J)*x%QT( DOF_BF(K,2) )*x%QT( DOF_BF(K,2) ) & + + p%AxRedBld(K,3,3,J)*x%QT( DOF_BE(K,1) )*x%QT( DOF_BE(K,1) ) & + + 2.0*p%AxRedBld(K,1,2,J)*x%QT( DOF_BF(K,1) )*x%QT( DOF_BF(K,2) ) & + + 2.0*p%AxRedBld(K,2,3,J)*x%QT( DOF_BF(K,2) )*x%QT( DOF_BE(K,1) ) & + + 2.0*p%AxRedBld(K,1,3,J)*x%QT( DOF_BF(K,1) )*x%QT( DOF_BE(K,1) ) ) )*CoordSys%j3(K,:) + RtHSdat%rQS (:,K,J) = RtHSdat%rS0S(:,K,J) + p%HubRad*CoordSys%j3(K,:) ! Position vector from apex of rotation (point Q) to the current node (point S(RNodes(J)). + RtHSdat%rS (:,K,J) = RtHSdat%rQS (:,K,J) + RtHSdat%rQ ! Position vector from inertial frame origin to the current node (point S(RNodes(J)). + + + END DO !J = 1,p%BldNodes ! Loop through the blade nodes / elements + + + + END DO !K = 1,p%NumBl + + + + + !---------------------------------------------------------------------------------------------------- + ! Get the tower element positions + !---------------------------------------------------------------------------------------------------- + RtHSdat%rZT (:,0) = RtHSdat%rZT0 + DO J = 1,p%TwrNodes ! Loop through the tower nodes / elements + + + ! Calculate the position vector of the current node: + + RtHSdat%rT0T(:,J) = ( p%TwrFASF(1,J,0)*x%QT(DOF_TFA1) + p%TwrFASF(2,J,0)*x%QT(DOF_TFA2) )*CoordSys%a1 & ! Position vector from base of flexible portion of tower (point T(0)) to current node (point T(J)). + + ( p%HNodes(J) - 0.5*( p%AxRedTFA(1,1,J)*x%QT(DOF_TFA1)*x%QT(DOF_TFA1) & + + p%AxRedTFA(2,2,J)*x%QT(DOF_TFA2)*x%QT(DOF_TFA2) & + + 2.0*p%AxRedTFA(1,2,J)*x%QT(DOF_TFA1)*x%QT(DOF_TFA2) & + + p%AxRedTSS(1,1,J)*x%QT(DOF_TSS1)*x%QT(DOF_TSS1) & + + p%AxRedTSS(2,2,J)*x%QT(DOF_TSS2)*x%QT(DOF_TSS2) & + + 2.0*p%AxRedTSS(1,2,J)*x%QT(DOF_TSS1)*x%QT(DOF_TSS2) ) )*CoordSys%a2 & + + ( p%TwrSSSF(1,J,0)*x%QT(DOF_TSS1) + p%TwrSSSF(2,J,0)*x%QT(DOF_TSS2) )*CoordSys%a3 + RtHSdat%rZT (:,J) = RtHSdat%rZT0 + RtHSdat%rT0T(:,J) ! Position vector from platform reference (point Z) to the current node (point T(HNodes(J)). + + + RtHSdat%rT(:,J) = RtHSdat%rZ + RtHSdat%rZT (:,J) ! Position vector from inertial frame origin to the current node (point T(HNodes(J)). + + END DO + + +END SUBROUTINE CalculatePositions + +``` \ No newline at end of file diff --git a/书籍/力学书籍/OpenFast/FASTKinetics/auto/Coeff.md b/书籍/力学书籍/OpenFast/FASTKinetics/auto/Coeff.md new file mode 100644 index 0000000..e1e595a --- /dev/null +++ b/书籍/力学书籍/OpenFast/FASTKinetics/auto/Coeff.md @@ -0,0 +1,678 @@ + +```fortran +SUBROUTINE Coeff(p,InputFileData, ErrStat, ErrMsg) +!.................................................................................................................................. + + ! Passed variables + + TYPE(ED_ParameterType), INTENT(INOUT) :: p !< Parameters of the structural dynamics module + TYPE(ED_InputFile), INTENT(IN) :: InputFileData !< all the data in the ElastoDyn input file + INTEGER(IntKi), INTENT(OUT) :: ErrStat !< Error status + CHARACTER(*), INTENT(OUT) :: ErrMsg !< Error message when ErrStat =/ ErrID_None + + + ! Local variables. + + REAL(ReKi) :: AxRdBld (3,3) ! Temporary result holding the current addition to the p%AxRedBld() array. + REAL(ReKi) :: AxRdBldOld(3,3) ! Previous AxRdBld (i.e., AxRdBld from the previous node) + REAL(ReKi) :: AxRdTFA (2,2) ! Temporary result holding the current addition to the AxRedTFA() array. + REAL(ReKi) :: AxRdTFAOld(2,2) ! Previous AxRdTFA (i.e., AxRdTFA from the previous node) + REAL(ReKi) :: AxRdTSS (2,2) ! Temporary result holding the current addition to the AxRedTSS() array. + REAL(ReKi) :: AxRdTSSOld(2,2) ! Previous AxRdTSS (i.e., AxRdTSS from the previous node) + REAL(ReKi) :: TmpDist ! Temporary distance used in the calculation of the aero center locations. + REAL(ReKi) :: TmpDistj1 ! Temporary distance used in the calculation of the aero center locations. + REAL(ReKi) :: TmpDistj2 ! Temporary distance used in the calculation of the aero center locations. + REAL(ReKi) :: ElmntStff ! (Temporary) stiffness of an element. + REAL(ReKi) :: ElStffFA ! (Temporary) tower fore-aft stiffness of an element + REAL(ReKi) :: ElStffSS ! (Temporary) tower side-to-side stiffness of an element + REAL(ReKi) :: FMomAbvNd (p%NumBl,p%BldNodes) ! FMomAbvNd(K,J) = portion of the first moment of blade K about the rotor centerline (not root, like FirstMom(K)) associated with everything above node J (including tip brake masses). + REAL(ReKi) :: KBECent (p%NumBl,1,1) ! Centrifugal-term of generalized edgewise stiffness of the blades. + REAL(ReKi) :: KBFCent (p%NumBl,2,2) ! Centrifugal-term of generalized flapwise stiffness of the blades. + REAL(ReKi) :: KTFAGrav (2,2) ! Gravitational-term of generalized fore-aft stiffness of the tower. + REAL(ReKi) :: KTSSGrav (2,2) ! Gravitational-term of generalized side-to-side stiffness of the tower. + REAL(ReKi) :: MBE (p%NumBl,1,1) ! Generalized edgewise mass of the blades. + REAL(ReKi) :: MBF (p%NumBl,2,2) ! Generalized flapwise mass of the blades. + REAL(ReKi) :: MTFA (2,2) ! Generalized fore-aft mass of the tower. + REAL(ReKi) :: MTSS (2,2) ! Generalized side-to-side mass of the tower. + REAL(ReKi) :: Shape ! Temporary result holding a value from the SHP function + REAL(ReKi) :: Shape1 ! Temporary result holding a value from the SHP function + REAL(ReKi) :: Shape2 ! Temporary result holding a value from the SHP function + REAL(ReKi) :: TMssAbvNd (p%TwrNodes) ! Portion of the tower mass associated with everything above node J (including tower-top effects) + REAL(ReKi) :: TwstdSF (2,3,0:1) ! Temperory result holding the current addition to the TwistedSF() array. + REAL(ReKi) :: TwstdSFOld(2,3,0:1) ! Previous TwstdSF (i.e., TwstdSF from the previous node) + + INTEGER(IntKi) :: I ! Generic index. + INTEGER(IntKi) :: J ! Loops through nodes / elements. + INTEGER(IntKi) :: K ! Loops through blades. + INTEGER(IntKi) :: L ! Generic index + + + ErrStat = ErrID_None + ErrMsg = '' + + !............................................................................................................................... + ! Calculate the distances from point S on a blade to the aerodynamic center in the j1 and j2 directions: + !............................................................................................................................... + + DO K = 1,p%NumBl ! Loop through the blades + + DO J = 1,p%BldNodes ! Loop through the blade nodes / elements + + TmpDist = ( 0.25 - p%PitchAxis(K,J) )*p%Chord(J) ! Distance along the chordline from point S (25% chord) to the aerodynamic center of the blade element J--positive towards the trailing edge. + TmpDistj1 = TmpDist*p%SAeroTwst(J) ! Distance along the j1-axis from point S (25% chord) to the aerodynamic center of the blade element J + TmpDistj2 = TmpDist*p%CAeroTwst(J) ! Distance along the j2-axis from point S (25% chord) to the aerodynamic center of the blade element J + p%rSAerCenn1(K,J) = TmpDistj1*p%CThetaS(K,J) - TmpDistj2*p%SThetaS(K,J) + p%rSAerCenn2(K,J) = TmpDistj1*p%SThetaS(K,J) + TmpDistj2*p%CThetaS(K,J) + + ENDDO ! J - Blade nodes / elements + + ENDDO ! K - Blades + + + !............................................................................................................................... + ! Calculate the structure that furls with the rotor inertia term: + !............................................................................................................................... + + p%RrfaIner = InputFileData%RFrlIner - p%RFrlMass*( (p%rVDxn**2 )*( 1.0 - p%CRFrlSkw2*p%CRFrlTlt2 ) & + + (p%rVDzn**2 )* p%CRFrlTlt2 & + + (p%rVDyn**2 )*( 1.0 - p%SRFrlSkw2*p%CRFrlTlt2 ) & + - 2.0*p%rVDxn*p%rVDzn* p%CRFrlSkew*p%CSRFrlTlt & + - 2.0*p%rVDxn*p%rVDyn* p%CSRFrlSkw*p%CRFrlTlt2 & + - 2.0*p%rVDzn*p%rVDyn* p%SRFrlSkew*p%CSRFrlTlt ) + IF ( p%RrfaIner < 0.0 ) THEN + ErrStat = ErrID_Fatal + ErrMsg = ' RFrlIner must not be less than RFrlMass*( perpendicular distance between rotor-furl'// & + ' axis and CM of the structure that furls with the rotor [not including rotor] )^2.' + RETURN + END IF + + !............................................................................................................................... + ! Calculate the tail boom inertia term: + !............................................................................................................................... + + p%AtfaIner = p%TFrlIner - p%BoomMass*( p%rWIxn*p%rWIxn*( 1.0 - p%CTFrlSkw2*p%CTFrlTlt2 ) & + + p%rWIzn*p%rWIzn* p%CTFrlTlt2 & + + p%rWIyn*p%rWIyn*( 1.0 - p%STFrlSkw2*p%CTFrlTlt2 ) & + - 2.0*p%rWIxn*p%rWIzn* p%CTFrlSkew*p%CSTFrlTlt & + - 2.0*p%rWIxn*p%rWIyn* p%CSTFrlSkw*p%CTFrlTlt2 & + - 2.0*p%rWIzn*p%rWIyn* p%STFrlSkew*p%CSTFrlTlt ) + IF ( p%AtfaIner < 0.0 ) THEN + ErrStat = ErrID_Fatal + ErrMsg = ' TFrlIner must not be less than BoomMass*( perpendicular distance between tail-furl'// & + ' axis and tail boom CM )^2.' + RETURN + ENDIF + + !............................................................................................................................... + ! Calculate the nacelle inertia terms: + !............................................................................................................................... + + p%Nacd2Iner = InputFileData%NacYIner - p%NacMass*( p%NacCMxn**2 + p%NacCMyn**2 ) ! Nacelle inertia about the d2-axis + IF ( p%Nacd2Iner < 0.0 ) THEN + ErrStat = ErrID_Fatal + ErrMsg = ' NacYIner must not be less than NacMass*( NacCMxn^2 + NacCMyn^2 ).' + RETURN + END IF + + ! Calculate hub inertia about its centerline passing through its c.g.. + ! This calculation assumes that the hub for a 2-blader is essentially + ! a uniform cylinder whose centerline is transverse through the cylinder + ! passing through its c.g.. That is, for a 2-blader, Hubg1Iner = + ! Hubg2Iner is the inertia of the hub about both the g1- and g2- axes. For + ! 3-bladers, Hubg1Iner is simply equal to HubIner and Hubg2Iner is zero. + ! Also, Initialize RotMass and RotIner to associated hub properties: + + IF ( p%NumBl == 2 ) THEN ! 2-blader + p%Hubg1Iner = ( InputFileData%HubIner - p%HubMass*( ( p%UndSling - p%HubCM )**2 ) )/( p%CosDel3**2 ) + p%Hubg2Iner = p%Hubg1Iner + IF ( p%Hubg1Iner < 0.0 ) THEN + ErrStat = ErrID_Fatal + ErrMsg = ' HubIner must not be less than HubMass*( UndSling - HubCM )^2 for 2-blader.' + RETURN + END IF + ELSE ! 3-blader + p%Hubg1Iner = InputFileData%HubIner + p%Hubg2Iner = 0.0 + ENDIF + + p%RotMass = p%HubMass + p%RotIner = p%Hubg1Iner + + + !............................................................................................................................... + + ! Initialize several variables to 0.0: + + p%KBF = 0.0 + p%KBE = 0.0 + KBFCent = 0.0 + KBECent = 0.0 + + p%TwrMass = 0.0 + p%KTFA = 0.0 + p%KTSS = 0.0 + KTFAGrav = 0.0 + KTSSGrav = 0.0 + + + + DO K = 1,p%NumBl ! Loop through the blades + + + ! Initialize BldMass(), FirstMom(), and SecondMom() using TipMass() effects: + + p%BldMass (K) = p%TipMass(K) + p%FirstMom (K) = p%TipMass(K)*p%BldFlexL + p%SecondMom(K) = p%TipMass(K)*p%BldFlexL*p%BldFlexL + + + DO J = p%BldNodes,1,-1 ! Loop through the blade nodes / elements in reverse + + + ! Calculate the mass of the current element + + p%BElmntMass(J,K) = p%MassB(K,J)*p%DRNodes(J) ! Mass of blade element J + + + ! Integrate to find some blade properties which will be output in .fsm + + p%BldMass (K) = p%BldMass (K) + p%BElmntMass(J,K) + p%FirstMom (K) = p%FirstMom (K) + p%BElmntMass(J,K)*p%RNodes(J) + p%SecondMom(K) = p%SecondMom(K) + p%BElmntMass(J,K)*p%RNodes(J)*p%RNodes(J) + + + ! Integrate to find FMomAbvNd: + + FMomAbvNd (K,J) = ( 0.5*p%BElmntMass(J,K) )*( p%HubRad + p%RNodes(J ) + 0.5*p%DRNodes(J ) ) + + IF ( J == p%BldNodes ) THEN ! Outermost blade element + ! Add the TipMass() effects: + + FMomAbvNd(K,J) = FmomAbvNd(K,J) + p%TipMass(K)*p%TipRad + ELSE ! All other blade elements + ! Add to FMomAbvNd(K,J) the effects from the (not yet used) portion of element J+1 + + FMomAbvNd(K,J) = FMomAbvNd(K,J) + FMomAbvNd(K,J+1) & + + ( 0.5*p%BElmntMass(J+1,K) )*( p%HubRad + p%RNodes(J+1) - 0.5*p%DRNodes(J+1) ) + ENDIF + + + ENDDO ! J - Blade nodes / elements in reverse + + IF (.NOT. p%BD4Blades) THEN + ! Calculate BldCG() using FirstMom() and BldMass(); and calculate + ! RotMass and RotIner: + + p%BldCG (K) = p%FirstMom (K) / p%BldMass (K) + p%RotMass = p%RotMass + p%BldMass (K) + p%RotIner = p%RotIner + ( p%SecondMom(K) + p%BldMass (K)*p%HubRad*( 2.0*p%BldCG(K) + p%HubRad ) )*( p%CosPreC(K)**2 ) + END IF + + ENDDO ! K - Blades + + + + DO K = 1,p%NumBl ! Loop through the blades + + + ! Initialize the generalized blade masses using tip mass effects: + + MBF(K,1,1) = p%TipMass(K) + MBF(K,2,2) = p%TipMass(K) + MBE(K,1,1) = p%TipMass(K) + + + DO J = 1,p%BldNodes ! Loop through the blade nodes / elements + + + ! Integrate to find the generalized mass of the blade (including tip mass effects). + ! Ignore the cross-correlation terms of MBF (i.e. MBF(i,j) where i /= j) since + ! these terms will never be used. + + Shape1 = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl1Sh(:,K), 0, ErrStat, ErrMsg ) + Shape2 = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl2Sh(:,K), 0, ErrStat, ErrMsg ) + MBF (K,1,1) = MBF (K,1,1) + p%BElmntMass(J,K)*Shape1*Shape1 + MBF (K,2,2) = MBF (K,2,2) + p%BElmntMass(J,K)*Shape2*Shape2 + + Shape = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldEdgSh(:,K), 0, ErrStat, ErrMsg ) + MBE (K,1,1) = MBE (K,1,1) + p%BElmntMass(J,K)*Shape *Shape + + + ! Integrate to find the generalized stiffness of the blade (not including centrifugal + ! effects). + + ElmntStff = p%StiffBF(K,J)*p%DRNodes(J) ! Flapwise stiffness of blade element J + Shape1 = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl1Sh(:,K), 2, ErrStat, ErrMsg ) + Shape2 = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl2Sh(:,K), 2, ErrStat, ErrMsg ) + p%KBF (K,1,1) = p%KBF (K,1,1) + ElmntStff*Shape1*Shape1 + p%KBF (K,1,2) = p%KBF (K,1,2) + ElmntStff*Shape1*Shape2 + p%KBF (K,2,1) = p%KBF (K,2,1) + ElmntStff*Shape2*Shape1 + p%KBF (K,2,2) = p%KBF (K,2,2) + ElmntStff*Shape2*Shape2 + + ElmntStff = p%StiffBE(K,J)*p%DRNodes(J) ! Edgewise stiffness of blade element J + Shape = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldEdgSh(:,K), 2, ErrStat, ErrMsg ) + p%KBE (K,1,1) = p%KBE (K,1,1) + ElmntStff*Shape *Shape + + + ! Integrate to find the centrifugal-term of the generalized flapwise and edgewise + ! stiffness of the blades. Ignore the cross-correlation terms of KBFCent (i.e. + ! KBFCent(i,j) where i /= j) since these terms will never be used. + + ElmntStff = FMomAbvNd(K,J)*p%DRNodes(J)*p%RotSpeed**2 ! Centrifugal stiffness of blade element J + + Shape1 = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl1Sh(:,K), 1, ErrStat, ErrMsg ) + Shape2 = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl2Sh(:,K), 1, ErrStat, ErrMsg ) + KBFCent(K,1,1) = KBFCent(K,1,1) + ElmntStff*Shape1*Shape1 + KBFCent(K,2,2) = KBFCent(K,2,2) + ElmntStff*Shape2*Shape2 + + Shape = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldEdgSh(:,K), 1, ErrStat, ErrMsg ) + KBECent(K,1,1) = KBECent(K,1,1) + ElmntStff*Shape *Shape + + + ! Calculate the 2nd derivatives of the twisted shape functions: + + Shape = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl1Sh(:,K), 2, ErrStat, ErrMsg ) + p%TwistedSF(K,1,1,J,2) = Shape*p%CThetaS(K,J) ! 2nd deriv. of Phi1(J) for blade K + p%TwistedSF(K,2,1,J,2) = -Shape*p%SThetaS(K,J) ! 2nd deriv. of Psi1(J) for blade K + + Shape = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldFl2Sh(:,K), 2, ErrStat, ErrMsg ) + p%TwistedSF(K,1,2,J,2) = Shape*p%CThetaS(K,J) ! 2nd deriv. of Phi2(J) for blade K + p%TwistedSF(K,2,2,J,2) = -Shape*p%SThetaS(K,J) ! 2nd deriv. of Psi2(J) for blade K + + Shape = SHP( p%RNodesNorm(J), p%BldFlexL, p%BldEdgSh(:,K), 2, ErrStat, ErrMsg ) + p%TwistedSF(K,1,3,J,2) = Shape*p%SThetaS(K,J) ! 2nd deriv. of Phi3(J) for blade K + p%TwistedSF(K,2,3,J,2) = Shape*p%CThetaS(K,J) ! 2nd deriv. of Psi3(J) for blade K + + + ! Integrate to find the 1st derivatives of the twisted shape functions: + + DO I = 1,2 ! Loop through Phi and Psi + DO L = 1,3 ! Loop through all blade DOFs + TwstdSF ( I,L, 1) = p%TwistedSF(K,I,L,J,2)*0.5*p%DRNodes(J) + p%TwistedSF (K,I,L,J,1) = TwstdSF ( I,L, 1) + ENDDO ! L - All blade DOFs + ENDDO ! I - Phi and Psi + + IF ( J /= 1 ) THEN ! All but the innermost blade element + ! Add the effects from the (not yet used) portion of element J-1 + + DO I = 1,2 ! Loop through Phi and Psi + DO L = 1,3 ! Loop through all blade DOFs + p%TwistedSF(K,I,L,J,1) = p%TwistedSF(K,I,L,J,1) + p%TwistedSF(K,I,L,J-1,1) & + + TwstdSFOld( I,L, 1) + ENDDO ! L - All blade DOFs + ENDDO ! I - Phi and Psi + ENDIF + + + ! Integrate to find the twisted shape functions themselves (i.e., their zeroeth derivative): + + DO I = 1,2 ! Loop through Phi and Psi + DO L = 1,3 ! Loop through all blade DOFs + TwstdSF ( I,L, 0) = p%TwistedSF(K,I,L,J,1)*0.5*p%DRNodes(J) + p%TwistedSF (K,I,L,J,0) = TwstdSF ( I,L, 0) + ENDDO ! L - All blade DOFs + ENDDO ! I - Phi and Psi + + IF ( J /= 1 ) THEN ! All but the innermost blade element + ! Add the effects from the (not yet used) portion of element J-1 + + DO I = 1,2 ! Loop through Phi and Psi + DO L = 1,3 ! Loop through all blade DOFs + p%TwistedSF(K,I,L,J,0) = p%TwistedSF(K,I,L,J,0) + p%TwistedSF(K,I,L,J-1,0) & + + TwstdSFOld( I,L, 0) + ENDDO ! L - All blade DOFs + ENDDO ! I - Phi and Psi + ENDIF + + + ! Integrate to find the blade axial reduction shape functions: + + DO I = 1,3 ! Loop through all blade DOFs + DO L = 1,3 ! Loop through all blade DOFs + AxRdBld ( I,L ) = 0.5*p%DRNodes(J)*( & + p%TwistedSF(K,1,I,J,1)*p%TwistedSF(K,1,L,J,1) & + + p%TwistedSF(K,2,I,J,1)*p%TwistedSF(K,2,L,J,1) ) + p%AxRedBld (K,I,L,J) = AxRdBld(I,L) + ENDDO ! L - All blade DOFs + ENDDO ! I - All blade DOFs + + IF ( J /= 1 ) THEN ! All but the innermost blade element + ! Add the effects from the (not yet used) portion of element J-1 + + DO I = 1,3 ! Loop through all blade DOFs + DO L = 1,3 ! Loop through all blade DOFs + p%AxRedBld(K,I,L,J) = p%AxRedBld(K,I,L,J) + p%AxRedBld(K,I,L,J-1) & + + AxRdBldOld(I,L) + ENDDO ! L - All blade DOFs + ENDDO ! I - All blade DOFs + ENDIF + + + ! Store the TwstdSF and AxRdBld terms of the current element (these will be used for the next element) + + TwstdSFOld = TwstdSF + AxRdBldOld = AxRdBld + + + ENDDO ! J - Blade nodes / elements + + + + + IF (p%BD4Blades) THEN + + !p%KBF ( K,:,: ) = 0.0_ReKi + + ! the 1st and zeroeth derivatives of the twisted shape functions at the blade root: + p%TwistedSF(K,:,:,:,1) = 0.0_ReKi + p%TwistedSF(K,:,:,:,0) = 0.0_ReKi + p%AxRedBld( K,:,:,: ) = 0.0_ReKi + ELSE + + ! Apply the flapwise modal stiffness tuners of the blades to KBF(): + + DO I = 1,2 ! Loop through flap DOFs + DO L = 1,2 ! Loop through flap DOFs + p%KBF(K,I,L) = SQRT( p%FStTunr(K,I)*p%FStTunr(K,L) )*p%KBF(K,I,L) + ENDDO ! L - Flap DOFs + ENDDO ! I - Flap DOFs + + ! Calculate the blade natural frequencies: + + DO I = 1,NumBF ! Loop through flap DOFs + p%FreqBF(K,I,1) = Inv2Pi*SQRT( p%KBF(K,I,I) /( MBF(K,I,I) - p%TipMass(K) ) ) ! Natural blade I-flap frequency w/o centrifugal stiffening nor tip mass effects + p%FreqBF(K,I,2) = Inv2Pi*SQRT( p%KBF(K,I,I) / MBF(K,I,I) ) ! Natural blade I-flap frequency w/o centrifugal stiffening, but w/ tip mass effects + p%FreqBF(K,I,3) = Inv2Pi*SQRT( ( p%KBF(K,I,I) + KBFCent(K,I,I) )/ MBF(K,I,I) ) ! Natural blade I-flap frequency w/ centrifugal stiffening and tip mass effects + ENDDO ! I - Flap DOFs + + p%FreqBE (K,1,1) = Inv2Pi*SQRT( p%KBE(K,1,1) /( MBE(K,1,1) - p%TipMass(K) ) ) ! Natural blade 1-edge frequency w/o centrifugal stiffening nor tip mass effects + p%FreqBE (K,1,2) = Inv2Pi*SQRT( p%KBE(K,1,1) / MBE(K,1,1) ) ! Natural Blade 1-edge frequency w/o centrifugal stiffening, but w/ tip mass effects + p%FreqBE (K,1,3) = Inv2Pi*SQRT( ( p%KBE(K,1,1) + KBECent(K,1,1) )/ MBE(K,1,1) ) ! Natural Blade 1-edge frequency w/ centrifugal stiffening and tip mass effects + + + ! Calculate the generalized damping of the blades: + + DO I = 1,NumBF ! Loop through flap DOFs + DO L = 1,NumBF ! Loop through flap DOFs + p%CBF(K,I,L) = ( 0.01*p%BldFDamp(K,L) )*p%KBF(K,I,L)/( Pi*p%FreqBF(K,L,1) ) + ENDDO ! L - Flap DOFs + ENDDO ! I - Flap DOFs + + p%CBE (K,1,1) = ( 0.01*p%BldEDamp(K,1) )*p%KBE(K,1,1)/( Pi*p%FreqBE(K,1,1) ) + + + ! Calculate the 2nd derivatives of the twisted shape functions at the blade root: + + Shape = SHP( 0.0_ReKi, p%BldFlexL, p%BldFl1Sh(:,K), 2, ErrStat, ErrMsg ) + p%TwistedSF(K,1,1,0,2) = Shape*p%CThetaS(K,0) ! 2nd deriv. of Phi1(0) for blade K + p%TwistedSF(K,2,1,0,2) = -Shape*p%SThetaS(K,0) ! 2nd deriv. of Psi1(0) for blade K + + Shape = SHP( 0.0_ReKi, p%BldFlexL, p%BldFl2Sh(:,K), 2, ErrStat, ErrMsg ) + p%TwistedSF(K,1,2,0,2) = Shape*p%CThetaS(K,0) ! 2nd deriv. of Phi2(0) for blade K + p%TwistedSF(K,2,2,0,2) = -Shape*p%SThetaS(K,0) ! 2nd deriv. of Psi2(0) for blade K + + Shape = SHP( 0.0_ReKi, p%BldFlexL, p%BldEdgSh(:,K), 2, ErrStat, ErrMsg ) + p%TwistedSF(K,1,3,0,2) = Shape*p%SThetaS(K,0) ! 2nd deriv. of Phi3(0) for blade K + p%TwistedSF(K,2,3,0,2) = Shape*p%CThetaS(K,0) ! 2nd deriv. of Psi3(0) for blade K + + + ! Calculate the 2nd derivatives of the twisted shape functions at the tip: + + Shape = SHP( 1.0_ReKi, p%BldFlexL, p%BldFl1Sh(:,K), 2, ErrStat, ErrMsg ) + p%TwistedSF(K,1,1,p%TipNode,2) = Shape*p%CThetaS(K,p%TipNode) ! 2nd deriv. of Phi1(p%TipNode) for blade K + p%TwistedSF(K,2,1,p%TipNode,2) = -Shape*p%SThetaS(K,p%TipNode) ! 2nd deriv. of Psi1(p%TipNode) for blade K + + Shape = SHP( 1.0_ReKi, p%BldFlexL, p%BldFl2Sh(:,K), 2, ErrStat, ErrMsg ) + p%TwistedSF(K,1,2,p%TipNode,2) = Shape*p%CThetaS(K,p%TipNode) ! 2nd deriv. of Phi2(p%TipNode) for blade K + p%TwistedSF(K,2,2,p%TipNode,2) = -Shape*p%SThetaS(K,p%TipNode) ! 2nd deriv. of Psi2(p%TipNode) for blade K + + Shape = SHP( 1.0_ReKi, p%BldFlexL, p%BldEdgSh(:,K), 2, ErrStat, ErrMsg ) + p%TwistedSF(K,1,3,p%TipNode,2) = Shape*p%SThetaS(K,p%TipNode) ! 2nd deriv. of Phi3(p%TipNode) for blade K + p%TwistedSF(K,2,3,p%TipNode,2) = Shape*p%CThetaS(K,p%TipNode) ! 2nd deriv. of Psi3(p%TipNode) for blade K + + + ! Integrate to find the 1st and zeroeth derivatives of the twisted shape functions + ! at the tip: + + DO I = 1,2 ! Loop through Phi and Psi + DO L = 1,3 ! Loop through all blade DOFs + p%TwistedSF(K,I,L,p%TipNode,1) = p%TwistedSF(K,I,L,p%BldNodes,1) + TwstdSFOld(I,L,1) + p%TwistedSF(K,I,L,p%TipNode,0) = p%TwistedSF(K,I,L,p%BldNodes,0) + TwstdSFOld(I,L,0) + ENDDO ! L - All blade DOFs + ENDDO ! I - Phi and Psi + + ! the 1st and zeroeth derivatives of the twisted shape functions at the blade root: + p%TwistedSF(K,:,:,0,1) = 0.0_ReKi + p%TwistedSF(K,:,:,0,0) = 0.0_ReKi + p%AxRedBld( K,:,:,0 ) = 0.0_ReKi + + ! Integrate to find the blade axial reduction shape functions at the tip: + + DO I = 1,3 ! Loop through all blade DOFs + DO L = 1,3 ! Loop through all blade DOFs + p%AxRedBld(K,I,L,p%TipNode) = p%AxRedBld(K,I,L,p%BldNodes) + AxRdBldOld(I,L) + ENDDO ! L - All blade DOFs + ENDDO ! I - All blade DOFs + END IF ! p%BD4Blades + + + ENDDO ! K - Blades + + + + ! Calculate the tower-top mass: + + p%TwrTpMass = p%RotMass + p%RFrlMass + p%BoomMass + p%TFinMass + p%NacMass + p%YawBrMass + + + DO J = p%TwrNodes,1,-1 ! Loop through the tower nodes / elements in reverse + + + ! Calculate the mass of the current element + + p%TElmntMass(J) = p%MassT(J)*abs(p%DHNodes(J)) ! Mass of tower element J + + + ! Integrate to find the tower mass which will be output in .fsm + + p%TwrMass = p%TwrMass + p%TElmntMass(J) + + + ! Integrate to find TMssAbvNd: + + TMssAbvNd (J) = 0.5*p%TElmntMass(J) + + IF ( J == p%TwrNodes ) THEN ! Uppermost tower element + ! Add the TwrTpMass effects: + + TMssAbvNd(J) = TMssAbvNd(J) + p%TwrTpMass + ELSE ! All other tower elements + ! Add to TMssAbvNd(J) the effects from the (not yet used) portion of element J+1 + + TMssAbvNd(J) = 0.5*p%TElmntMass(J+1) + TMssAbvNd(J) + TMssAbvNd(J+1) + ENDIF + + + ENDDO ! J - Tower nodes / elements in reverse + + + + ! Initialize the generalized tower masses using tower-top mass effects: + + DO I = 1,2 ! Loop through all tower modes in a single direction + MTFA(I,I) = p%TwrTpMass + MTSS(I,I) = p%TwrTpMass + ENDDO ! I - All tower modes in a single direction + + ! set values for tower base (note that we haven't corrctly defined the values for (:,0,2) in the arrays below): + p%TwrFASF( :,0,0:1) = 0.0_ReKi + p%TwrSSSF( :,0,0:1) = 0.0_ReKi + p%AxRedTFA(:,:,0) = 0.0_ReKi + p%AxRedTSS(:,:,0) = 0.0_ReKi + + DO J = 1,p%TwrNodes ! Loop through the tower nodes / elements + + + ! Calculate the tower shape functions (all derivatives): + + p%TwrFASF(1,J,2) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwFAM1Sh(:), 2, ErrStat, ErrMsg ) + p%TwrFASF(2,J,2) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwFAM2Sh(:), 2, ErrStat, ErrMsg ) + p%TwrFASF(1,J,1) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwFAM1Sh(:), 1, ErrStat, ErrMsg ) + p%TwrFASF(2,J,1) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwFAM2Sh(:), 1, ErrStat, ErrMsg ) + p%TwrFASF(1,J,0) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwFAM1Sh(:), 0, ErrStat, ErrMsg ) + p%TwrFASF(2,J,0) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwFAM2Sh(:), 0, ErrStat, ErrMsg ) + + p%TwrSSSF(1,J,2) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwSSM1Sh(:), 2, ErrStat, ErrMsg ) + p%TwrSSSF(2,J,2) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwSSM2Sh(:), 2, ErrStat, ErrMsg ) + p%TwrSSSF(1,J,1) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwSSM1Sh(:), 1, ErrStat, ErrMsg ) + p%TwrSSSF(2,J,1) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwSSM2Sh(:), 1, ErrStat, ErrMsg ) + p%TwrSSSF(1,J,0) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwSSM1Sh(:), 0, ErrStat, ErrMsg ) + p%TwrSSSF(2,J,0) = SHP( p%HNodesNorm(J), p%TwrFlexL, InputFileData%TwSSM2Sh(:), 0, ErrStat, ErrMsg ) + + + ! Integrate to find the generalized mass of the tower (including tower-top mass effects). + ! Ignore the cross-correlation terms of MTFA (i.e. MTFA(i,j) where i /= j) and MTSS + ! since these terms will never be used. + + + DO I = 1,2 ! Loop through all tower DOFs in one direction + MTFA (I,I) = MTFA (I,I) + p%TElmntMass(J)*p%TwrFASF(I,J,0)**2 + MTSS (I,I) = MTSS (I,I) + p%TElmntMass(J)*p%TwrSSSF(I,J,0)**2 + ENDDO ! I - through all tower DOFs in one direction + + + ! Integrate to find the generalized stiffness of the tower (not including gravitational + ! effects). + + ElStffFA = p%StiffTFA(J)*abs(p%DHNodes(J)) ! Fore-aft stiffness of tower element J + ElStffSS = p%StiffTSS(J)*abs(p%DHNodes(J)) ! Side-to-side stiffness of tower element J + + DO I = 1,2 ! Loop through all tower DOFs in one direction + DO L = 1,2 ! Loop through all tower DOFs in one direction + p%KTFA (I,L) = p%KTFA (I,L) + ElStffFA *p%TwrFASF(I,J,2)*p%TwrFASF(L,J,2) + p%KTSS (I,L) = p%KTSS (I,L) + ElStffSS *p%TwrSSSF(I,J,2)*p%TwrSSSF(L,J,2) + ENDDO ! L - All tower DOFs in one direction + ENDDO ! I - through all tower DOFs in one direction + + + ! Integrate to find the gravitational-term of the generalized stiffness of the tower. + ! Ignore the cross-correlation terms of KTFAGrav (i.e. KTFAGrav(i,j) where i /= j) + ! and KTSSGrav since these terms will never be used. + + ElmntStff = -TMssAbvNd(J)*abs(p%DHNodes(J))*p%Gravity ! Gravitational stiffness of tower element J + + DO I = 1,2 ! Loop through all tower DOFs in one direction + KTFAGrav(I,I) = KTFAGrav(I,I) + ElmntStff*p%TwrFASF(I,J,1)**2 + KTSSGrav(I,I) = KTSSGrav(I,I) + ElmntStff*p%TwrSSSF(I,J,1)**2 + ENDDO + + + ! Integrate to find the tower axial reduction shape functions: + + DO I = 1,2 ! Loop through all tower DOFs in one direction + DO L = 1,2 ! Loop through all tower DOFs in one direction + AxRdTFA (I,L) = 0.5*p%DHNodes(J)*p%TwrFASF(I,J,1)*p%TwrFASF(L,J,1) + AxRdTSS (I,L) = 0.5*p%DHNodes(J)*p%TwrSSSF(I,J,1)*p%TwrSSSF(L,J,1) + + p%AxRedTFA(I,L,J) = AxRdTFA(I,L) + p%AxRedTSS(I,L,J) = AxRdTSS(I,L) + ENDDO ! L - All tower DOFs in one direction + ENDDO + + IF ( J /= 1 ) THEN ! All but the lowermost tower element + ! Add the effects from the (not yet used) portion of element J-1 + + DO I = 1,2 ! Loop through all tower DOFs in one direction + DO L = 1,2 ! Loop through all tower DOFs in one direction + p%AxRedTFA(I,L,J) = p%AxRedTFA(I,L,J) + p%AxRedTFA(I,L,J-1)+ AxRdTFAOld(I,L) + p%AxRedTSS(I,L,J) = p%AxRedTSS(I,L,J) + p%AxRedTSS(I,L,J-1)+ AxRdTSSOld(I,L) + ENDDO ! L - All tower DOFs in one direction + ENDDO + ENDIF + + + ! Store the AxRdTFA and AxRdTSS terms of the current element (these will be used for the next element) + + AxRdTFAOld = AxRdTFA + AxRdTSSOld = AxRdTSS + + + ENDDO ! J - Tower nodes / elements + + + ! Apply the modal stiffness tuners of the tower to KTFA() and KTSS(): + + DO I = 1,2 ! Loop through all tower DOFs in one direction + DO L = 1,2 ! Loop through all tower DOFs in one direction + p%KTFA(I,L) = SQRT( InputFileData%FAStTunr(I)*InputFileData%FAStTunr(L) )*p%KTFA(I,L) + + p%KTSS(I,L) = SQRT( InputFileData%SSStTunr(I)*InputFileData%SSStTunr(L) )*p%KTSS(I,L) + ENDDO ! L - All tower DOFs in one direction + ENDDO ! I - through all tower DOFs in one direction + + + ! Calculate the tower natural frequencies: + + DO I = 1,2 ! Loop through all tower DOFs in one direction + if ( EqualRealNos(( MTFA(I,I) - p%TwrTpMass ), 0.0_ReKi) ) then + p%FreqTFA(I,1) = NaN ! Avoid creating a divide by zero signal, but set p%FreqTFA(I,1) = NaN + else + p%FreqTFA(I,1) = Inv2Pi*SQRT( p%KTFA(I,I)/( MTFA(I,I) - p%TwrTpMass ) ) ! Natural tower I-fore-aft frequency w/o gravitational destiffening nor tower-top mass effects + end if + if ( EqualRealNos(( MTSS(I,I) - p%TwrTpMass ), 0.0_ReKi) ) then + p%FreqTSS(I,1) = NaN ! Avoid creating a divide by zero signal, but set p%FreqTFS(I,1) = NaN + else + p%FreqTSS(I,1) = Inv2Pi*SQRT( p%KTSS(I,I)/( MTSS(I,I) - p%TwrTpMass ) ) ! Natural tower I-side-to-side frequency w/o gravitational destiffening nor tower-top mass effects + end if + p%FreqTFA(I,2) = Inv2Pi*SQRT( ( p%KTFA(I,I) + KTFAGrav(I,I) )/MTFA(I,I) ) ! Natural tower I-fore-aft frequency w/ gravitational destiffening and tower-top mass effects + p%FreqTSS(I,2) = Inv2Pi*SQRT( ( p%KTSS(I,I) + KTSSGrav(I,I) )/MTSS(I,I) ) ! Natural tower I-side-to-side frequency w/ gravitational destiffening and tower-top mass effects + ENDDO ! I - All tower DOFs in one direction + + + ! Calculate the generalized damping of the tower: + + DO I = 1,2 ! Loop through all tower DOFs in one direction + DO L = 1,2 ! Loop through all tower DOFs in one direction + p%CTFA(I,L) = ( 0.01*InputFileData%TwrFADmp(L) )*p%KTFA(I,L)/( Pi*p%FreqTFA(L,1) ) + + p%CTSS(I,L) = ( 0.01*InputFileData%TwrSSDmp(L) )*p%KTSS(I,L)/( Pi*p%FreqTSS(L,1) ) + ENDDO ! L - All tower DOFs in one direction + ENDDO ! I - All tower DOFs in one direction + + + ! Calculate the tower shape functions (all derivatives) at the tower-top: + + p%TwrFASF(1,p%TTopNode,2) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwFAM1Sh(:), 2, ErrStat, ErrMsg ) + p%TwrFASF(2,p%TTopNode,2) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwFAM2Sh(:), 2, ErrStat, ErrMsg ) + p%TwrFASF(1,p%TTopNode,1) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwFAM1Sh(:), 1, ErrStat, ErrMsg ) + p%TwrFASF(2,p%TTopNode,1) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwFAM2Sh(:), 1, ErrStat, ErrMsg ) + p%TwrFASF(1,p%TTopNode,0) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwFAM1Sh(:), 0, ErrStat, ErrMsg ) + p%TwrFASF(2,p%TTopNode,0) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwFAM2Sh(:), 0, ErrStat, ErrMsg ) + + p%TwrSSSF(1,p%TTopNode,2) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwSSM1Sh(:), 2, ErrStat, ErrMsg ) + p%TwrSSSF(2,p%TTopNode,2) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwSSM2Sh(:), 2, ErrStat, ErrMsg ) + p%TwrSSSF(1,p%TTopNode,1) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwSSM1Sh(:), 1, ErrStat, ErrMsg ) + p%TwrSSSF(2,p%TTopNode,1) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwSSM2Sh(:), 1, ErrStat, ErrMsg ) + p%TwrSSSF(1,p%TTopNode,0) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwSSM1Sh(:), 0, ErrStat, ErrMsg ) + p%TwrSSSF(2,p%TTopNode,0) = SHP( 1.0_ReKi, p%TwrFlexL, InputFileData%TwSSM2Sh(:), 0, ErrStat, ErrMsg ) + + + ! Integrate to find the tower axial reduction shape functions at the tower-top: + + DO I = 1,2 ! Loop through all tower DOFs in one direction + DO L = 1,2 ! Loop through all tower DOFs in one direction + p%AxRedTFA(I,L,p%TTopNode) = p%AxRedTFA(I,L,p%TwrNodes)+ AxRdTFAOld(I,L) + p%AxRedTSS(I,L,p%TTopNode) = p%AxRedTSS(I,L,p%TwrNodes)+ AxRdTSSOld(I,L) + ENDDO ! L - All tower DOFs in one direction + ENDDO + + + ! Calculate the turbine mass: + + p%TurbMass = p%TwrTpMass + p%TwrMass + + + RETURN +END SUBROUTINE Coeff +``` \ No newline at end of file diff --git a/多体求解器debug/明阳机型转速发散 叶片变形量为负.md b/多体求解器debug/明阳机型转速发散 叶片变形量为负.md index 18dc1be..f153055 100644 --- a/多体求解器debug/明阳机型转速发散 叶片变形量为负.md +++ b/多体求解器debug/明阳机型转速发散 叶片变形量为负.md @@ -4,4 +4,16 @@ 关闭叶片柔性，转速不发散 -思路： 判断是否由于叶片变形量为负导致转速发散，尝试从叶片方程探究为什么变形量为负 \ No newline at end of file +思路： 判断是否由于叶片变形量为负导致转速发散，尝试从叶片方程探究为什么变形量为负 + + +p.kbf +p.cbf +数太大了 +明阳机型是nrel 5mw的100倍 + +p.ktfa +p.ctfa +p.ktss +p.ctss +明阳机型是nrel 5mw的1000倍 diff --git a/未命名 1.base b/未命名 1.base new file mode 100644 index 0000000..f165474 --- /dev/null +++ b/未命名 1.base @@ -0,0 +1,3 @@ +views: + - type: table + name: 表格