obsidian_backup/20260120_093252__[[meshmap]]分析SUBROUTINE_MeshMapCreate的实现方法and预弯.md at 525fdf5fcba81db21d1fb2ce43688b98bd05e5a5

gyz/obsidian_backup

Fork 0

yize ca3dc9f85a vault backup: 2026-01-23 15:42:30

2026-01-23 15:42:35 +08:00

49 KiB

Raw Blame History

epoch, modelKey, tags

epoch

modelKey

1. 运动和标量映射 (Motion and Scalar Mapping)

当 Src 和 Dest 都具有运动场 (HasMotionFields(Src) .AND. HasMotionFields(Dest)) 并且它们都是 Line2 元素网格时，MeshMapCreate 会执行以下步骤：

确定映射大小:
- PointsInMap 将被设置为 Dest%Nnodes（目标网格的节点数）。这意味着 MeshMap%MapMotions 数组将为目标网格中的每个节点分配一个条目。
- ElementNodes 将设置为 2，因为 Line2 元素有两个节点。
分配 MeshMap%MapMotions:
- 子程序会尝试分配 MeshMap%MapMotions(PointsInMap)。如果分配失败，将报告致命错误。
调用 CreateMotionMap_L2_to_L2:
- 由于 Dest%ElemTable(ELEMENT_LINE2)%nelem > 0 (目标网格有 Line2 元素) 且 Src%ElemTable(ELEMENT_LINE2)%nelem > 0 (源网格也有 Line2 元素)，MeshMapCreate 会调用 CreateMotionMap_L2_to_L2 子程序。
CreateMotionMap_L2_to_L2 的实现:
- 这个子程序的主要任务是调用 CreateMapping_ProjectToLine2。
- 调用方式为：CALL CreateMapping_ProjectToLine2(Dest, Src, MeshMap%MapMotions, ELEMENT_LINE2, ErrStat2, ErrMsg2)。
  - 这里，Mesh1 参数是 Dest 网格，Mesh2 参数是 Src 网格。
  - NodeMap 是 MeshMap%MapMotions。
  - Mesh1_TYPE 是 ELEMENT_LINE2。
  - Mesh2_TYPE 在 CreateMapping_ProjectToLine2 中被硬编码为 ELEMENT_LINE2。
CreateMapping_ProjectToLine2 的实现 (L2-to-L2 运动映射):
- 这个子程序的目标是将 Mesh1 (即 Dest 网格) 的每个节点投影到 Mesh2 (即 Src 网格) 的 Line2 元素上，以找到最佳的映射关系。
- 初始化: NodeMap 中的每个条目 NodeMap(i)%OtherMesh_Element 被初始化为 NODE_NOT_MAPPED，表示尚未映射。
- 遍历目标网格节点:
  - 它会遍历 Mesh1 (目标网格 Dest) 中所有 ELEMENT_LINE2 类型的元素 (iElem)。
  - 对于每个元素，它会遍历该元素的所有节点 (iNode)。
  - 获取当前目标节点 i 的三维坐标 Mesh1_xyz = Mesh1%Position(:, i)。
- 寻找最近的源 Line2 元素:
  - 对于每个目标节点 i，它会遍历 Mesh2 (源网格 Src) 中所有 ELEMENT_LINE2 类型的元素 (jElem)。
  - 对于每个源 Line2 元素 jElem，它会获取其两个节点 n1 和 n2 的坐标。
  - 投影计算:
    - 计算从 n1 到 n2 的向量 n1_n2_vector。
    - 计算从 n1 到目标节点 i 的向量 n1_Point_vector。
    - 通过点积计算目标节点 i 在 n1_n2_vector 方向上的投影位置 elem_position。
    - elem_position = DOT_PRODUCT(n1_n2_vector, n1_Point_vector) / DOT_PRODUCT(n1_n2_vector, n1_n2_vector)。
  - 判断是否在元素上:
    - 如果 elem_position 在 [0.0, 1.0] 之间（包括边界上的数值容差），则认为目标节点投影到了源 Line2 元素上。
    - 如果不在元素上，但非常接近元素端点（例如，小于 7.5 毫米），则会记录为“最近的元素”，并可能在后续作为映射目标。
  - 计算距离并选择最佳映射:
    - 如果目标节点投影在源 Line2 元素上，则计算目标节点到该 Line2 元素的垂直距离（实际上是距离的平方）。
    - 如果当前 Line2 元素的距离小于之前找到的最小距离 min_dist，则更新 min_dist，并将当前源 Line2 元素 jElem 及其投影位置 elem_position 记录为最佳映射。
    - NodeMap(i)%OtherMesh_Element 存储源 Line2 元素的索引 jElem。
    - NodeMap(i)%shape_fn(1) 和 NodeMap(i)%shape_fn(2) 存储插值形函数值，分别为 1.0 - elem_position 和 elem_position。这些形函数用于在源 Line2 元素上进行线性插值，以获取目标节点处的运动值。
- 错误处理: 如果一个目标节点无法投影到任何源 Line2 元素上（即使考虑了数值容差），则会报告致命错误。

总结 L2-to-L2 运动映射: 对于目标网格 Dest 中的每个节点，系统会找到源网格 Src 中离它最近的 Line2 元素，并将目标节点映射到该源 Line2 元素上的一个投影点。映射信息包括源元素的索引和在该元素上的插值权重。

2. 载荷映射 (Load Mapping)

当 Src 和 Dest 都具有载荷场 (HasLoadFields(Src) .AND. HasLoadFields(Dest)) 并且它们都是 Line2 元素网格时，MeshMapCreate 会执行以下步骤：

兼容性检查:
- 子程序会检查载荷场的兼容性。例如，如果源网格有力 (Src%FieldMask(MASKID_Force))，则目标网格必须有力 (Dest%FieldMask(MASKID_Force)) 和力矩 (Dest%FieldMask(MASKID_Moment))。如果源网格有力矩 (Src%FieldMask(MASKID_Moment))，则目标网格必须有力矩。这些是确保载荷可以正确传输的必要条件。
确定映射大小:
- PointsInMap 将被设置为 Src%Nnodes（源网格的节点数）。这意味着 MeshMap%MapLoads 数组将为源网格中的每个节点分配一个条目。
- ElementNodes 将设置为 2，因为 Line2 元素有两个节点。
分配 MeshMap%MapLoads:
- 子程序会尝试分配 MeshMap%MapLoads(PointsInMap)。如果分配失败，将报告致命错误。
调用 CreateLoadMap_L2_to_L2:
- 由于 Dest%ElemTable(ELEMENT_LINE2)%nelem > 0 (目标网格有 Line2 元素) 且 Src%ElemTable(ELEMENT_LINE2)%nelem > 0 (源网格也有 Line2 元素)，MeshMapCreate 会调用 CreateLoadMap_L2_to_L2 子程序。
CreateLoadMap_L2_to_L2 的实现 (未在提供的代码中，但可推断):
- 根据 CreateMotionMap_L2_to_L2 的模式，CreateLoadMap_L2_to_L2 也很可能会调用 CreateMapping_ProjectToLine2。
- 然而，对于载荷映射，源和目标网格的角色通常是相反的。载荷是从源网格传递到目标网格。因此，在载荷映射中，我们通常需要将源网格的节点映射到目标网格的元素上。
- 因此，推断的调用方式可能是：CALL CreateMapping_ProjectToLine2(Src, Dest, MeshMap%MapLoads, ELEMENT_LINE2, ErrStat2, ErrMsg2)。
  - 这里，Mesh1 参数是 Src 网格，Mesh2 参数是 Dest 网格。
  - NodeMap 是 MeshMap%MapLoads。
  - Mesh1_TYPE 是 ELEMENT_LINE2。
  - Mesh2_TYPE 在 CreateMapping_ProjectToLine2 中被硬编码为 ELEMENT_LINE2。
CreateMapping_ProjectToLine2 的实现 (L2-to-L2 载荷映射):
- 与运动映射类似，但角色互换。它将 Mesh1 (即 Src 网格) 的每个节点投影到 Mesh2 (即 Dest 网格) 的 Line2 元素上。
- 遍历源网格节点:
  - 它会遍历 Mesh1 (源网格 Src) 中所有 ELEMENT_LINE2 类型的元素 (iElem)。
  - 对于每个元素，它会遍历该元素的所有节点 (iNode)。
  - 获取当前源节点 i 的三维坐标 Mesh1_xyz = Mesh1%Position(:, i)。
- 寻找最近的目标 Line2 元素:
  - 对于每个源节点 i，它会遍历 Mesh2 (目标网格 Dest) 中所有 ELEMENT_LINE2 类型的元素 (jElem)。
  - 执行与运动映射中相同的投影计算和距离判断逻辑。
  - NodeMap(i)%OtherMesh_Element 存储目标 Line2 元素的索引 jElem。
  - NodeMap(i)%shape_fn(1) 和 NodeMap(i)%shape_fn(2) 存储插值形函数值，用于在目标 Line2 元素上进行线性插值，以确定源节点处的载荷如何分布到目标元素上。
- 错误处理: 同样，如果一个源节点无法投影到任何目标 Line2 元素上，则会报告致命错误。

总结 L2-to-L2 载荷映射: 对于源网格 Src 中的每个节点，系统会找到目标网格 Dest 中离它最近的 Line2 元素，并将源节点映射到该目标 Line2 元素上的一个投影点。映射信息包括目标元素的索引和在该元素上的插值权重。这些权重将用于将源节点上的载荷（力、力矩）分配到目标 Line2 元素的两个节点上。

3. `DisplacedPosition` 的分配

无论载荷或运动映射是否发生，如果 MeshMap%DisplacedPosition 尚未分配，它都会被分配。

PointsInTmpMap 在 Line2-to-Line2 情况下，对于运动映射是 Dest%Nnodes，对于载荷映射是 Src%Nnodes。因此，PointsInTmpMap 将取两者中的最大值。
ElementNodes 在 Line2 网格情况下为 2。
MeshMap%DisplacedPosition 将被分配为一个 3 x PointsInTmpMap x ElementNodes 的数组，用于存储位移后的位置信息。

总结

当您有两套 Line2 网格进行映射时：

运动映射 (L2-to-L2)：目标网格 Dest 中的每个节点会被投影到源网格 Src 中最近的 Line2 元素上。映射结果存储在 MeshMap%MapMotions 中，包含源元素的索引和插值形函数，用于从源网格插值运动到目标网格节点。
载荷映射 (L2-to-L2)：源网格 Src 中的每个节点会被投影到目标网格 Dest 中最近的 Line2 元素上。映射结果存储在 MeshMap%MapLoads 中，包含目标元素的索引和插值形函数，用于将源网格节点上的载荷分配到目标网格元素节点。

这两种映射都依赖于 CreateMapping_ProjectToLine2 子程序，该子程序通过几何投影和距离计算来确定节点与 Line2 元素之间的对应关系。关键的区别在于 CreateMapping_ProjectToLine2 被调用的参数，特别是 Mesh1 和 Mesh2 的角色，这决定了是“目标节点找源元素”还是“源节点找目标元素”。 [Timestamp: 2026/01/20 09:46:27]

user: meshmap CalculatePositions我现在的两套网格，aero部分是有风电叶片预弯的网格，elast部分叶片网格是直的，两套网格是否可以匹配 [Timestamp: 2026/01/20 10:11:20]

ai: 您提出的问题是关于两套网格的匹配性：一套是气动（aero）部分使用的预弯叶片网格，另一套是弹性（elast）部分使用的直叶片网格。在 [[meshmap]] 和 [[CalculatePositions]] 的上下文中，这两套网格是否可以匹配，以及如何匹配。

核心观点：

从 MeshMapCreate 的实现来看，它主要依赖于网格节点的几何位置来建立映射关系。因此，如果两套网格的节点在初始（未变形）状态下的几何位置差异很大，那么直接进行映射可能会导致不准确或失败。

然而，OpenFAST（以及类似的耦合仿真工具）通常会处理这种几何差异，通过在仿真过程中动态更新网格位置来实现匹配。

`MeshMapCreate` 的作用和限制

MeshMapCreate 子程序负责在仿真开始时建立初始的网格映射关系。它会：

检查网格是否已提交: 确保网格定义完整。
分配映射数据结构: 为运动和载荷分配 MeshMap%MapMotions 和 MeshMap%MapLoads。
调用映射算法: 根据网格类型（例如，Line2-to-Line2）调用 CreateMapping_ProjectToLine2 或 CreateMapping_NearestNeighbor。

这些映射算法的核心是：

CreateMapping_ProjectToLine2: 将一个网格的节点投影到另一个网格的 Line2 元素上，找到最近的投影点和插值权重。
CreateMapping_NearestNeighbor: 找到一个网格的节点在另一个网格中的最近邻节点。

对于您的情况：

Aero 网格 (预弯叶片)：这代表了叶片的初始几何形状，包含了预弯信息。
Elast 网格 (直叶片)：这代表了叶片的未变形参考构型，通常是直的。

如果 MeshMapCreate 直接尝试将预弯网格的节点映射到直网格的节点或元素上，由于几何形状的显著差异（预弯与直），初始映射的准确性会很差，甚至可能失败。CreateMapping_ProjectToLine2 会寻找最近的投影点，但如果叶片形状差异太大，投影可能落在元素之外，或者找到的“最近”元素在物理上并不对应。

`CalculatePositions` 的作用和动态匹配

[[CalculatePositions]] 子程序在 OpenFAST 的弹性动力学 (Elastodynamics) 模块中扮演着关键角色。它的主要功能是：

计算各个点的瞬时位置: 根据系统的连续状态 x（包括位移和旋转）和参数 p，计算风力机各个部件（包括塔架和叶片）在当前时刻的瞬时位置。
考虑变形: 特别是对于叶片，CalculatePositions 会根据弹性变形状态 (x%QT(DOF_BF(K,1)) 等) 计算叶片节点的实际位置 RtHSdat%rS。这意味着，即使弹性部分的参考网格是直的，在仿真过程中，CalculatePositions 会根据叶片的变形（包括预弯和柔性变形）计算出其实际的、变形后的位置。

因此，匹配的关键在于：

初始几何差异: MeshMapCreate 在建立初始映射时，会使用网格的初始几何位置。如果 aero 网格是预弯的，而 elast 网格是直的，那么在初始时刻，这两套网格的几何形状是不同的。
动态变形: 在仿真过程中，弹性部分的直叶片网格会根据其变形状态（包括预弯效应，如果预弯被建模为一种初始变形）计算出其实际的、变形后的位置。
映射的更新: 理想情况下，映射关系应该在每个时间步或在变形显著时进行更新，以反映网格的实际几何形状。MeshMapCreate 的注释 bjj: maybe the MeshMapCreate routine shouldn't actually allocate arrays; allocate them in the "IF (RemapFlag)" sections so that if people add nodes during the simulation, the structures get reallocated to correct size? 也暗示了这种动态更新的需求。

结论和建议

两套网格是否可以匹配？

直接匹配（基于初始几何）: 如果 aero 网格是预弯的，而 elast 网格是直的，并且 MeshMapCreate 仅基于它们的初始几何位置进行映射，那么匹配会非常困难，甚至可能失败。CreateMapping_ProjectToLine2 可能会因为几何差异过大而找不到合适的投影点，或者找到的映射关系不准确。
通过动态变形匹配: 在 OpenFAST 这样的耦合仿真中，通常会通过以下方式处理：
1. 弹性网格的参考构型: 弹性部分的直叶片网格是其未变形的参考构型。
2. 预弯的建模: 预弯通常被建模为弹性叶片在初始状态下的一个变形，或者在坐标系定义中考虑。CalculatePositions 会将这种预弯效应（以及其他柔性变形）计算到叶片的实际位置中。
3. 气动网格的变形: 气动网格（预弯叶片）在仿真过程中也会根据弹性部分的变形进行更新，使其与弹性叶片的实际变形形状保持一致。
4. 映射的动态更新: 关键在于，在每个时间步，当弹性叶片发生变形时，气动网格的实际位置会根据弹性叶片的变形进行更新。然后，映射算法（或其后续的数据传输逻辑）会使用这些更新后的、变形后的位置来传递载荷和运动。

为了实现准确的匹配，通常需要确保：

几何一致性: 在进行映射时，两套网格的实际几何形状应该尽可能接近。这意味着，如果 aero 网格是预弯的，那么在进行映射时，elast 网格也应该被“变形”到预弯状态，或者 aero 网格的节点位置是基于 elast 网格的变形来计算的。
坐标系和参考点: 确保两套网格的坐标系定义和参考点是一致的，或者有明确的转换关系。
映射算法的鲁棒性: CreateMapping_ProjectToLine2 包含了一个 7.5 毫米的容差，用于处理节点非常接近元素端点的情况。这在一定程度上增加了鲁棒性，但如果初始几何差异过大，这个容差可能不足以解决问题。

具体建议：

检查 OpenFAST 的具体实现: 了解 OpenFAST 如何处理预弯叶片的气动-弹性耦合。通常，气动网格的节点位置会根据弹性叶片的变形（包括预弯）进行更新。这意味着，在进行映射时，气动网格的节点位置实际上是变形后的位置。
确保 Src 和 Dest 在映射时具有相似的几何:
- 如果 Src 是 aero 网格（预弯），Dest 是 elast 网格（直），那么在调用 MeshMapCreate 之前，可能需要将 Dest 网格的节点位置预先变形到与 Src 相似的预弯状态，或者 Src 网格的初始位置是基于 Dest 网格的参考构型加上预弯变形来定义的。
- 更常见的情况是，在每个时间步，弹性模块会计算出叶片的实际变形位置（包括预弯），然后气动模块会使用这些位置来构建其气动网格。映射发生在这些瞬时变形位置之间。
关注错误信息: 如果 MeshMapCreate 报告致命错误，例如“Node ... does not project onto any line2 element”，这通常意味着几何差异过大，导致映射失败。

总之，虽然 MeshMapCreate 依赖于几何位置，但 OpenFAST 的整体框架会通过 CalculatePositions 等子程序在仿真过程中动态处理叶片的变形。因此，关键在于确保在进行网格映射时，所使用的网格位置（无论是源网格还是目标网格）能够准确反映叶片的实际几何形状，包括预弯效应。如果初始几何差异过大，直接的静态映射会存在问题，需要通过动态更新网格位置来解决。 [Timestamp: 2026/01/20 10:11:36]

user: CalculatePositions CalculateAngularPosVelPAcc如果我在结构上使用预弯的叶片模型和网格，如何使用模态计算叶片变形量 [Timestamp: 2026/01/21 14:35:09]

ai: 您的问题是关于在结构上使用预弯叶片模型和网格时，如何使用模态计算叶片变形量，结合 [[CalculatePositions]] 和 [[CalculateAngularPosVelPAcc]] 子程序。

在 OpenFAST 这样的多体动力学仿真工具中，叶片的变形通常通过模态叠加法来表示。这意味着叶片的总变形是其各个振动模态乘以相应的模态坐标（或广义坐标）的线性组合。

1. 预弯叶片模型和网格的结构表示

当您说“结构上使用预弯的叶片模型和网格”时，这通常意味着：

初始几何: 您的叶片网格在未变形状态下就已经包含了预弯的几何形状。这与之前讨论的“直叶片网格”作为参考构型不同。在这种情况下，叶片的“零变形”位置已经是预弯的。
模态形状: 弹性模态（如挥舞、摆振、扭转模态）是基于这个预弯几何来计算的。模态形状函数 p%TwistedSF 和 p%AxRedBld 包含了叶片沿展向的几何和结构特性，它们会反映预弯的影响。

2. `CalculatePositions` 如何计算叶片变形量

CalculatePositions 子程序的核心任务是计算风力机各个部件在当前时刻的瞬时位置。对于叶片，它通过模态叠加来计算其变形。

我们来看 CalculatePositions 中计算叶片节点位置的关键部分：

   DO K = 1,p%NumBl ! Loop through all blades
      ...
      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) ) &  ! 展向位移分量 (j1方向)
                                + 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) ) &  ! 弦向位移分量 (j2方向)
                                + 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* &                           ! 展向位置分量 (j3方向)
                              (      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,:)
         RtHSdat%rS  (:,K,J) = RtHSdat%rQS (:,K,J) + RtHSdat%rQ
      END DO
   END DO

解释：

x%QT(...): x 是连续状态变量，x%QT 包含了广义坐标（模态坐标）。
- DOF_BF(K,1) 和 DOF_BF(K,2) 通常代表叶片 K 的第一和第二阶挥舞 (Flapwise) 模态的模态坐标。
- DOF_BE(K,1) 通常代表叶片 K 的第一阶摆振 (Edgewise) 模态的模态坐标。
- （注意：这里没有直接看到扭转模态的广义坐标，但通常也会有类似的表示，可能在 DOF_BE 中或通过其他方式处理。）
p%TwistedSF(K,dim,mode,J,0): 这是模态形状函数 (Mode Shape Function)。
- K: 叶片编号。
- dim: 位移分量方向（1代表 j1 方向，2代表 j2 方向）。
- mode: 模态索引（例如，1代表第一阶挥舞，2代表第二阶挥舞，3代表第一阶摆振）。
- J: 叶片沿展向的节点索引。
- 0: 表示位移形状函数（而不是旋转形状函数）。
- 这些形状函数描述了在给定模态下，叶片沿展向不同位置的单位变形量。它们是预先通过有限元分析（FEA）等方法计算得到的，并且已经包含了叶片的初始预弯几何。
CoordSys%j1(K,:), CoordSys%j2(K,:), CoordSys%j3(K,:): 这些是叶片 K 局部坐标系的单位向量。j3 通常沿叶片展向，j1 和 j2 垂直于 j3，分别代表挥舞和摆振方向。
p%RNodes(J): 这是叶片节点 J 沿展向的未变形径向位置。
p%AxRedBld(...): 这些项是轴向缩短 (Axial Shortening) 效应的系数，表示由于大变形引起的展向长度缩短。它们是模态坐标的二次项，用于修正叶片沿展向的实际位置。

如何使用模态计算叶片变形量：

叶片上任意一点 J 的变形位移可以表示为：

\vec{u}_J = \sum_{m} q_m \vec{\phi}_{m,J}

其中：

\vec{u}_J 是节点 J 的变形位移向量。
q_m 是第 m 阶模态的模态坐标（即 x%QT(DOF_BF(K,1)) 等）。
\vec{\phi}_{m,J} 是第 m 阶模态在节点 J 处的模态形状向量。

在 CalculatePositions 中，RtHSdat%rS0S(:,K,J) 实际上计算的是从叶片根部 (S(0)) 到当前节点 (S(RNodes(J))) 的相对位置向量。这个向量包含了：

未变形的展向位置: p%RNodes(J) 沿 CoordSys%j3(K,:) 方向。
模态引起的变形位移:
- j1 方向的变形：p%TwistedSF(K,1,mode,J,0) * x%QT(DOF_BF/BE)
- j2 方向的变形：p%TwistedSF(K,2,mode,J,0) * x%QT(DOF_BF/BE)
轴向缩短效应: 0.5 * (p%AxRedBld(...) * x%QT(...) * x%QT(...)) 沿 CoordSys%j3(K,:) 方向。

因此，叶片变形量就是当前计算出的位置与叶片在初始预弯状态下的未变形位置之间的差异。由于 p%TwistedSF 已经包含了预弯几何，所以 x%QT 乘以 p%TwistedSF 得到的位移就是相对于预弯构型的柔性变形。

3. `CalculateAngularPosVelPAcc` 的作用

CalculateAngularPosVelPAcc 子程序计算各个部件的角位置、角速度和角加速度。对于叶片，它也使用模态叠加来计算其弹性旋转。

关键部分如下：

   DO K = 1,p%NumBl ! Loop through all blades
      DO J = 0,p%TipNode ! Loop through the blade nodes / elements
         ...
         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,:)
         ...
      END DO
   END DO

解释：

p%TwistedSF(K,dim,mode,J,1): 这里的 TwistedSF 带有 1 作为最后一个索引，表示这是旋转形状函数（或模态斜率/扭转形状函数），而不是位移形状函数。它描述了在给定模态下，叶片沿展向不同位置的单位旋转量。
x%QDT(...): x%QDT 包含了广义速度（模态坐标的时间导数）。
AngVelHM: 这是叶片节点 J 相对于轮毂 (H) 的弹性角速度。它通过模态速度 x%QDT 乘以旋转形状函数来计算。
AngPosHM: 这是叶片节点 J 相对于轮毂 (H) 的弹性角位置。它通过模态坐标 x%QT 乘以旋转形状函数来计算。

叶片变形量（位移和旋转）的获取：

如果您想获取叶片上某个节点 J 的总变形位移和总变形旋转，您需要：

位移:
- 从 RtHSdat%rS(:,K,J) 中减去叶片 K 节点 J 的初始预弯位置（即在 x%QT 都为零时的位置）。
- RtHSdat%rS(:,K,J) 已经是相对于惯性坐标系的总位置。
- 相对变形位移：RtHSdat%rS0S(:,K,J) 已经包含了相对于叶片根部的变形位移（包括柔性变形和轴向缩短）。如果您需要相对于预弯构型的纯柔性变形，那么 RtHSdat%rS0S(:,K,J) 的计算结果就是这个。
旋转:
- RtHSdat%AngPosHM(:,K,J) 直接提供了叶片节点 J 相对于轮毂的弹性角位置（即弹性旋转变形量）。这通常是相对于预弯构型的额外旋转。
- RtHSdat%AngVelEM(:,J,K) 提供了叶片节点 J 在惯性坐标系中的总角速度，其中包含了轮毂的刚体运动和叶片的弹性变形引起的角速度。

总结：

当结构上使用预弯的叶片模型和网格时，OpenFAST 已经将预弯几何信息编码在模态形状函数 p%TwistedSF 和 p%AxRedBld 中。

CalculatePositions 通过将模态坐标 x%QT 与位移模态形状函数 p%TwistedSF(...,0) 叠加，并考虑轴向缩短效应 p%AxRedBld，来计算叶片节点相对于其根部的变形位移。这个位移是相对于预弯构型的柔性变形。
CalculateAngularPosVelPAcc 通过将模态坐标 x%QT 与旋转模态形状函数 p%TwistedSF(...,1) 叠加，来计算叶片节点相对于轮毂的弹性角位置（旋转变形量）。

因此，您不需要额外处理预弯，因为模态形状函数已经包含了预弯的影响。x%QT 中的模态坐标直接驱动了相对于预弯构型的柔性变形。 [Timestamp: 2026/01/21 14:35:31]

user: 我当前模态叠加的技术路线，支持叶片预弯变形，在模态形状的使用方法上是否有区别 [Timestamp: 2026/01/21 14:42:16]

ai: 如果您当前的模态叠加技术路线支持叶片预弯变形，那么在模态形状的使用方法上，核心原理没有区别，但模态形状的生成方式和物理含义会包含预弯信息。

让我们详细分析一下：

1. 模态叠加法的核心原理

模态叠加法的基本思想是将结构的任意变形表示为一系列预定义的模态形状（或振型）的线性组合。每个模态形状乘以一个随时间变化的模态坐标（或广义坐标），然后将这些乘积相加，得到结构的瞬时变形。

\vec{u}(s, t) = \sum_{i=1}^{N} q_i(t) \vec{\phi}_i(s)

其中：

\vec{u}(s, t) 是叶片沿展向位置 s 在时间 t 的变形位移。
q_i(t) 是第 i 阶模态的模态坐标（广义坐标），由仿真求解得到。
\vec{\phi}_i(s) 是第 i 阶模态在展向位置 s 处的模态形状向量。

2. 模态形状中包含预弯信息

当您的结构模型（例如，用于有限元分析生成模态的几何模型）本身就包含了叶片的预弯几何时，生成的模态形状 \vec{\phi}_i(s) 会自然地反映这种预弯。

这意味着：

参考构型: 此时，叶片的“未变形”或“零模态坐标”状态，就是其预弯后的几何形状。
模态形状的物理含义: 模态形状 \vec{\phi}_i(s) 描述的是叶片在预弯构型基础上的柔性变形。例如，第一阶挥舞模态描述的是叶片在预弯状态下，如何进一步向上或向下弯曲。
p%TwistedSF 的作用: 在 OpenFAST 的 CalculatePositions 子程序中，p%TwistedSF 数组就是这些模态形状函数。如果您的输入模型是预弯的，那么这些 p%TwistedSF 值就是从预弯叶片模型中提取或计算得到的。它们已经包含了预弯对模态形状的影响。

3. 模态形状使用方法的区别（概念上而非代码实现上）

从代码实现的角度来看，CalculatePositions 和 CalculateAngularPosVelPAcc 子程序中模态形状 (p%TwistedSF) 的使用方式没有区别。它们仍然是将模态坐标 (x%QT) 乘以相应的模态形状函数来计算位移和旋转。

然而，在概念和数据准备上存在关键区别：

传统直叶片模型:
- 模态形状 \vec{\phi}_i(s) 是基于直叶片的参考构型计算的。
- 预弯效应通常需要通过额外的几何变换或初始位移来引入，或者在气动网格中直接使用预弯几何。
- x%QT 乘以 \vec{\phi}_i(s) 得到的是相对于直叶片构型的变形。
预弯叶片模型:
- 模态形状 \vec{\phi}_i(s) 是基于预弯叶片的参考构型计算的。
- 预弯效应已经内嵌在模态形状中，不需要额外的几何变换来引入预弯。
- x%QT 乘以 \vec{\phi}_i(s) 得到的是相对于预弯叶片构型的柔性变形。

这意味着：

模态形状的生成: 您需要确保用于生成 p%TwistedSF 和 p%AxRedBld 等参数的有限元模型（或等效梁模型）本身就包含了叶片的预弯几何。这是最关键的一步。
“零变形”状态: 当所有模态坐标 x%QT 都为零时，叶片的位置将是其预弯后的初始几何形状，而不是完全笔直的。
变形的解释: 此时，通过模态叠加计算出的位移和旋转，是叶片在预弯构型基础上的额外柔性变形。

总结

如果您当前的模态叠加技术路线支持叶片预弯变形，并且您已经使用了包含预弯几何的结构模型来生成模态形状（p%TwistedSF 等），那么：

代码中的模态形状使用方法不变：CalculatePositions 和 CalculateAngularPosVelPAcc 中的计算公式保持不变。
模态形状的物理含义不同：这些模态形状本身已经包含了预弯的影响，它们描述的是相对于预弯构型的柔性变形。
数据准备是关键：确保输入的模态形状数据 (p%TwistedSF, p%AxRedBld) 是从预弯叶片模型中正确提取的。

这种方法是处理预弯叶片的一种常见且有效的方式，因为它将预弯视为结构初始几何的一部分，从而简化了后续的动力学计算。 [Timestamp: 2026/01/21 14:42:27]

49 KiB Raw Blame History Unescape Escape