2025-06-27 16:08:58 +08:00
# 状态空间方程基础
## 状态空间描述
状态空间描述（内部描述）：通过建立系统内部状态和系统的输入以及输出之间的数学关系，来描述系统的行为。
系统输入引起系统状态的变化，输入及系统状态共同导致系统输出的变化
状态空间描述是对系统完全的描述，能表征系统的一切动力学
## 状态
系统的状态：能够完全表征系统时间域行为的最小内部变量组
$$
x(t) = {[x_1(t), x_2(t), ..., x_n(t)]}^T
$$
状态变量：构成系统状态的每一个变量。状态变量构成的列向量为状态向量
## 状态空间
状态向量的取值空间称为状态空间
状态空间：
系统在任意时刻的状态在状态空间中都可以用一点表示。
状态轨线：随着时间的推移，系统状态$x(t)$在状态空间所留下的轨迹称为状态轨线或状态轨迹
## 状态空间描述的形式
1、状态方程：
**表征系统中输入所引起的内部状态变化的过程。**
连续时间系统：
$$
\dot{x}(t) = f[x(t), u(t), t]
$$
离散时间系统
$$
{x}(t_{k+1}) = f[x(t_k), u(t_k), t_k]
$$
2、输出方程：
连续时间系统：
$$
{y}(t) = g[x(t), u(t), t]
$$
离散时间系统
$$
{y}(t_{k}) = g[x(t_k), u(t_k), t_k]
$$
3、状态空间表达式： ，
对于一个n个状态变量， ， ：
$$
\begin{array}{r}
{\dot{x}(t) = f[x(t), u(t), t] }\\
{{y}(t) = g[x(t), u(t), t]}
\end{array}
$$
连续线性系统：
$$
\begin{array}{r}
{\dot{x}(t) = A(t)x(t) + B(t)u(t) }\\
{{y}(t) = C(t)x(t)+ D(t)u(t)}
\end{array}
$$
其中：
- $A(t)$系统矩阵
- $B(t)$控制矩阵
- $C(t)$观测矩阵
- $D(t)$前馈矩阵
2025-06-25 16:42:20 +08:00
# 1. 线性化分析背景
线性分析计算将气弹模型简化为用户请求计算的**每个工作点处的线性系统**。**以状态空间形式表示的线性系统方程组为**：
$$
\begin{array}{r}{\underline{{\dot{\mathbf{x}}}}=\mathbf{A}\underline{{\mathbf{x}}}+\mathbf{B}\underline{{\mathbf{u}}}}\\ {\underline{{\mathbf{y}}}=\mathbf{C}\underline{{\mathbf{x}}}+\mathbf{D}\underline{{\mathbf{u}}}}\end{array}
$$
$$
\underline{{\mathbf{x}}}=\mathbf{x}-\mathbf{x_{0}},\quad\underline{{\mathbf{y}}}=\mathbf{y}-\mathbf{y_{0}},\quad\mathrm{~and~}\quad\underline{{\mathbf{u}}}=\mathbf{u}-\mathbf{u_{0}}
$$
其中 $\mathbf{x}$ 是状态向量，代表系统状态；$\mathbf{u}$ 是系统输入向量；$\mathbf{y}$ 是系统输出向量。归一化向量 $\underline{{\mathbf{x}}}、\underline{{\mathbf{y}}}$ 和 $\underline{{\mathbf{u}}}$ 代表偏离平衡态的量。
矩阵 $\mathbf{A}$、$\mathbf{B}$、$\mathbf{C}$ 和 $\mathbf{D}$ 代表这些向量之间的线性化关系。这代表对完整模型的一个简化，该模型使用一组完全非线性方程来计算状态导数和输出。
$$
\begin{array}{r}{\dot{\mathbf{x}}=f(t,\mathbf{x},\mathbf{u})}\\ {\mathbf{y}=h(t,\mathbf{x},\mathbf{u}).}\end{array}
$$
需要注意的是，为了能够建立适当的线性化风轮动力学系统，需要考虑以下模型准备原则：
- 需移除方位角依赖性，包括风切变、偏航、风轮不平衡等。
- 需移除无法线性化的物理效应，例如风湍流、黏滞滑动等。
为了进行线性分析， ， （ ） ，
对于每个输入或状态， ， ，
然后， ， ， ，
图 1：
# 2. 多叶片坐标系转换
为了线性化计算或坎贝尔图的绘制，建议选择多叶片坐标变换，该变换会生成与非旋转坐标系耦合的模态，包括风轮的后向和前向挥舞模态。这基于 (Bir, 2008) 和 (Hansen, 2003) 中提出的理论。线性化后的模型对方位角高度依赖， ， ， ，
## 2.1. 单个模态转换
一个三叶片系统，其方位角为 $\psi_{1}$ 到 $\psi_{3}$，从非旋转坐标系到旋转坐标系的位移变换矩阵如下：
$$
\begin{align}
\begin{bmatrix}
q_{1} \\
q_{2} \\
q_{3} \\
\end{bmatrix} & = \widetilde{{t}}_{NR \rightarrow R}
\begin{bmatrix}
q_{0} \\
q_{c} \\
q_{s} \\
\end{bmatrix},
\end{align}
$$
$$
\tilde{\mathbf{t}}_{N R\to R}=\left[\begin{array}{l l l}{1}& {\cos\psi_{1}}& {\sin\psi_{1}}\\ {1}& {\cos\psi_{2}}& {\sin\psi_{2}}\\ {1}& {\cos\psi_{3}}& {\sin\psi_{3}}\end{array}\right].
$$
多叶片坐标变换通常以上述形式引用， ， ， ， ：
$$
\begin{align}
{t}_{R \rightarrow NR} = \frac{1}{3}\begin{bmatrix}
1 & 1 & 1 \\
2\cos\psi_{1} & 2\cos\psi_{2} & 2\cos\psi_{3} \\
2\sin\psi_{1} & 2\sin\psi_{2} & 2\sin\psi_{3} \\
\end{bmatrix}.
\end{align}
$$
对于具有任意数量叶片 $(n)$ 的风轮的通用变换矩阵
$$
\mathbf{t}_{R\rightarrow N R}=\frac{1}{n}\left[\begin{array}{c c c c c}{1}& {1}& {1}& {\cdots}& {1}\\ {2\cos\psi_{1}}& {2\cos\psi_{2}}& {2\cos\psi_{3}}& {\cdots}& {2\cos\psi_{n}}\\ {2\sin\psi_{1}}& {2\sin\psi_{2}}& {2\sin\psi_{3}}& {\cdots}& {2\sin\psi_{n}}\\ {2\cos j\psi_{1}}& {2\cos j\psi_{2}}& {2\cos j\psi_{3}}& {\cdots}& {2\cos j\psi_{n}}\\ {2\sin j\psi_{1}}& {2\sin j\psi_{2}}& {2\sin j\psi_{3}}& {\cdots}& {2\sin j\psi_{n}}\\ {\vdots}& {\vdots}& {\vdots}& {\ddots}& {\vdots}\\ {1}& {-1}& {1}& {\cdots}& {(-1)^{n}}\end{array}\right],
$$
其中最后一行是变换到微分模态，且仅在叶片数量为偶数时存在。对于叶片数量为奇数的风轮，最后一行将是一个正弦循环行。计数器 $j$ 从 1 到 $(n-1)/2$ (当 $n$ 为奇数时)，从 2 到 $(n-2)/2$ (当 ${n}$ 为偶数时)。
省略矩阵表示，非旋转坐标可以计算如下：
$$
\begin{array}{l}{\displaystyle q_{0}=\frac{1}{n}\sum_{i=1}^{n}q_{i}}\\ {\displaystyle q_{c j}=\frac{2}{n}\sum_{i=1}^{n}q_{i}\cos{\left(j\psi_{i}\right)}}\\ {\displaystyle q_{s j}=\frac{2}{n}\sum_{i=1}^{n}q_{i}\sin{\left(j\psi_{i}\right)}}\\ {\displaystyle q_{d}=\frac{1}{n}\sum_{i=1}^{n}q_{i}(-1)^{n}}\end{array}
$$
以三叶风轮为例，现在计算导数变换矩阵。每个方位角 $\psi_{i}$ 可以表示为（假设恒定的）风轮转速 $\Omega$ 和初始方位角 $\Psi_{i}$ 的线性关系。
$$
\psi_{i}=\Omega t+\Psi_{i}.
$$
对变换矩阵求时间导数，得到：
$$
\begin{align*}
\dot{{t}}_{R \rightarrow NR} =
\frac{\Omega}{3}
\begin{bmatrix}
0 & 0 & 0 \\
-2\sin\psi_{1} & - 2\sin\psi_{2} & - 2\sin\psi_{3} \\
2\cos\psi_{1} & 2\cos\psi_{2} & 2\cos\psi_{3} \\
\end{bmatrix}
\end{align*}
$$
$$
\begin{align*}
\ddot{{t}}_{R \rightarrow NR} =
-\frac{\Omega^{2}}{3}
\begin{bmatrix}
0 & 0 & 0 \\
2\cos\psi_{1} & 2\cos\psi_{2} & 2\cos\psi_{3} \\
2\sin\psi_{1} & 2\sin\psi_{2} & 2\sin\psi_{3} \\
\end{bmatrix}
\end{align*}
$$
## 2.2. 系统变换矩阵
需要一个包含位移和速度状态的整个状态列表的变换矩阵。正如我们在公式(1)中已经确立的，对于位移状态，
$$
{\bf q}_{N R}={\bf t}_{R\rightarrow N R}{\bf q}_{R}
$$
对公式(1)求时间导数，得到：
$$
\dot{\bf q}_{N R}={\bf t}_{R\rightarrow N R}\dot{\bf q}_{R}+\dot{\bf t}_{R\rightarrow N R}{\bf q}_{R}
$$
将 ${\bf q}_{N R}$ 和 $\dot{\bf q}_{N R}$ 组合成一个包含所有状态（位移和速度）的向量，允许我们定义一个与 $A$ 具有相同维度的公共变换矩阵 $\mathbf{T}$。 我们定义
$$
\begin{align}
{T} & := \begin{bmatrix}
{t}_{R \rightarrow NR} & 0 \\
\dot{{t}}_{R \rightarrow NR} & {t}_{R \rightarrow NR} \\
\end{bmatrix}
\end{align}
$$
需要注意的是，一般来说，位移和速度状态并非以这种方式排列，系统变换矩阵 $\mathbf{T}$ 会发生置换。系统变换矩阵非奇异，可以计算其逆矩阵。
系统变换矩阵的导数可以很容易地推导出来。
$$
\begin{align}
\dot{{T}}=
\begin{bmatrix}
\dot{{t}}_{R \rightarrow NR} & 0 \\
\ddot{{t}}_{R \rightarrow NR} & \dot{{t}}_{R \rightarrow NR} \\
\end{bmatrix}.
\end{align}
$$
## 2.3. 变换A, B, C, D矩阵
我们考虑在旋转参考系下的线性模型方程，并定义：
$$
\begin{align}
{x}_{R} := \begin{bmatrix}
{q}_{R} \\ \dot{{q}}_{R}
\end{bmatrix}
\end{align}
$$
为了表达主要系统为
$$
\begin{array}{r}{\dot{\mathbf{x}}_{R}=\mathbf{A}_{R}\mathbf{x}_{R}+\mathbf{B}_{R}\mathbf{u}}\\ {\mathbf{y}=\mathbf{C}_{R}\mathbf{x}_{R}+\mathbf{D}_{R}\mathbf{u}}\end{array}
$$
相对于旋转叶片坐标系。
状态向量从旋转坐标系变换到非旋转坐标系的表达式如下：
$$
\mathbf{x}_{N R}=\mathbf{T}\mathbf{x}_{R}
$$
导数为：
$$
\dot{{\bf x}}_{N R}={\bf T}\dot{\bf x}_{R}+\dot{\bf T}{\bf x}_{R}.
$$
结合方程17和15，
$$
\begin{array}{r l}& {\dot{\mathbf{x}}_{N R}=\mathbf{T}\left(\mathbf{A}_{R}\mathbf{x}_{R}+\mathbf{B}_{R}\mathbf{u}\right)+\dot{\mathbf{T}}\mathbf{x}_{R}}\\ & {\qquad=\Big(\mathbf{T}\mathbf{A}_{R}+\dot{\mathbf{T}}\Big)\mathbf{x}_{R}+\mathbf{T}\mathbf{B}_{R}\mathbf{u}}\\ & {\qquad=\Big(\mathbf{T}\mathbf{A}_{R}+\dot{\mathbf{T}}\Big)\mathbf{T}^{-1}\mathbf{x}_{N R}+\mathbf{T}\mathbf{B}_{R}\mathbf{u}}\end{array}
$$
进而得到
$$
\begin{array}{r l}& {\mathbf{A}_{N R}=\Big(\mathbf{TA}_{R}+\dot{\mathbf{T}}\Big)\mathbf{T}^{-1}}\\ & {\mathbf{B}_{N R}=\mathbf{TB}_{R}}\end{array}
$$
类似的，
$$
\begin{array}{r}{\mathbf{y}=\mathbf{C}_{R}\mathbf{x}_{R}+\mathbf{D}_{R}\mathbf{u}}\\ \\{=\mathbf{C}_{R}\mathbf{T}^{-1}\mathbf{x}_{N R}+\mathbf{D}_{R}\mathbf{u}}\end{array}
$$
对于线性化系统的输出，我们定义
$$
\mathbf{C}_{N R}:=\mathbf{C}_{R}\mathbf{T}^{-1},\quad\mathrm{~and~}\mathbf{D}_{N R}:=\mathbf{D}_{R}
$$
对于与线性模型输出相关的矩阵，本部分内容完成了对非旋转参考系下的线性模型推导。
旋转变换仅应用于状态，这些状态代表了定义在所有叶片上的数学模型中的自由度。这些状态包括叶片模态和动态失速状态，而其他任何单个叶片状态，例如桨距角位置、速率、执行器内部状态等，均不进行变换。矩阵 T 在未进行变换的状态和状态导数对应的行和列上具有单位对角线元素。对于其他行和列，
# 3. 耦合模态
## 3.1. 计算耦合模态
坎贝尔图和叶片稳定性分析是对在每个指定工作点处的**矩阵 A** 的分析。**每个耦合模态对应一个特征值及其特征向量**。**给定矩阵 A 的一个（复数）特征值 λ**，
Figure 1: Argand Diagram阿根图
每个耦合模态的未耦合模态贡献由其特征向量决定。如果耦合模态具有二阶状态（结构状态）的贡献，这些状态由状态向量中的两个状态表示，则使用位移状态来确定贡献。
以原始形式，这些特征向量贡献代表每个模态的相对位移，可用于构建耦合模态形状。然而，在坎贝尔图中，这些贡献已被归一化。这是通过修改特征向量矩阵来实现的，使得每一行和每一列都具有单位和。这会增加具有高质量和刚度的模态的百分比贡献，这些模态在位移中贡献很小，但在能量方面贡献很大。
每个贡献的相位，$\phi_{i\,,}$ 由相应复特征向量元素的参数决定，即 $v_{i},$ 即。
$$
\phi_{i}=\arctan\left(\frac{\operatorname{Im}\left(v_{i}\right)}{\operatorname{Re}\left(v_{i}\right)}\right).
$$
## 3.2. 命名耦合模态
在计算耦合模态，例如在坎贝尔图分析中，以下章节将详细介绍命名规范。重点关注多叶片坐标变换使用时的行为。
### 3.2.1. 支撑结构模态
对于支撑结构模态，耦合模态以对贡献最大的整个塔架模态命名。整个塔架模态通过后续带有固定-自由边界条件的特征分析来计算， （ ） ，
### 3.2.2. 旋转坐标系风轮模态
如果未使用多叶片坐标变换(MBC)来描述风轮模态，则以下逻辑适用于命名耦合风轮模态：
- 如果单个叶片模态贡献超过75%，则将耦合风轮模态命名为该叶片模态。换句话说，该模态被称为“叶片”模态，而不是“风轮”模态。
- ， ， ，
### 3.2.3. 非旋转坐标系风轮模态
**如果应用了 MBC 变换，则各个叶片模态将被转换成一组风轮模态**。对于三叶风轮，
在转换和重命名各个叶片模态之后，命名耦合风轮模态。旋转模态的识别遵循下表中的逻辑。
| Coupled mode name | 1st uncoupled mode | 2nd uncoupled mode | Phase angle ($\phi_{2}-\phi_{1})$) |
| ----------------- | ------------------ | ------------------ | ---------------------------------- |
| Forward whirl
| < 0.0 |
| Backward whirl < 0.0 |
|
如果一个耦合模态不符合颤动模态的标准，则该模态会被命名为其主要贡献者的名称。这与风轮在旋转参考系中的模态以及支撑结构模态的命名逻辑类似。
# 4. 连接不同工况下的耦合模态
**坎贝尔图显示了不同耦合位移模态相对于风轮转速的频率**， （ ） ，
## 4.1. 连接过程
给定一组在每个工作点耦合的模态，**创建最终的坎贝尔图的关键步骤是识别相邻工作点中相似的模态**， ，
**通过比较在扩展模态保证准则（ ） （ ） ， ， （ & Shapley, 1962） ，
为了确保最终坎贝尔图中的耦合模态序列主要涉及结构动力学，在第一个工作点执行一次仅包含结构状态的耦合模态计算。然后，将这些仅结构模态与如上所述的第一个工作点的耦合模态连接起来，这有效地排除了主要来自气动状态贡献的耦合模态。
![[images/Pasted image 20250625151451.png]]
图 1： ，
接下来的连接过程如下：
1. 计算在第一个工作点的一组仅结构模态。这些模态也将构成模态序列的基础，代表最终的坎贝尔图中的线。
2. 将仅结构模态与在第一个工作点的耦合模态连接起来。需要注意的是，耦合模态的数量通常大于模态序列的数量，这意味着并非所有耦合模态都包含在模态序列中。
3. 将在当前点（第一个工作点）包含在模态序列中的耦合模态与下一个点（第二个工作点）的耦合模态连接起来。重复此过程直至到达最后一个工作点。
### 4.1.1. 耦合模态序列的命名与排除
一种耦合模态级数，仅根据第一个工作点处的结构-模态耦合贡献来命名（关于耦合模态命名的更多细节请参见此处）。耦合模态的贡献，因此其形状，在工作点范围的不同情况下可能会发生显著变化，因此仅凭名称无法确定模态的特性。
一种耦合模态级数，其中包含所有工作点处具有**实特征值**（因此没有振荡行为）的耦合模态，将被排除在最终的坎贝尔图之外。这是因为此类模态无法引起共振行为，**而坎贝尔图的主要目的是检测共振行为。**
如果所有工作点处的耦合模态频率超过用户定义的上限值，则该耦合模态级数也会被排除。
# 5. 与风轮中心轴对准风向
图1所示， ， ， ， （ ） ， ， ， ， （ ） ， ， ， （ ） ， ，
图1： ， ： ， ： ， ：
为了避免倾斜轴向导致的载荷不平衡， ， ， ， ，
图2： ， ： ， ：
图3： ， ： ， ：
初始条件程序中确定风轮中心轴， ， ， ， ， ， ， ，
