obsidian_backup/多体+耦合求解器/CaseToolBox项目学习/casetoolbox项目分析.md

CASEToolBox 是由DTU M. H. Hansen开发的一款用于执行风电机组计算气动-伺服-弹性分析的软件包。该代码仍在开发中。该软件包包含用于进行功率和稳定性分析的工具 CASEStab，以及用于分析二维翼型的推力变形和扭角的空气动力学阻尼工具 CASEDamp。目前，CASEStab 仅可用于计算具有相同叶片的风轮在均匀气流中无重力作用下的稳态，以及计算结构的叶片模态频率和模态。

casetoolbox
├── casestab/
└── casedamp/

casestab包分析

casetoolbox/casestab/
├── __init__.py
├── aerodynamics.py
├── casestab.py
├── corotbeam_precompiled_functions.py
├── corotbeam.py
├── generic_model_components.py
├── HAWC2_blade_translator.py
├── math_functions.py
├── model_assembler.py
├── model_precompiled_functions.py
├── rigidbody.py
├── timoshenko_beam_section.py
├── wake_model.py
├── wind_model.py
├── docs/
└── tests/

• casetoolbox/casestab/__init__.py: 这是一个空的 __init__.py 文件，表明 casestab 是一个 Python 包。

• casetoolbox/casestab/aerodynamics.py: 该文件定义了风力涡轮机叶片的空气动力学模型和状态，包括 aero_blade 和 aero_calc_point 类，以及用于读取和插值翼型极线数据的函数。

• casetoolbox/casestab/casestab.py: 这是 CASEStab 项目的核心文件，负责风力涡轮机转子的稳态计算和变桨曲线调优。它通过 rotor_models 类组装整个模型，并执行稳态分析和优化。

• casetoolbox/casestab/corotbeam_precompiled_functions.py: 该文件包含使用 Numba 进行 JIT 编译的底层数学和力学函数，用于处理共旋转梁单元的复杂运动学和动力学计算，以提高性能。

• casetoolbox/casestab/corotbeam.py: 该文件实现了共旋转梁单元的有限元公式，定义了梁单元的运动学、静力学和动力学行为，并提供了从结构输入文件创建单元、更新单元状态、计算内部力和刚度、更新惯性以及处理气动弹性耦合的完整框架。

• casetoolbox/casestab/generic_model_components.py: 该文件定义了用于构建和管理复杂系统（如风力涡轮机模型）的通用模型组件，特别是与惯性、气动弹性耦合和轴承相关的组件。

• casetoolbox/casestab/HAWC2_blade_translator.py: 该文件负责读取 HAWC2 模型文件，并将其结构和气动数据转换为 CASEStab 所需的格式，实现了与 HAWC2 的互操作性。

• casetoolbox/casestab/math_functions.py: 该文件提供了一系列用于三维旋转、向量操作和数据插值的数学函数，其中一些使用 Numba 进行 JIT 编译以提高性能。

• casetoolbox/casestab/model_assembler.py: 该文件负责将风力涡轮机模型的各个部分（子结构、叶片、转子、尾流、风）组装成一个完整的系统模型，并提供进行稳态分析和结果可视化的功能。

• casetoolbox/casestab/model_precompiled_functions.py: 该文件包含使用 Numba 进行 JIT 编译的函数，主要用于计算模型中的离心力、陀螺矩阵和离心刚度矩阵，以提高性能。

• casetoolbox/casestab/rigidbody.py: 该文件定义了一个简单的刚体子结构模型，主要用于表示模型中不发生变形的部分，例如风力涡轮机的轮毂或塔基。

• casetoolbox/casestab/timoshenko_beam_section.py: 该文件主要提供了用于处理 Timoshenko 梁截面属性的函数，特别是将截面属性转换为 6x6 柔度矩阵，以及在不同参考点之间转换矩阵。

• casetoolbox/casestab/wake_model.py: 该文件主要负责风力涡轮机转子上的尾流模型，用于计算诱导速度，包括无诱导模型和轴对称诱导模型。

• casetoolbox/casestab/wind_model.py: 该文件主要负责定义风力涡轮机仿真中使用的风场模型，目前支持均匀风场。

casestab/aerodynamics.py

主要负责风力涡轮机叶片的空气动力学建模。

主要组成部分：

1. aero_blade 类：

   • 目的： 表示风力涡轮机叶片的空气动力学模型和状态。

   • 初始化 (__init__)： 从指定的几何文件 (geo_file) 和翼型极线文件 (pro_file) 读取数据。它支持 HAWC2 和 Flex 两种极线文件格式。它对几何和极线数据进行插值，以创建空气动力学计算点 (zaero)。

   • 输入参数 (para 字典)： 包含几何文件路径、几何设置编号、极线文件路径、空气动力学计算点的分布类型、以及各种插值方法等。

   • 方法：

     • chord_coordinate_system(z)：计算叶片上给定 z 位置的弦坐标系 (CCS)。
     • states_and_forces()：收集每个空气动力学计算点 (ACP) 的空气动力学状态和力。
     • geometry()：返回用于可视化的几何数据。
     • get_AoA_stall_margins(max_rel_thickness)：计算攻角失速裕度。
     • plot_stall_margins(...)：绘制失速裕度。

2. aero_calc_point 类：

   • 目的： 表示叶片上的单个空气动力学计算点。
   • 初始化 (__init__)： 存储初始和更新后的几何信息、翼型极线数据（攻角、升力系数、阻力系数、力矩系数），并初始化空气动力学状态变量（攻角、相对速度、力、力矩）。它还为升力系数、阻力系数和力矩系数创建插值函数。
   • 方法：
     • update_steady_aero_forces(rho)：根据空气密度 (rho) 计算 ACP 处的稳态空气动力学力和力矩。

3. 内部函数：

   • change_reference_curve(...)：修改现有空气动力学输入文件的参考曲线并保存新文件。
   • read_and_interpolate_HAWC2_polars(...)：从 HAWC2 格式文件读取和插值翼型极线。
   • read_and_interpolate_Flex_polars(...)：从 Flex 格式文件读取和插值翼型极线。

该文件提供了定义、初始化和计算风力涡轮机叶片空气动力学特性和力的核心功能，包括处理叶片几何、翼型极线和动态失速建模（尽管 unsteady_airfoil_dynamics 类被注释掉，表明它是一个正在进行的工作或未来功能）。它对于在 CASEStab 框架内模拟风力涡轮机的空气动力学行为至关重要。

casestab/casestab.py

主要负责风力涡轮机转子的稳态计算和变桨曲线调优。

主要组成部分：

1. rotor_models 类：

   • 目的： 计算轴对称转子变形的稳态，并提供变桨曲线调优功能。

   • 初始化 (__init__)：

     • 读取一个 JSON 格式的输入文件 (filename)，该文件定义了风力涡轮机模型的各种参数。
     • 根据 JSON 输入文件中的数据，定义并初始化模型的各个子结构（如轮毂 hub 和叶片 blade）。
     • 定义叶片的空气动力学参数（如几何文件、极线文件、插值方法等）。
     • 定义尾流模型 (wake) 和风模型 (wind)。
     • 定义轴对称转子模型 (rotor)。
     • 根据输入文件中的 deflection 标志设置是否包含变形计算。
     • 使用 model_assembler 模块中的 ma.model 类组装整个模型。
     • 如果提供了操作数据文件 (operation 字段)，则读取该文件并为每个操作点创建模型的深拷贝。

   • 方法：

     • steady_state_computation(ops_list=[])：对指定的操作点列表执行稳态计算。它会迭代每个操作点，设置相应的风速、转速和变桨角，然后调用模型的 compute_rotor_stationary_steady_state 方法来计算稳态。

     • save_steady_state_results(prefix='')：将稳态计算结果（包括结构变形和气动 BEM 结果）保存到文件中。文件命名会包含风速、变桨角和转速信息。

     • tune_pitch_curve(...)：这是一个核心功能，用于调优风力涡轮机的变桨曲线，以满足额定功率 (Prated)、推力限制 (Tlimit) 和失速裕度 (StallMargin) 等要求。

       • 它通过迭代的方式调整变桨角，每次迭代都会进行扰动运行和验证运行，以计算变桨角对功率、推力和失速裕度的梯度，并据此更新变桨角。
       • 支持绘制变桨、功率和推力曲线的迭代过程。
       • 最终将调优后的操作数据保存到 .opt 文件中，并附带调优原因和详细信息。

casestab.py 是 CASEStab 应用程序的入口点和主要控制器。它负责从配置文件中读取模型定义，构建复杂的风力涡轮机模型，执行稳态分析，并提供高级功能，如自动调优变桨曲线以优化性能和满足运行限制。这使得用户能够对风力涡轮机在不同运行条件下的行为进行详细的模拟和优化。

casestab/corotbeam_precompiled_functions.py

主要包含使用 Numba njit 装饰器进行即时 (JIT) 编译的函数，以提高性能。这些函数是用于处理共旋转梁单元（corotational beam elements）的底层数学和力学操作，这在结构动力学和有限元分析中非常常见。

主要组成部分和功能：

1. Numba 编译设置：

   • from numba import njit,float64,int32 和 from numba.pycc import CC：导入 Numba 库，用于 JIT 编译和提前编译 (AOT)。
   • cc = CC('corotbeam_precompiled_functions')：创建一个 CC 对象，用于将这些函数编译成一个可导入的模块，从而在 Python 解释器之外运行，进一步提高性能。
   • @njit(...) 装饰器：应用于每个函数，指示 Numba 对其进行编译，并指定输入参数的类型，以实现最佳性能。

2. 基本数学运算：

   • innerproduct(v1, v2)：计算两个三维向量的内积。
   • crossproduct(v1, v2)：计算两个三维向量的叉积。
   • skewmul(v, e)：计算向量 v 的反对称矩阵与向量 e 的乘积。
   • matvec33(A, b)：计算 3x3 矩阵 A 与 3x1 向量 b 的乘积。
   • matmul33(A, B)：计算两个 3x3 矩阵 A 和 B 的乘积。
   • matvec67(A, b)：计算 6x7 矩阵 A 与 7x1 向量 b 的乘积（用于特定于梁单元的转换）。
   • Skew(v)：返回给定向量 v 的反对称矩阵。

3. Rodrigues 旋转矩阵及其导数：

   • Rmat(v)：根据 Rodrigues 参数伪向量 v 计算 Rodrigues 旋转矩阵。这是一种表示三维旋转的有效方法。
   • dRmat(i, v)：计算 Rodrigues 旋转矩阵对第 i 个 Rodrigues 参数的一阶导数。
   • ddRmat(i, j, v)：计算 Rodrigues 旋转矩阵对第 i 和第 j 个 Rodrigues 参数的二阶导数。
   • RotateTriad(T, v)：使用 Rodrigues 旋转矩阵旋转一个三维正交基（triad）。
   • dRotateTriad(i, T, v)：计算旋转三维正交基对第 i 个 Rodrigues 参数的一阶导数。
   • ddRotateTriad(i, j, T, v)：计算旋转三维正交基对第 i 和第 j 个 Rodrigues 参数的二阶导数。

4. 辅助矩阵和索引：

   • Smat：用于 Rodrigues 旋转矩阵及其导数的常数辅助矩阵。
   • Imat：3x3 单位矩阵。
   • ij3sym, ij6sym, ij12sym：对称矩阵的索引映射，用于将对称矩阵的唯一元素映射到一维数组中，这在处理对称刚度矩阵和质量矩阵时很常见，可以节省存储和计算。

5. 梁单元运动学和动力学相关函数：

   • c_function(r)：一个辅助函数，可能用于多项式积分的系数计算。
   • ek_vT(v, k)：计算单位向量 e_k 与向量 v 的转置的乘积矩阵。
   • scalar_G_operator(...) 和 matrix_G_operator(...)：用于惯性状态计算的通用运算符，可能与质量矩阵的形成有关。
   • update_nodal_triads_and_position(...)：更新节点正交基和位置，并计算其一阶和二阶导数。这是共旋转梁单元公式中的关键步骤，它将全局位移转换为局部变形。
   • compute_element_triad_and_position(...)：计算单元正交基和位置，并计算其一阶和二阶导数。
   • update_local_nodal_rotations_elongation(...)：计算局部节点旋转（假定为小量）和伸长量。
   • update_first_derivative_local_nodal_rotations_elongation(...)：计算局部节点旋转和伸长量的一阶导数。
   • update_second_derivative_local_nodal_rotations_elongation(...)：计算局部节点旋转和伸长量的二阶导数。
   • update_element_deflection_subvectors_and_derivatives(...)：计算每个形函数阶次的变形子向量及其导数。
   • update_element_inertia(...)：更新单元惯性属性，包括质量中心和惯性张量及其导数。
   • update_forcing_point_position_and_moment_arm_vectors(...)：更新力作用点位置和力矩臂向量及其导数。
   • compute_element_total_force_matrix(...)：计算单元总力矩阵。
   • compute_element_total_moment_matrix(...)：计算单元总力矩矩阵。
   • compute_element_generalized_force_matrix(...)：计算单元广义力矩阵。
   • compute_element_stiffness_generalized_force_matrix(...)：从稳态广义力计算单元刚度矩阵。

该文件是 CASEStab 结构分析模块的性能关键部分。它提供了高效、预编译的函数，用于处理共旋转梁单元的复杂运动学和动力学计算，包括大位移和小应变。这些函数是构建和求解风力涡轮机结构模型的基础，确保了计算的准确性和效率。

casestab/corotbeam.py

实现了共旋转梁单元 (co-rotational beam element) 的公式，这是一种用于结构分析的非线性有限元方法，能够处理大位移和小应变。

主要组成部分和功能：

1. 辅助函数：

   • element_coordinate_system(r0, angle)：根据初始节点位置和平均旋转角度计算初始单元坐标系 (ECS)。
   • c_function(r)：一个辅助函数，可能用于多项式积分的系数计算（与 corotbeam_precompiled_functions.py 中的同名函数功能类似）。

2. corotbeam_substructure 类：

   • 目的： 表示由共旋转平衡梁单元描述的结构体（例如风力涡轮机叶片）。

   • 初始化 (__init__)：

     • 根据输入的节点分布 (znode) 和结构数据文件 (fname) 创建梁单元。它支持从 HAWC2 元素或通用结构文件读取数据。
     • 为每个梁单元创建 corotbeam_element_kinematics（运动学）和 equilibrium_beam_element（模型）实例。
     • 设置内部自由度向量 (q, dqdt)、内部力向量 (Fint) 和刚度矩阵 (K)。
     • 初始化惯性状态变量 (inertia)。
     • 处理与气动弹性耦合相关的力点和运动点数据。

   • 方法：

     • update_substructure()：更新所有梁单元的运动学状态，包括节点正交基、位置、局部变形及其导数。
     • update_elastic_internal_forces_and_stiffness()：计算内部弹性力和弹性刚度矩阵。
     • update_inertia()：更新子结构的惯性状态，包括质量中心、惯性张量及其导数。
     • update_node_position_and_rotation(inode)：获取指定节点的当前位置和旋转矩阵。
     • node_rotations()：计算所有节点的旋转伪向量。
     • initiate_aeroelastic_coupling(...)：初始化气动弹性耦合，识别作用在每个梁单元上的力点和运动点，并存储相关数据。
     • update_aeroelastic_coupling()：更新气动弹性耦合，包括从分布式力矩到节点力的转换矩阵，以及气动中心和共定位点的运动状态。
     • compute_rotation_peudo_vector()：计算旋转伪向量。

3. elements_created_from_file 类：

   • 目的： 从结构输入文件创建梁单元数据。
   • 初始化 (__init__)： 读取结构数据文件（支持各向同性截面和完整的 6x6 柔度矩阵格式），并根据节点分布创建梁单元的初始几何、长度、初始位置和初始单元坐标系。它还会对柔度矩阵和惯性参数进行多项式拟合。

4. equilibrium_beam_element 类：

   • 目的： 定义 Krenk 和 Couturier 平衡梁单元的形函数和刚度矩阵。

   • 初始化 (__init__)： 保存柔度矩阵 (C) 和单元长度 (l)。计算 H 矩阵及其逆矩阵 Hinv，以及 GT 矩阵。计算单元刚度矩阵 K 和局部单元刚度矩阵 Kl。计算应变函数和形函数的多项式分量。插入惯性参数并计算单元质量。

   • 方法：

     • shape_function_matrix(zeta)：在给定 zeta（无量纲单元坐标）处评估完整的 6x12 形函数矩阵。
     • local_shape_function_matrix(zeta)：在给定 zeta 处评估局部 6x7 形函数矩阵。
     • local_beam_element_deflection(ql, zeta)：计算多个 zeta 点处的局部梁单元变形和旋转。
     • print_element_properties()：打印单元属性（长度和质量）。

5. corotbeam_element_kinematics 类：

   • 目的： 包含共旋转梁单元的运动学信息。

   • 初始化 (__init__)： 存储初始长度、位置和方向。初始化动态值（当前长度、节点正交基、单元坐标系等）及其一阶和二阶导数。

   • 方法： 这些方法主要调用 corotbeam_precompiled_functions 中对应的预编译函数，以更新各种运动学量及其导数，例如：

     • update_nodal_triads_and_position(q)
     • compute_element_triad_and_position()
     • update_local_nodal_rotations_elongation()
     • update_first_derivative_local_nodal_rotations_elongation()
     • update_second_derivative_local_nodal_rotations_elongation()
     • update_element_deflection_subvectors_and_derivatives(Nl)

6. acp_motion 类：

   • 目的： 包含单元运动的耦合反馈信息。
   • 初始化 (__init__)： 保存关联的单元编号、无量纲单元坐标、运动点在单元坐标系中的位置和力坐标系。

7. forced_element 类：

   • 目的： 处理单元节点力与外部力分布的耦合。

   • 初始化 (__init__)： 保存力点信息（索引、无量纲坐标、位置、方向），计算插值系数。

   • 方法： 这些方法主要调用 corotbeam_precompiled_functions 中对应的预编译函数，以计算各种力矩阵及其导数，例如：

     • compute_element_total_force_matrix(elem_model)
     • compute_element_total_moment_matrix(elem_model)
     • compute_element_generalized_force_matrix(elem_model)
     • compute_element_stiffness_generalized_force_matrix(elem_model)
     • update_forcing_point_position_and_moment_arm_vectors(elem_state, elem_model)

corotbeam.py 文件是 CASEStab 中实现共旋转梁单元有限元分析的核心。它定义了梁单元的运动学、静力学和动力学行为，并提供了从结构输入文件创建单元、更新单元状态、计算内部力和刚度、更新惯性以及处理气动弹性耦合的完整框架。通过与预编译函数的结合，它旨在提供高效且准确的结构分析能力，特别适用于风力涡轮机叶片等柔性结构的建模。

casestab/generic_model_components.py

定义了用于构建和管理复杂系统（如风力涡轮机模型）的通用模型组件，特别是与惯性、气动弹性耦合和轴承相关的组件。

主要组成部分和功能：

1. substructure_inertia 类：

   • 目的： 封装子结构的惯性状态。在有限元分析中，子结构通常指模型中的一个独立部分（例如，风力涡轮机中的叶片、轮毂或塔架）。

   • 初始化 (__init__)：

     • ndofs：子结构的自由度数量。

     • 初始化各种惯性相关的属性，包括：

       • mass：总质量。
       • Ibase：基坐标系中的旋转惯量张量。
       • M11：局部质量矩阵。
       • G11：局部陀螺和离心刚度矩阵（初始化为零，可能在后续计算中填充）。
       • Kc11：离心刚度矩阵（初始化为零）。
       • Fc1：离心力向量（初始化为零）。
       • H111：局部非线性陀螺矩阵（初始化为零）。
       • m_rcg：质量乘以质心向量。
       • m_drcg_dqi、m_ddrcg_dqidqj：质心向量对自由度的一阶和二阶导数。
       • Abase_i、Abase1_ij、Abase2_ij：用于惯性计算的辅助矩阵及其导数。
       • jcol_nonzero_irow：非零惯性矩阵元素的列索引和行索引（用于稀疏矩阵操作）。

     • nsym 和 ijsym：用于处理对称矩阵的辅助变量，通过 math_functions.generate_ijNsym 生成对称索引映射。

   • 方法：

     • reset_inertia()：将所有惯性相关的属性重置为零，以便在每次更新或迭代之前进行清理。
     • compute_local_centrifugal_forces_and_matrix(...)：计算局部离心力和离心刚度矩阵。这个方法调用了 model_precompiled_functions 中的预编译函数来执行实际的计算，以提高性能。

2. initiate_acp_motion_state 类：

   • 目的： 封装气动计算点 (ACP) 的运动状态信息。ACP 是叶片上用于计算气动力的特定点。

   • 初始化 (__init__)：

     • idofs：与该 ACP 相关的自由度索引。

     • ndofs：相关自由度的数量。

     • 初始化各种运动状态相关的属性，包括：

       • rtp：扭转点 (Torsional Point) 的位置向量。
       • rcp：共定位点 (Coallocation Point) 的位置向量。
       • Ec：CCS（弦坐标系）在子结构坐标系中的方向矩阵。
       • drtp_dqi、drcp_dqi、dec1_dqi、dec2_dqi：这些位置和方向对自由度的一阶导数。
       • ddrtp_dqidqj、ddrcp_dqidqj：这些位置对自由度的二阶导数。

3. bearing 类：

   • 目的： 定义不同类型的轴承行为。轴承用于连接子结构并定义它们之间的相对运动。

   • 初始化 (__init__)：

     • bear_flag：指示是否定义了轴承。
     • bear_dof：轴承的自由度。
     • 根据 bearing_text 字符串解析轴承类型（free、constant_angle、constant_speed）、轴、名称和单位，并创建相应的内部状态对象。

   • 内部类：

     • free_bearing：表示自由轴承，允许绕指定轴自由旋转。
     • constant_angle_bearing：表示固定角度轴承，保持指定轴上的角度不变。
     • constant_speed_bearing：表示恒定速度轴承，以恒定角速度绕指定轴旋转。
     • 这些内部类都包含 B（旋转矩阵）、dB（旋转矩阵的一阶导数）、ddB（旋转矩阵的二阶导数）以及 BTdB 和 BTddB（辅助矩阵）等属性，并提供 update(t) 方法来更新其状态。

generic_model_components.py 文件提供了构建复杂多体动力学模型所需的通用、可重用的组件。它抽象了惯性属性、气动弹性耦合点运动学和轴承行为的细节，使得在 CASEStab 中可以更模块化地组装和分析风力涡轮机等系统。通过将计算密集型操作委托给预编译函数，它也确保了性能。

casestab/HAWC2_blade_translator.py

主要负责读取 HAWC2（一种风力涡轮机仿真软件）模型文件，并将其结构和气动数据转换为 CASEStab 所需的格式。这对于在 CASEStab 中使用 HAWC2 定义的风力涡轮机模型进行分析至关重要。

主要组成部分和功能：

1. HAWC2_elements 类：

   • 目的： 从 HAWC2 模型文件读取并插入结构元素数据。

   • 初始化 (__init__)：

     • 读取 HAWC2 结构文件 (fname)，并根据 mbdy_name 提取特定主体的元素位置 (elempos)、元素旋转 (elemrot) 和元素数据 (elemdat)。
     • 根据 HAWC2 的各向同性梁单元数据（例如，弹性模量 E、剪切模量 G、惯性矩 Ix, Iy, Iz、剪切修正因子 kx, ky 等），使用 isotropic_to_6x6_compliance_matrix 函数计算并存储每个元素的 6x6 柔度矩阵 (Cs)。
     • 存储每个元素的长度 (l)、z 坐标 (z)、节点位置 (pos)、初始节点位置 (r0s) 和初始单元坐标系 (E0s)。
     • 计算并存储每个元素的惯性参数 (iner_pars)。

2. chord_coordinate_system_in_HAWC2(c2def, z) 函数：

   • 目的： 根据 HAWC2 的 c2_def 定义和给定的 z 坐标，计算叶片的弦坐标系 (CCS)。

   • 功能：

     • 使用 Akima 插值器计算半弦曲线的斜率 (xp, yp)。
     • 插值气动扭转角 (twist)。
     • 迭代每个截面，计算其在叶片坐标系中的 CCS 矩阵 (Ec)，以及从叶片坐标系到 CCS 的有限旋转伪向量 (phi)。
     • 返回半弦点的位置 (rc2)、CCS 矩阵 (ccs) 和旋转伪向量 (phi)。

3. read_HAWC2_main_body 类：

   • 目的： 读取 HAWC2 HTC 文件，提取指定主体的结构设置编号和 c2_def 数据。

   • 初始化 (__init__)：

     • 解析 HTC 文件 (htcfile)，查找与 mbdy_name 匹配的主体。
     • 提取 timoschenko_input 部分的 set 和 subset 编号。
     • 提取 c2_def 部分的几何定义数据 (c2def)。
     • 读取结构文件 (st_filename)，并根据提取的 set_no 和 subset_no 获取结构截面数据 (stset)。

4. read_HAWC2_ae_set 类：

   • 目的： 读取 HAWC2 HTC 文件，提取指定叶片的气动设置编号和气动数据。

   • 初始化 (__init__)：

     • 解析 HTC 文件 (htcfile)，查找 aero 部分的 ae_sets 编号。
     • 读取气动文件 (ae_filename)，并根据提取的 set_no 获取气动截面数据 (aeset)。

   • 方法：

     • pc_set_nr(z)：根据 z 坐标插值获取翼型极线设置编号。

5. read_HAWC2_pc_file 类：

   • 目的： 读取 HAWC2 翼型极线文件。

   • 初始化 (__init__)：

     • 读取翼型极线文件 (pc_filename)。
     • 解析文件内容，将翼型极线数据组织成字典 pcsets，其中包含不同设置编号和厚度下的极线数据。

6. translate_HAWC2_blade_model(...) 函数：

   • 目的： 将 HAWC2 叶片输入模型文件（HTC、气动文件、极线文件、结构文件）转换为 CASEStab 所需的 SDU 格式。

   • 功能：

     • 使用上述 read_HAWC2_main_body、read_HAWC2_ae_set 和 read_HAWC2_pc_file 类读取 HAWC2 模型的原始数据。

     • 创建结构输入文件：

       • 根据 HAWC2 结构数据和 c2_def 中的几何信息，计算 CASEStab 结构文件所需的各项参数（例如，参考点位置、角度、质量属性、惯性半径、弹性轴、剪切中心、刚度属性等）。
       • 支持 ISO（各向同性）和 6x6（完整柔度矩阵）两种输出格式。
       • 将转换后的结构数据保存到名为 stru_ + st_format + .dat 的文件中。

     • 创建气动输入文件：

       • 根据 HAWC2 气动数据和 c2_def 中的几何信息，计算 CASEStab 气动文件所需的各项参数（例如，弦坐标系位置、旋转伪向量、弦长、相对厚度、气动压力中心位置、极线设置编号等）。
       • 将转换后的气动数据保存到名为 aero.dat 的文件中。

     • 返回转换后的结构和气动数据数组。

7. HAWCStab2_blade 类：

   • 目的： 读取 HAWCStab2 结果文件，提取气动中心位置和 CCS 方向。

   • 初始化 (__init__)：

     • 解析 HAWCStab2 结果文件 (fname)。
     • 提取气动中心位置 (rac) 和 CCS 方向矩阵 (Ec)。
     • 提取弦长和厚度信息 (ct)。

casestab/math_functions.py

提供了一系列用于三维旋转、向量操作和数据插值的数学函数，这些函数在结构动力学和气动弹性分析中是基础性的。其中一些函数使用 Numba njit 装饰器进行即时 (JIT) 编译，以提高性能。

主要组成部分和功能：

1. 常量：

   • Smat：用于 Rodrigues 旋转矩阵及其导数的常数辅助矩阵（反对称矩阵）。
   • Imat：3x3 单位矩阵。

2. 旋转矩阵函数：

   • rotmat(v, phi)：根据旋转轴单位向量 v 和旋转角度 phi 计算三维旋转矩阵（使用 Rodrigues 公式）。该函数使用 njit 编译。
   • rotmat_from_pseudovec(p)：根据 Rodrigues 参数伪向量 p 计算旋转矩阵。如果伪向量的模很小，则返回单位矩阵。
   • rotmat_to_quaternion(R)：将旋转矩阵 R 转换为四元数。
   • quaternion_to_vector_and_angle(q)：将四元数 q 转换为旋转轴单位向量 v 和旋转角度 phi。
   • interpolate_rotmat(T, Q, ratio)：在两个三维正交基 T 和 Q 之间进行旋转插值。
   • R1(phi)、R2(phi)、R3(phi)：分别计算绕第一、第二、第三轴旋转 phi 角度的旋转矩阵。
   • Ri(phi, iaxis)：计算绕指定轴 iaxis（可以是负数，表示反方向）旋转 phi 角度的旋转矩阵。
   • Skew(v)：返回给定向量 v 的反对称矩阵。该函数使用 njit 编译。
   • deskew(S)：从反对称矩阵 S 中提取原始向量。
   • Rmat(v)：根据 Rodrigues 参数伪向量 v 计算 Rodrigues 旋转矩阵。
   • dRmat(i, v)：计算 Rodrigues 旋转矩阵对第 i 个 Rodrigues 参数的一阶导数。
   • ddRmat(i, j, v)：计算 Rodrigues 旋转矩阵对第 i 和第 j 个 Rodrigues 参数的二阶导数。
   • pseudo_vector_from_Rodrigues(q)：从 Rodrigues 参数计算伪向量。
   • small_rotation_pseudo_vector_from_Rodrigues(q)：从 Rodrigues 参数计算小旋转伪向量（处理小角度旋转的特殊情况）。

3. 向量和矩阵操作：

   • innerproduct(v1, v2)：计算两个三维向量的内积。该函数使用 njit 编译。
   • transposed_innerproduct(v1, v2)：计算向量 v1 与向量 v2 转置的乘积矩阵。
   • inner_matrix_product(A1, A2)：计算两个 3x3 矩阵的内积（Frobenius 内积）。
   • vector_length(v)：计算向量的长度。
   • unit_vector(v)：将向量归一化为单位向量。
   • crossproduct(v1, v2)：计算两个三维向量的叉积。

4. 索引和插值函数：

   • generate_ijNsym(N)：生成一个 N x N 对称矩阵的索引映射，用于将对称矩阵的唯一元素映射到一维数组中。

   • curve_interpolate 类：

     • 目的： 提供不同类型的曲线插值功能。
     • 初始化 (__init__)： 根据 itype（pchip、akima 或 linear）选择相应的 SciPy 插值器来创建函数 (fcn) 和其导数 (der)。

   • piecewise_linear_function(x, xbreak, coeffs)：实现分段线性函数。

   • quick_interpolation_periodic_function 类：

     • 目的： 为等间距数据点提供快速周期性线性插值。

     • 初始化 (__init__)： 存储数据点 x 和 y，以及步长 dx 和范围 xrange。检查函数是否周期性。

     • 方法：

       • fcn(x)：计算给定 x 值的插值结果。
       • der(x)：计算给定 x 值的导数。

math_functions.py 文件是 CASEStab 项目的数学工具箱。它提供了处理三维几何、旋转和数据插值所需的基本和高级数学运算。通过使用 Numba 进行编译，这些函数旨在提供高性能的计算能力，这对于风力涡轮机结构和气动弹性分析中的复杂数值模拟至关重要。

casestab/model_assembler.py

负责将风力涡轮机模型的各个部分（子结构、叶片、转子、尾流、风）组装成一个完整的系统模型，并提供进行稳态分析和结果可视化的功能。

主要组成部分和功能：

1. substructure 类：

   • 目的： 封装不同类型的子结构模型（例如刚性向量或共旋转梁）。

   • 初始化 (__init__)：

     • 保存子结构参数，包括唯一编号 (isubs)、名称 (name)、连接信息 (isubs_connection, inode_connection) 和初始方向 (Sbase)。
     • 处理轴承 (bearing) 的定义。
     • 根据 para['type'] 创建具体的子结构实例（rigidbody.rigidbody_substructure 或 corotbeam.corotbeam_substructure）。
     • 初始化与子结构在地面固定坐标系中的位置、方向及其时间导数和自由度导数相关的变量。

   • 方法：

     • compute_modes()：计算子结构的结构模态（仅适用于柔性梁子结构）。它会设置质量和刚度矩阵，求解特征值问题，并对模态进行排序和命名。
     • plot_substructure_modes(...)：绘制子结构的模态形状。
     • create_data_for_deflection_state()：创建用于结构变形状态输出的数据，包括叶片坐标系和转子平面内的变形以及旋转。

2. rotor 类：

   • 目的： 封装不同类型的转子模型（例如轴对称转子）。

   • 初始化 (__init__)：

     • 关联叶片 (blades) 和子结构 (substructures)。
     • 检查叶片的气动计算点 (ACP) 数量是否一致。
     • 关联风模型 (wind)。
     • 保存转子中心的连接点信息。
     • 计算转子平面在地面固定坐标系中的法向量 (nvec)。
     • 初始化转子速度 (omega)、中心位置 (rc0) 和方向矩阵 (Rc0)。
     • 初始化功率 (power)、扭矩 (torque)、推力 (thrust) 和功率系数 (CP)。
     • 计算 ACP 的径向位置 (radii) 和实度 (sigma)。
     • 创建尾流模型 (wake_model.wake)。

   • 方法：

     • update_rotor_kinematic_state()：更新转子的运动学状态，包括中心位置、方向和转子速度。
     • update_steady_blade_forces()：计算每个叶片上的稳态气动力，假设尾流处于平衡状态。这涉及到使用 scipy.optimize.root 求解 BEM（叶素动量理论）方程。
     • create_data_for_BEM_results(iblade=0)：创建用于 BEM 结果输出的数据，包括气动计算点的各种气动参数和力。

3. model 类：

   • 目的： 组装整个风力涡轮机模型。

   • 初始化 (__init__)：

     • 根据 subs_para_set 创建所有子结构实例，并进行初步更新。
     • 处理子结构之间的连接链，并识别具有恒定速度轴承的支撑子结构。
     • 根据 blade_para_set 创建所有叶片实例，并将叶片与相应的子结构关联起来，初始化气动弹性耦合。
     • 初始化子结构在地面固定坐标系中的初始位置和方向。
     • 根据 wind_para_set 创建风场模型。
     • 根据 rotor_para_set 创建转子模型，并将叶片、子结构、尾流和风模型传递给转子实例。

   • 方法：

     • update_all_substructures(t)：在给定时间 t 更新所有子结构的状态，包括轴承、弹性内力、刚度、惯性状态、离心力和气动弹性耦合。
     • update_steady_state_all_rotors()：更新所有转子的稳态气动力。
     • compute_rotor_stationary_steady_state(irotor, Rnorm_limit=1.0, include_deform=True)：计算指定转子的稳态（静态外部力和内力之间的平衡）。这是一个迭代过程，它会反复更新气动力和结构变形，直到收敛。
     • compute_substructure_steady_state_deformation(isubs, Rnorm_limit=1.0)：计算单个子结构的稳态变形。
     • plot_aerodynamic_points(isubs, ielem, iblade)：绘制气动计算点和结构节点的位置，用于可视化。
     • plot_input_data_blade(iblade, fn='')：绘制叶片的气动和结构输入数据，包括气动中心、质心、弹性轴、剪切中心等的位置。

model_assembler.py 文件是 CASEStab 仿真框架的集成器。它将所有独立的物理模型（结构、气动、风、尾流、轴承）组合成一个连贯的系统，并提供了执行稳态分析和可视化结果的接口。这使得 CASEStab 能够对风力涡轮机进行全面的多物理场仿真。

casestab/model_precompiled_functions.py

包含使用 Numba njit 装饰器进行即时 (JIT) 编译的函数，以及使用 numba.pycc.CC 进行提前 (AOT) 编译的设置。这些函数主要用于计算模型中的离心力、陀螺矩阵和离心刚度矩阵，这些是多体动力学仿真中的关键组成部分。

主要组成部分和功能：

1. Numba 编译设置：

   • from numba import njit,float64 和 from numba.pycc import CC：导入 Numba 库，用于 JIT 编译和 AOT 编译。
   • cc = CC('model_precompiled_functions')：创建一个 CC 对象，用于将这些函数编译成一个可导入的模块。
   • @njit(...) 装饰器：应用于每个函数，指示 Numba 对其进行编译，并指定输入参数的类型，以实现最佳性能。

2. 基本数学运算：

   • innerproduct(v1, v2)：计算两个三维向量的内积。
   • inner_matrix_product(A1, A2)：计算两个 3x3 矩阵的内积（Frobenius 内积）。

3. 离心力、陀螺矩阵和离心刚度矩阵计算：

   • compute_local_centrifugal_forces_and_matrix(...)：

     • 目的： 计算局部离心力向量 (Fc)、局部陀螺矩阵 (Gc) 和局部离心刚度矩阵 (Kc)。这些量在处理旋转系统中的惯性效应时非常重要。

     • 输入参数：

       • ndofs：自由度数量。
       • jcol_nonzero_irow：非零元素的列索引和行索引（用于稀疏矩阵操作）。
       • ijsym：对称矩阵的索引映射。
       • m_drcg_dqi、m_ddrcg_dqidqj：质心向量对自由度的一阶和二阶导数。
       • Abase_i、Abase1_ij、Abase2_ij：用于惯性计算的辅助矩阵及其导数。
       • R0Tddr0：基坐标系中位置向量的二阶时间导数。
       • R0TddR0：基坐标系中旋转矩阵的二阶时间导数。
       • R0TdR0：基坐标系中旋转矩阵的一阶时间导数。

     • 计算逻辑： 该函数通过遍历自由度，利用输入的各种导数和辅助矩阵，计算 Fc、Gc 和 Kc 的各个分量。它考虑了对称性，以填充矩阵的上下三角形。

model_precompiled_functions.py 文件是 CASEStab 动力学分析模块的性能关键部分。它提供了高效、预编译的函数，用于计算旋转系统中由惯性引起的力和刚度效应。这些函数是构建和求解风力涡轮机动力学方程的基础，确保了计算的准确性和效率。通过提前编译，这些函数可以在 Python 环境之外运行，进一步优化了性能。

casestab/rigidbody.py

定义了一个简单的刚体子结构模型，主要用于表示模型中不发生变形的部分，例如风力涡轮机的轮毂或塔基。

主要组成部分和功能：

1. rigidbody_substructure 类：

   • 目的： 表示一个刚体子结构。

   • 初始化 (__init__)：

     • type：设置为 'rigid'，表明这是一个刚体。
     • ndofs：自由度数量，初始化为 0，因为刚体子结构本身不引入内部变形自由度。其运动由其连接点（在 model_assembler 中处理）决定。
     • rnode：存储节点的定义（位置）。
     • aero_part：标志，指示该子结构是否与气动部分相关，初始化为 False。
     • iblade：关联的叶片编号，初始化为 -1。
     • inertia：初始化一个 gmc.substructure_inertia 对象，但由于是刚体且 ndofs 为 0，其内部惯性矩阵和向量将是空的或零。

   • 方法：

     • update_node_position_and_rotation(inode)：返回指定节点的当前位置和单位旋转矩阵（因为是刚体，内部不发生旋转）。
     • update_substructure()：一个空方法 (pass)，因为刚体子结构内部不发生变形，不需要更新其内部状态。
     • update_elastic_internal_forces_and_stiffness()：一个空方法 (pass)，因为刚体没有弹性力或刚度。
     • update_inertia()：一个空方法 (pass)，因为刚体子结构的惯性通常在更高级别的模型中处理，或者其内部惯性状态不随时间变化。

rigidbody.py 文件提供了一个轻量级的刚体模型，用于 CASEStab 中不需要详细结构变形分析的部分。它简化了模型，提高了计算效率，因为它避免了对这些部分进行复杂的有限元计算。

casestab/timoshenko_beam_section.py

主要提供了用于处理 Timoshenko 梁截面属性的函数，特别是将截面属性转换为 6x6 柔度矩阵，以及在不同参考点之间转换矩阵。

主要组成部分和功能：

1. transform_reference_point_of_matrix(S, dx, dy, theta) 函数：

   • 目的： 转换一个 6x6 的位移和旋转矩阵（例如柔度矩阵或刚度矩阵）的参考点。

   • 输入：

     • S：要转换的矩阵，定义在系统 2 中。
     • dx, dy：系统 2 原点在系统 1 x、y 轴上的位置。
     • theta：系统 2 相对于系统 1 的旋转角度。

   • 输出： 在系统 1 中定义的旋转后的矩阵。

   • 功能： 该函数构建了两个转换矩阵 T12 和 T21，它们包含了平移和旋转效应，然后通过 T21 @ S @ T12 的形式进行矩阵转换。

2. isotropic_to_6x6_compliance_matrix(rref, rea, rsc, theta, E, G, A, Ix, Iy, Iz, kx, ky) 函数：

   • 目的： 将传统的 Timoshenko 梁截面属性（各向同性材料）转换为 6x6 的柔度矩阵。

   • 输入：

     • rref：参考点在截面参考系中的 2D 位置。
     • rea：弹性轴（弯曲形心）在截面参考系中的 2D 位置。
     • rsc：剪切中心在截面参考系中的 2D 位置。
     • theta：主轴在截面参考系中的结构扭转角。
     • E：平均弹性杨氏模量。
     • G：平均剪切模量。
     • A：截面面积。
     • Ix：绕主 x 轴的弯曲惯性矩。
     • Iy：绕主 y 轴的弯曲惯性矩。
     • Iz：扭转惯性矩。
     • kx, ky：主 x 和 y 方向的剪切修正因子。

   • 输出： 6x6 的截面柔度矩阵 C。该矩阵将广义力向量（两个面内剪力、轴向正向力、两个弯矩和扭转力矩）与广义应变向量（两个面内剪切应变、轴向伸长应变、两个面外曲率和截面扭转率）关联起来。所有力和位移都定义在参考点。

   • 功能：

     • 首先计算弹性轴 (EA) 到剪切中心 (SC) 在主轴坐标系中的距离 (ex, ey)。
     • 根据 Hodges (2006) 的理论，在主轴坐标系中构建 6x6 的刚度矩阵 S。
     • 使用 transform_reference_point_of_matrix 函数将刚度矩阵 S 从主轴坐标系转换到参考点坐标系。
     • 最后，通过对刚度矩阵求逆来计算柔度矩阵 C。

timoshenko_beam_section.py 文件在 CASEStab 中扮演着关键角色，它使得程序能够处理和转换不同表示形式的梁截面属性。这对于构建准确的梁单元模型至关重要，因为它允许用户输入常见的工程参数（如杨氏模量、惯性矩等），并将其转换为有限元公式所需的柔度矩阵形式，同时还能处理不同参考点之间的转换。

casestab/wake_model.py

主要负责风力涡轮机转子上的尾流模型，用于计算诱导速度。

主要组成部分和功能：

1. wake 类：

   • 目的： 根据参数选择并实例化具体的尾流模型。

   • 初始化 (__init__)：

     • 保存参数 para 和尾流半径 wake_radii。
     • 根据 para['type'] 的值，实例化 no_induction（无诱导）或 axissym_induction（轴对称诱导）模型。

2. no_induction 类：

   • 目的： 表示没有诱导速度的尾流模型（即，不考虑尾流对气动力的影响）。

   • 初始化 (__init__)：

     • 初始化轴向诱导因子 a、切向诱导因子 ap、诱导速度 vi 及其时间导数 dvidt 为零。
     • 保存半径 radii、叶尖修正标志 tip_correction、叶片数量 nblades 和转子半径 R。
     • 初始化推力系数 CT、扭矩系数 CQ、入流角正弦 sinphi、局部叶尖速比 local_TSR 和叶尖损失因子 ftip 为零或一。

   • momentum_balance_point 内部类：

     • 目的： 在没有诱导的情况下，计算单个气动计算点 (ACP) 的动量平衡。

     • 初始化 (__init__)：

       • 计算局部叶尖速比 local_TSR。
       • 计算 ACP 处的相对速度分量 vx0, vy0。
       • 计算与推力 (aTx, aTy) 和扭矩 (aQx, aQy) 相关的气动系数。
       • 计算并保存升力系数 CL、阻力系数 CD 和力矩系数 CM。
       • 计算并保存推力系数 CT 和扭矩系数 CQ。

     • f(self, x) 和 fprime(self, x) 方法： 在此 no_induction 模型中，这些方法返回零，表示没有需要求解的非线性方程。

3. axissym_induction 类：

   • 目的： 表示轴对称诱导尾流模型，这是更实际的风力涡轮机气动模型。

   • 初始化 (__init__)：

     • 初始化轴向诱导因子 a、切向诱导因子 ap、诱导速度 vi 及其时间导数 dvidt。
     • 保存半径 radii、叶尖修正标志 tip_correction、叶片数量 nblades 和转子半径 R。
     • 根据 para['a_of_CT_model'] 选择轴向诱导因子与推力系数的关系模型（例如 HAWC2_a_of_CT）。
     • 初始化推力系数 CT、扭矩系数 CQ、入流角正弦 sinphi、局部叶尖速比 local_TSR 和叶尖损失因子 ftip。

   • momentum_balance_point 内部类：

     • 目的： 在轴对称诱导下，计算单个气动计算点 (ACP) 的动量平衡。

     • 初始化 (__init__)：

       • 计算局部叶尖速比 local_TSR。
       • 计算相对速度分量 vx0, vx1, vx2, vy0, vy1, vy2，这些分量依赖于诱导因子 a 和 ap。
       • 计算与推力 (aTx, aTy) 和扭矩 (aQx, aQy) 相关的气动系数。
       • 计算叶尖损失因子 beta。
       • 实例化 HAWC2_a_of_CT 和 HAWC2_ap_of_CQ 模型。

     • f(self, x) 方法： 定义了需要求解的非线性方程组。它根据当前的 a 和 ap 值计算相对速度、攻角、升阻力系数、推力系数 CT 和扭矩系数 CQ，并应用叶尖修正。然后，它返回两个残差值，表示 a 和 ap 是否满足动量平衡方程。

     • fprime(self, x) 方法： 定义了 f 方法的雅可比矩阵，用于非线性方程组的求解（例如牛顿法）。它计算 CT 和 CQ 对 a 和 ap 的偏导数，以及叶尖损失因子对 a 和 ap 的偏导数。

4. HAWC2_a_of_CT 类：

   • 目的： 实现了 HAWC2 中轴向诱导因子 a 与推力系数 CT 之间的分段多项式关系。

   • 初始化 (__init__)： 定义了多项式系数 k1, k2, k3 和分段点 CT1, CT2。

   • 方法：

     • inner_fcn(CT)：计算内部多项式函数。
     • inner_der(CT)：计算内部多项式函数的导数。
     • fcn(CT)：根据 CT 值返回轴向诱导因子 a，考虑了分段线性限制。
     • der(CT)：计算 a 对 CT 的导数。

5. HAWC2_ap_of_CQ 类：

   • 目的： 实现了 HAWC2 中切向诱导因子 ap 与扭矩系数 CQ 和轴向诱导因子 a 之间的关系。

   • 初始化 (__init__)： 根据局部叶尖速比 local_TSR 计算一个因子。

   • 方法：

     • fcn(CQ, a)：根据 CQ 和 a 返回切向诱导因子 ap，考虑了 a 大于 0.9 的特殊情况。
     • der(CQ, a)：计算 ap 对 CQ 和 a 的偏导数。

wake_model.py 文件是 CASEStab 中风力涡轮机气动弹性分析的关键组成部分。它提供了不同复杂程度的尾流模型，特别是轴对称诱导模型，该模型通过求解非线性方程组来计算叶片上的诱导速度和气动力。这些模型对于准确预测风力涡轮机的性能和载荷至关重要。

casestab/wind_model.py

主要负责定义风力涡轮机仿真中使用的风场模型。

主要组成部分和功能：

1. wind 类：

   • 目的： 根据参数选择并实例化具体的风场模型。

   • 初始化 (__init__)：

     • 保存参数 para。
     • 保存空气密度 rho，这是进行气动计算所必需的。
     • 根据 para['windtype'] 的值，实例化具体的风场查找对象。目前只支持 uniform（均匀风）。

2. uniform_wind 类：

   • 目的： 表示一个均匀风场，即风速在整个空间和时间上都是恒定的。

   • 初始化 (__init__)：

     • 保存平均风速 umean。

   • uvw_at_xyzt(self, pos, t) 方法：

     • 目的： 在给定位置 pos (x, y, z) 和时间 t 处返回风速向量 (u, v, w)。
     • 功能： 对于均匀风场，它返回一个固定方向（在此实现中是 y 方向）和大小为 umean 的风速向量。其他分量为零。

wind_model.py 文件为 CASEStab 仿真提供了风场输入。它允许用户定义不同类型的风场（目前仅支持均匀风），并提供一个接口来查询在模型中任意位置和时间点的风速。这对于计算风力涡轮机上的气动力和进行气动弹性分析至关重要。

casedamp包分析

casetoolbox/casedamp/
├── __init__.py
├── casedamp_precompiled_functions.py
├── casedamp.py
└── test/

casedamp/casedamp_precompiled_functions.py

主要包含使用 Numba njit 装饰器进行即时 (JIT) 编译的函数，以及使用 numba.pycc.CC 进行提前 (AOT) 编译的设置。这些函数专门用于计算风力涡轮机叶片的气动阻尼系数。

主要组成部分和功能：

1. Numba 编译设置：

   • from numba.pycc import CC：导入 Numba 库的提前编译模块。
   • cc = CC('casedamp_precompiled_functions')：创建一个 CC 对象，用于将这些函数编译成一个可导入的模块，从而在 Python 解释器之外运行，进一步提高性能。
   • @cc.export(...) 装饰器：用于将函数导出为编译模块的一部分，并指定输入参数的类型，以实现最佳性能。

2. 阻尼系数计算函数：

   • compute_damping_terms(aoas, psis, clcd_clpcdp)：

     • 目的： 计算气动阻尼的四个主要分量：W_tran1 (平动阻尼1)、W_tran2 (平动阻尼2)、W_tors1 (扭转阻尼1) 和 W_tors2 (扭转阻尼2)。

     • 输入参数：

       • aoas：攻角数组。
       • psis：方位角数组。
       • clcd_clpcdp：包含升力系数 (cl)、阻力系数 (cd) 及其对攻角的导数 (clp, cdp) 的数组。

     • 计算逻辑： 该函数遍历所有攻角和方位角的组合，根据气动系数及其导数，以及角度 v = psi + aoa 的三角函数，计算出四个阻尼分量。这些分量是气动阻尼模型的基础。

   • compute_damping_eta(ured, gama, beta, phi, W_tran1, W_tran2, W_tors1, W_tors2)：

     • 目的： 计算总的气动阻尼系数 eta。

     • 输入参数：

       • ured：无量纲减缩速度。
       • gama：一个与叶片几何相关的参数。
       • beta：一个与叶片几何相关的参数。
       • phi：入流角。
       • W_tran1, W_tran2, W_tors1, W_tors2：由 compute_damping_terms 函数计算出的阻尼分量。

     • 计算逻辑： 该函数将四个基本阻尼分量与无量纲速度、几何参数和入流角结合起来，计算出最终的综合气动阻尼系数 eta。

casedamp_precompiled_functions.py 文件是 CASEDamp 模块的性能关键部分。它提供了高效、预编译的函数，用于执行气动阻尼计算中涉及的复杂数学运算。这些函数是评估风力涡轮机叶片气动阻尼特性的核心，对于理解和预测叶片的动态响应至关重要。通过提前编译，这些函数可以在 Python 环境之外运行，从而显著提高了计算效率。

casedamp/casedamp.py

提供了一个交互式工具，用于分析和编辑翼型极线（升力系数 CL 和阻力系数 CD 曲线），并计算其气动阻尼特性。

主要组成部分和功能：

1. curve_interpolate 类：

   • 目的： 提供不同类型的曲线插值功能，特别是针对攻角 (AoA) 和升力/阻力系数。
   • 初始化 (__init__)： 根据 itype（pchip、akima 或 linear）选择相应的 SciPy 插值器来创建函数 (fcn) 和其导数 (der)。

2. airfoil_class 类：

   • 目的： 封装单个翼型的原始极线数据和厚度信息。
   • 初始化 (__init__)： 存储原始数据 (data) 和厚度 (thickness)，并初始化 fcn 为 None（稍后由 read_HAWC2_pc_file 或 read_Flex_pro_file 填充）。

3. read_HAWC2_pc_file 类：

   • 目的： 从 HAWC2 格式的翼型极线文件读取数据。
   • 初始化 (__init__)： 解析 HAWC2 极线文件，将数据组织成 set 字典，其中包含不同设置编号和翼型编号下的 airfoil_class 实例，并为每个翼型创建 curve_interpolate 对象。
   • save_file() 方法： 将当前编辑的极线数据保存回 HAWC2 格式的新文件。

4. read_Flex_pro_file 类：

   • 目的： 从 Flex 格式的翼型极线文件读取数据。
   • 初始化 (__init__)： 解析 Flex 极线文件，将数据组织成 set 字典（Flex 文件通常只有一个设置），并为每个翼型创建 curve_interpolate 对象。
   • save_file() 方法： 将当前编辑的极线数据保存回 Flex 格式的新文件。

5. aero_damp_analyzer 类：

   • 目的： 核心类，用于分析和交互式编辑翼型极线，并计算气动阻尼。

   • 初始化 (__init__)：

     • 保存原始和可编辑的 CL/CD 曲线数据。
     • 初始化用于点选择和曲线编辑的内部状态变量（例如 ipoint_selected, dc1, dc2）。
     • 设置插值后的攻角 (aoas) 和方位角 (psis) 范围。
     • 初始化气动阻尼的四个分量 (W_tran1, W_tran2, W_tors1, W_tors2) 和总阻尼 eta。
     • 设置气动阻尼计算所需的参数：beta、gama (伽马)、phi (入流角) 和 k (减缩频率)。
     • 调用 cpf.compute_damping_terms 和 cpf.compute_damping_eta 来计算阻尼项和总阻尼。
     • 设置 Matplotlib 图形界面：创建两个子图，一个用于显示 CL/CD 曲线，另一个用于显示阻尼系数的等高线图。
     • 事件连接：连接鼠标点击 (onpick) 和键盘按键 (key_input) 事件，实现交互式编辑和参数调整。

   • 方法：

     • update_interpolated_values()：根据当前编辑的 CL/CD 曲线更新插值后的值及其导数。

     • onpick(event)：处理鼠标点击事件，用于选择或取消选择 CL/CD 曲线上的点。

     • key_input(event)：处理键盘输入事件，允许用户：

       • 通过上下箭头键微调选定点的 CL/CD 值。
       • 按 'u' 键更新阻尼等高线图。
       • 按 'b'/'B'、't'/'T'、'f'/'F'、'k'/'K' 键调整 beta、gama、phi 和 k 参数。
       • 按 's' 键保存当前图为 PNG。
       • 按 'S' 键保存编辑后的极线数据到文件。

     • update_parameters_in_title()：更新图的标题，显示当前的参数值。

     • update_plot()：重新绘制 CL/CD 曲线和选定点。

     • update_damping_plot()：重新绘制阻尼等高线图。

6. casedamp(...) 函数：

   • 目的： 包的入口函数，用于启动气动阻尼分析器。
   • 功能： 尝试读取 HAWC2 或 Flex 格式的极线文件，然后实例化 aero_damp_analyzer 类并显示交互式绘图界面。

计算稳态运行状态叶片模态、频率

从structural_blade_modes_H2_elements.py入口，该程序计算了DTU 10MW风电机组在standstill、风速为0，转速10rpm、风速11m/s，转速9.6rpm三种工况下叶片的十阶模态频率和形状，并于hawc2软件计算结果对比。

脚本中调用的主要函数及其功能：

1. 初始化和设置：

• os.getcwd(): 获取当前工作目录，用于构建输入和输出文件的路径。
• os.path.join(work_folder, 'DTU10MW_H2_elements_few_ops.json'): 将工作目录路径与文件名连接起来，创建 JSON 配置文件的完整路径。
• casestab.rotor_models(...): 通过从指定的 JSON 文件加载风力涡轮机模型配置来初始化 rotor_models 对象。此对象包含整个涡轮机的定义，包括其子结构（如叶片）。

2. 静止叶片结构模态（第一个运行工况）：

• rotor_models.models[0].substructures[1]: 访问第一个运行模型（索引 0），然后选择叶片子结构（假设索引 1 代表叶片）。
• rotor_models.models[0].update_all_substructures(0.0): 更新第一个模型所有子结构的质量和刚度矩阵。0.0 表示静止状态，没有旋转或外部力。
• blade.compute_modes(): 使用 CASEStab 方法计算叶片子结构的固有频率和模态形状。结果存储在 casestab_sol 中。
• np.loadtxt(...): 从 blade.cmb 加载频率数据，从 blade_standstill.amp 加载模态形状数据（HAWCStab2 结果）到 NumPy 数组中进行比较。skiprows 用于跳过文件中的标题信息。
• np.zeros(...): 创建全零的 NumPy 数组，用于预分配存储 HAWCStab2 模态形状的振幅 (modeamp) 和相位 (modepha) 的空间。
• np.radians(...): 将从 HAWCStab2 结果中读取的相位角从度转换为弧度。
• blade.plot_substructure_modes(subs_mode_solutions, nmodes, False, fname): 绘制计算出的结构模态。它接受一个解决方案列表（CASEStab 和 HAWCStab2）、要绘制的模态数量 (nmodes)、一个布尔标志（可能用于显示选项）以及用于保存绘图的文件名。

3. 额定转速下的叶片结构模态（第二个运行工况）：

• rotor_models.models[1].substructures[1]: 从第二个运行模型（索引 1）中选择叶片子结构。
• rotor_models.models[1].compute_substructure_steady_state_deformation(1): 计算第二个模型子结构的稳态变形。此步骤考虑了额定转速下作用在旋转叶片上的离心力。
• blade.compute_modes(): 在额定转速条件下计算叶片的结构模态。
• np.loadtxt(...), np.zeros(...), np.radians(...): 与静止部分类似，这些函数用于加载、处理和准备 HAWCStab2 结果 (blade_10_rpm.amp) 以进行比较。
• blade.plot_substructure_modes(...): 绘制额定转速工况下的模态，比较 CASEStab 和 HAWCStab2 的结果。

4. 11 米/秒风速下偏转叶片的结构模态（第三个运行工况）：

• rotor_models.models[2].substructures[1]: 从第三个运行模型（索引 2）中选择叶片子结构。
• rotor_models.models[2].compute_rotor_stationary_steady_state(0): 计算整个转子模型在 11 米/秒风速下的静止稳态。这涉及计算转子在气动载荷和重力载荷共同作用下的平衡位置，这会影响叶片的结构行为。
• blade.compute_modes(): 在 11 米/秒运行条件下计算叶片的结构模态。
• np.loadtxt(...), np.zeros(...), np.radians(...): 这些函数再次用于加载、处理和准备 HAWCStab2 结果 (blade_11ms.amp) 以进行比较。
• blade.plot_substructure_modes(...): 绘制 11 米/秒运行工况下的模态，比较 CASEStab 和 HAWCStab2 的结果。

update_all_substructures函数

这个函数位于 casetoolbox/casestab/model_assembler.py 中的 model 类中。它的主要目的是在给定时间 t 的情况下，更新模型中所有子结构（substructures）的当前位置、方向、内部力和惯性状态。这对于模拟风力涡轮机在不同运行条件下的动态行为至关重要。

功能概述：

该函数遍历模型中的每一个子结构，并执行以下关键更新：

1. 子结构内部更新 (s.subs.update_substructure())：

   • 这个调用会更新子结构（例如，corotbeam_substructure 实例）内部的节点三联体（nodal triads）和位置。对于柔性梁，这意味着根据其当前的变形状态更新每个梁单元的几何形状。
   • 它还会更新局部变形及其一阶和二阶导数，以及每个形状函数阶次的挠度子向量及其导数。这些是计算内部力和惯性属性的基础。

2. 轴承状态更新 (s.bearing.state.update(t))：

   • 如果子结构包含轴承（s.bearing.bear_flag 为 True），则会更新轴承的状态。这可能包括根据时间 t 更新轴承的旋转角度或速度，特别是对于“恒定速度”类型的轴承。

3. 子结构在地面固定坐标系中的位置和方向更新：

   • 函数会计算每个子结构在地面固定坐标系中的当前位置向量 s.r0 和方向矩阵 s.R0。
   • 如果子结构连接到另一个子结构（s.isubs_connection > -1），则其位置和方向是根据其支撑子结构的位置、方向以及连接节点处的局部变形来计算的。
   • 如果子结构是地面固定的（s.isubs_connection == -1），则其位置通常为零向量，方向由其初始方向矩阵 s.S 和轴承的旋转 s.bearing.state.B 决定。
   • 还会计算这些位置和方向的时间导数（s.R0TdR0, s.R0TddR0, s.R0Tdr0, s.R0Tddr0），这些对于惯性力的计算是必需的。

4. 弹性内部力和刚度矩阵更新 (s.subs.update_elastic_internal_forces_and_stiffness())：

   • 这个调用会计算子结构内部的弹性恢复力 (s.Fint) 和弹性刚度矩阵 (s.K)。对于柔性梁，这涉及到根据其当前的变形计算梁单元的应变能和相应的力。

5. 惯性力更新 (s.subs.update_inertia())：

   • 这个调用会更新子结构的惯性状态，包括总质量、惯性矩以及与质量和惯性相关的导数。

6. 离心力、陀螺矩阵和离心刚度矩阵更新 (s.subs.inertia.compute_local_centrifugal_forces_and_matrix(...))：

   • 根据子结构的时间导数（速度和角速度），计算离心力 (Fc1)、陀螺矩阵 (G11) 和离心刚度矩阵 (Kc11)。这些是旋转体动力学中的重要组成部分。

7. 气动弹性耦合更新 (s.subs.update_aeroelastic_coupling())：

   • 如果子结构是叶片的一部分（s.subs.aero_part 为 True），则会更新气动弹性耦合。这包括更新气动计算点（ACP）的位置和方向，以及计算将分布式气动力和力矩转换为节点力所需的转换矩阵。

update_all_substructures 是一个核心函数，它在每个时间步或每次迭代中被调用，以确保模型中所有子结构的几何、力学和惯性状态都得到准确更新，为后续的动力学或模态分析提供正确的基础。

compute_rotor_stationary_steady_state 函数

这个函数位于 casetoolbox/casestab/model_assembler.py 中的 model 类中。它的主要目的是计算转子在静止稳态下的变形和力平衡。这通常涉及到迭代求解气动载荷和结构变形之间的耦合问题，直到达到收敛。

功能概述：

该函数通过一个迭代过程来寻找转子的静止稳态，即在给定风速和操作条件下，气动力和结构内部力达到平衡时的变形状态。

参数：

• irotor: 要计算的转子的索引。
• Rnorm_limit: 残差的收敛容差。当力平衡残差的范数低于此值时，迭代停止。
• include_deform: 布尔值，指示是否在迭代中考虑结构变形。如果为 False，则只计算气动力，不更新结构变形。

工作流程和关键步骤：

1. 初始更新 (self.update_all_substructures(0.0))：

   • 在开始迭代之前，首先调用 update_all_substructures(0.0)。这会初始化所有子结构在时间 t=0.0 时的位置、方向、内部力和惯性状态。这是稳态计算的起点。

2. 指向目标转子和叶片：

   • r = self.rotors[irotor]: 获取要分析的特定转子对象。
   • b = r.blades[0]: 由于稳态计算假设转子是轴对称的（axissym），因此通常只使用第一个叶片的数据进行计算，然后将结果乘以叶片数量。

3. 外层迭代（气动力收敛循环）：

   • while Anorm > Rnorm_limit:: 这是一个外层循环，用于确保气动力的分布在迭代过程中收敛。Anorm 代表气动力变化的范数。

   • 重置功率和推力： 在每次外层迭代开始时，将转子的功率 (r.power) 和推力 (r.thrust) 重置为零。

   • 更新稳态叶片力 (self.update_steady_state_all_rotors())：

     • 这个调用会更新所有转子上的气动力。对于每个叶片，它会调用 rotor.update_steady_blade_forces()。
     • rotor.update_steady_blade_forces() 内部会使用 BEM（叶素动量理论）方法，通过 scipy.optimize.root 求解每个气动计算点（ACP）处的轴向和切向诱导因子，从而计算出气动力和力矩。这个过程本身可能包含一个内部的求解循环。

4. 内层迭代（结构变形收敛循环）：

   • while Rnorm > Rnorm_limit:: 这是一个内层循环，用于确保结构变形达到力平衡。Rnorm 代表力平衡残差的范数。

   • 更新所有子结构 (self.update_all_substructures(0.0))： 在每次内层迭代中，再次调用此函数，以根据当前的变形状态更新所有子结构的几何、力和惯性。

   • 获取气动力和力矩： 从叶片中获取当前的气动力 (force) 和力矩 (moment)。

   • 计算广义力 (faero) 和节点力 (fnode)：

     • faero=self.substructures[isubs].subs.TQf@force + self.substructures[isubs].subs.TQm@moment: 将分布式气动力和力矩转换为子结构节点上的广义力。
     • fnode=faero-self.substructures[isubs].subs.inertia.Fc1: 计算总节点力，包括气动力和离心力（Fc1）。

   • 构建切线刚度矩阵 (K)：

     • K = self.substructures[isubs].subs.K+self.substructures[isubs].subs.inertia.Kc11+Ksf: 构建总切线刚度矩阵。它包括弹性刚度 (subs.K)、离心刚度 (subs.inertia.Kc11) 和由气动力引起的几何刚度 (Ksf)。Ksf 是通过 subs.TKQf 和 subs.TKQm 矩阵计算的，这些矩阵表示气动力对刚度的贡献。

   • 计算残差 (R) 和位移增量 (x)：

     • R=fnode-self.substructures[isubs].subs.Fint: 计算力平衡残差，即外部节点力与内部弹性力 (Fint) 之间的差值。
     • x=solve(K,R): 求解线性方程组 K * x = R，得到当前的位移增量 x。

   • 更新变形 (self.substructures[isubs].subs.q+=relaxation_factor*x)：

     • 将计算出的位移增量 x 应用到子结构的广义坐标 q 上，从而更新其变形。这里使用了 relaxation_factor 来控制步长，以帮助收敛，防止过大的增量导致不稳定。

   • 检查内层收敛： Rnorm=np.linalg.norm(R) 计算残差的范数，如果低于 Rnorm_limit，则内层循环收敛。

5. 更新功率和推力：

   • 在内层循环收敛后，根据当前的力分布计算转子的总功率 (r.power) 和推力 (r.thrust)。这涉及到将子结构坐标系中的力转换到地面固定坐标系，并与速度和角速度相乘。

6. 检查外层收敛：

   • Anorm=np.linalg.norm(b.f-force0): 计算当前气动力分布与上一次外层迭代开始时的气动力分布之间的差异范数。如果 include_deform 为 True，则此值用于判断外层循环是否收敛。
   • 如果 include_deform 为 False，则 Anorm 被设置为 0，表示只进行一次气动力计算，不考虑变形对气动力的影响。

7. 轴对称转子处理：

   • 如果转子是轴对称的（'axissym' in r.type），则最终的功率和推力会乘以叶片数量 r.Nb，因为计算是基于单个叶片的。
   • r.CP = r.power/Pkin: 计算功率系数 CP，这是风力涡轮机性能的关键指标。

compute_modes 函数

这个函数位于 casetoolbox/casestab/model_assembler.py 中的 substructure 类中。它的主要目的是计算子结构的结构模态（固有频率和模态形状）。

功能概述：

该函数主要针对柔性梁子结构（corotbeam 类型）执行模态分析：

1. 类型检查 (if subs.type == 'corotbeam':)：

   • 首先检查子结构类型是否为 corotbeam。如果不是，则返回空的模态解决方案，因为只有柔性梁才具有结构模态。

2. 获取梁单元信息：

   • zpos=subs.pos[2,1:]: 获取梁单元的纵向位置（z坐标），用于后续绘图。

3. 设置质量和刚度矩阵：

   • K = subs.K+subs.inertia.Kc11: 构建总刚度矩阵 K。它包括梁的弹性刚度矩阵 subs.K 和由离心力引起的几何刚度（或离心刚度）矩阵 subs.inertia.Kc11。
   • M = subs.inertia.M11: 获取梁的质量矩阵 M。

4. 求解特征值问题：

   • Minv = np.linalg.inv(M): 计算质量矩阵的逆。
   • A = np.concatenate(...): 构造一个状态空间矩阵 A。这是将二阶动力学方程 M * ddot(q) + G * dot(q) + K * q = F 转换为一阶状态空间形式 dot(x) = A * x 的标准方法，其中 x = [q; dot(q)]。这里，它考虑了质量矩阵 M、刚度矩阵 K 和陀螺矩阵 subs.inertia.G11（尽管在模态分析中通常假设无阻尼，但这里可能包含了陀螺效应）。
   • evals,evecs=np.linalg.eig(A): 求解广义特征值问题。evals 包含特征值（复数，其虚部与固有频率相关），evecs 包含对应的特征向量（模态形状）。

5. 排序和处理特征值/特征向量：

   • isort=np.argsort(np.imag(evals))[6*subs.nelem:]: 根据特征值的虚部（即固有频率）对特征值和特征向量进行排序。6*subs.nelem 可能是为了跳过一些零频率或不相关的模态（例如，刚体模态）。
   • freq=np.imag(evals[isort])/np.pi/2.0: 从排序后的特征值的虚部计算固有频率（Hz）。
   • modeamp 和 modepha：初始化用于存储模态振幅和相位的 NumPy 数组。
   • 循环遍历每个模态，从特征向量中提取位移和旋转的振幅和相位。
   • np.abs(...) 和 np.angle(...): 用于从复数特征向量中提取振幅和相位。
   • np.radians(...): 将旋转角度转换为弧度。
   • Rodrigues 参数到欧拉角的转换： 对于旋转自由度，代码将 Rodrigues 参数转换为欧拉角。这里特别提到“evecs 的尺度（范数）必须很小”，这暗示了 Rodrigues 参数在处理大旋转时可能遇到的非线性问题，但在小变形模态分析中通常是有效的。
   • 模态命名和归一化： 根据模态在三个方向（Inplane, Out-of-plane, Torsion）上的最大振幅，为模态命名（例如“1F”、“1E”、“1T”）。然后对模态形状进行归一化，通常是使最大振幅为 1。

6. 保存模态解决方案：

   • 将计算出的模态信息（z位置、频率、振幅、相位、名称）存储在一个字典 subs_mode_solution 中，并返回。

compute_modes 函数是进行结构动力学分析的核心，它通过求解特征值问题来确定柔性子结构的固有振动特性，这些特性对于理解结构响应和避免共振至关重要。