gyz/obsidian_backup

Fork 0

yize 9a824b2201 vault backup: 2025-06-27 16:08:57

2025-06-27 16:08:58 +08:00

22 KiB

Raw Blame History

状态空间方程基础

状态空间描述

状态空间描述（内部描述）：通过建立系统内部状态和系统的输入以及输出之间的数学关系，来描述系统的行为。

系统输入引起系统状态的变化，输入及系统状态共同导致系统输出的变化

状态空间描述是对系统完全的描述，能表征系统的一切动力学

状态

系统的状态：能够完全表征系统时间域行为的最小内部变量组


x(t) = {[x_1(t), x_2(t), ..., x_n(t)]}^T

状态变量：构成系统状态的每一个变量。状态变量构成的列向量为状态向量

状态空间

状态向量的取值空间称为状态空间

状态空间：以n个线性无关的状态向量作为基底所组成的n维空间称为状态空间$R^n$。系统在任意时刻的状态在状态空间中都可以用一点表示。状态轨线：随着时间的推移，系统状态$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个状态变量，p个输入变量，q个输出变量的连续时间系统：


\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)$前馈矩阵

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}

需要注意的是，为了能够建立适当的线性化风轮动力学系统，需要考虑以下模型准备原则：

需移除方位角依赖性，包括风切变、偏航、风轮不平衡等。
需移除无法线性化的物理效应，例如风湍流、黏滞滑动等。

为了进行线性分析，Bladed依次取每个工作点，并找到风电机组的稳态条件（如同时域运行中的初始条件）。这意味着风轮没有加速，模态变形使得弹性载荷平衡外部载荷。这定义了$\mathbf{x_{0}},,\mathbf{y_{0}}$和$\mathbf{u_{0}}$的值，即所有扰动围绕的主要平衡点。

对于每个输入或状态，Bladed然后在平衡点两侧进行一系列幅度逐渐增加的扰动。人为地增加或减少状态或输入的数值，用这些修改后的数值求解系统，并记录状态导数和输出。扰动的数量和最大扰动幅度可以由用户定义。

然后，可以通过对状态导数与所有扰动值及其平衡值之间的关系进行线性回归，来推导矩阵A、B、C和D的元素。线性回归的梯度给出元素的值。如果相关系数小于最小相关系数，则认为该关系无效，并给该元素赋予零值。

图 1：示例线性回归计算单元 ${\bf A}_{7,4}$，其值为 -1.315，相关系数为 0.9982。平衡点以绿色显示。

2. 多叶片坐标系转换

为了线性化计算或坎贝尔图的绘制，建议选择多叶片坐标变换，该变换会生成与非旋转坐标系耦合的模态，包括风轮的后向和前向挥舞模态。这基于 (Bir, 2008) 和 (Hansen, 2003) 中提出的理论。线性化后的模型对方位角高度依赖，但当变换到非旋转坐标系后，结构动力学模型应仅表现出微弱的方位角依赖性。然而，对于2叶片风轮，仍然存在较强的方位角依赖性。

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].

多叶片坐标变换通常以上述形式引用，但主要目的是反向变换，从旋转坐标系变换到非旋转坐标系。一个3叶片系统的位移从旋转坐标系变换到非旋转坐标系的变换矩阵，是上述矩阵的逆矩阵，由以下给出：



\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}

类似的，方程16的转换得到



\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 在未进行变换的状态和状态导数对应的行和列上具有单位对角线元素。对于其他行和列，T 的元素代表上述为每组模态定义的基变换。需要注意的是，连接状态和状态导数的元素也需要通过对基变换方程进行微分来定义，同时考虑到方位角的变化率等于风轮转速（出于本目的，假设其为常数）。模型输入和输出不进行变换。

3. 耦合模态

3.1. 计算耦合模态

坎贝尔图和叶片稳定性分析是对在每个指定工作点处的矩阵 A 的分析。每个耦合模态对应一个特征值及其特征向量。给定矩阵 A 的一个（复数）特征值 λ，Bladed 会根据图 1 报告无阻尼频率undamped frequency (ωn)、阻尼频率damped frequency (ωd) 和阻尼比 (ζ)。

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. 支撑结构模态

对于支撑结构模态，耦合模态以对贡献最大的整个塔架模态命名。整个塔架模态通过后续带有固定-自由边界条件的特征分析来计算，用于线性化计算。该分析考虑了RNA（叶片）以及位于远端节点上的其他质量的影响。如果多个耦合支撑结构模态共享相同的整塔模态作为其主要贡献者，则通过附加字母A、B、C等来使耦合模态名称唯一化。

3.2.2. 旋转坐标系风轮模态

如果未使用多叶片坐标变换(MBC)来描述风轮模态，则以下逻辑适用于命名耦合风轮模态：

如果单个叶片模态贡献超过75%，则将耦合风轮模态命名为该叶片模态。换句话说，该模态被称为“叶片”模态，而不是“风轮”模态。
否则，耦合风轮模态以其主要贡献者命名，并通过附加字母A、B、C等来使其独一无二，以防多个耦合风轮模态共享相同的未耦合叶片模态作为主要贡献者。

3.2.3. 非旋转坐标系风轮模态

如果应用了 MBC 变换，则各个叶片模态将被转换成一组风轮模态。对于三叶风轮，通常存在集变模态collective、余弦循环模态cosine-cyclic和正弦循环sine-cyclic模态。所有叶片的 1 阶挥舞模态将被重命名为风轮 1 阶挥舞集变模态、风轮 1 阶挥舞正弦循环模态和风轮 1 阶挥舞余弦循环模态。如果叶片数量为偶数，则还会存在摆振模态。

在转换和重命名各个叶片模态之后，命名耦合风轮模态。旋转模态的识别遵循下表中的逻辑。

| Coupled mode name | 1st uncoupled mode | 2nd uncoupled mode | Phase angle (\phi_{2}-\phi_{1})) |

| ----------------- | ------------------ | ------------------ | ---------------------------------- |

如果一个耦合模态不符合颤动模态的标准，则该模态会被命名为其主要贡献者的名称。这与风轮在旋转参考系中的模态以及支撑结构模态的命名逻辑类似。

4. 连接不同工况下的耦合模态

坎贝尔图显示了不同耦合位移模态相对于风轮转速的频率，以及以风轮频率（P）的倍数给出的最重要激励频率。除了耦合模态的频率，坎贝尔图还显示了相应的阻尼比，其中包括结构阻尼和气动阻尼的影响。这两个特征的计算方法如文章“Calculating Coupled Modes”中所述，可用于识别需要进一步分析的关键运行点。

4.1. 连接过程

给定一组在每个工作点耦合的模态，创建最终的坎贝尔图的关键步骤是识别相邻工作点中相似的模态，这使得能够用线段连接这些相似的模态，从而给用户留下频率随风轮转速（或风速）连续变化的印象。这种连接过程通常具有挑战性，因为耦合模态在工作点之间会演变和改变其贡献。

通过比较在扩展模态保证准则（MACX）（Vacher, Jacquier, 和 Buchales, 2010）下的复特征向量和频率，并采用频率加权，可以识别出两个相邻工作点中的相似耦合模态。更具体地说，会为两个工作点中所有耦合模态组合计算频率加权 MACX 数值，形成一个评分矩阵，然后使用该矩阵通过 GaleShapley 算法（Gale & Shapley, 1962）连接这些模态。所有工作点中一系列相似的模态构成一个耦合模态序列，它代表了最终坎贝尔图中的一条线。

为了确保最终坎贝尔图中的耦合模态序列主要涉及结构动力学，在第一个工作点执行一次仅包含结构状态的耦合模态计算。然后，将这些仅结构模态与如上所述的第一个工作点的耦合模态连接起来，这有效地排除了主要来自气动状态贡献的耦合模态。

图 1：连接过程示意图，展示了耦合模态在工作点处的连接方式。点代表耦合模态。彩色线条代表耦合模态级数。

接下来的连接过程如下：

计算在第一个工作点的一组仅结构模态。这些模态也将构成模态序列的基础，代表最终的坎贝尔图中的线。
将仅结构模态与在第一个工作点的耦合模态连接起来。需要注意的是，耦合模态的数量通常大于模态序列的数量，这意味着并非所有耦合模态都包含在模态序列中。
将在当前点（第一个工作点）包含在模态序列中的耦合模态与下一个点（第二个工作点）的耦合模态连接起来。重复此过程直至到达最后一个工作点。

4.1.1. 耦合模态序列的命名与排除

一种耦合模态级数，仅根据第一个工作点处的结构-模态耦合贡献来命名（关于耦合模态命名的更多细节请参见此处）。耦合模态的贡献，因此其形状，在工作点范围的不同情况下可能会发生显著变化，因此仅凭名称无法确定模态的特性。

一种耦合模态级数，其中包含所有工作点处具有实特征值（因此没有振荡行为）的耦合模态，将被排除在最终的坎贝尔图之外。这是因为此类模态无法引起共振行为，而坎贝尔图的主要目的是检测共振行为。

如果所有工作点处的耦合模态频率超过用户定义的上限值，则该耦合模态级数也会被排除。

5. 与风轮中心轴对准风向

图1所示，即使遵循线性分析背景中的考虑因素，在某些情况下仍然需要调用方位角依赖性。当计算中考虑倾角，且塔架柔度强烈影响风轮姿态时，这种情况会发生。图1（左）提供了一个方位角独立的系统示意图，表明无论风轮的方位角角度如何，解都是相同的。然而，当考虑倾角时，如图1（中间）所示，风轮上部和下部受到的载荷会发生变化。这会产生类似于在风轮平面上添加“虚拟风切变”的载荷不平衡。当塔架柔度强烈影响风轮姿态时，情况会更加糟糕，如图1（右）所示。在这里，可以看到轮毂姿态可能会倾斜得更大，从而在风轮上产生更强的载荷不平衡。

图1：线性化计算中，风轮轴向相对于迎风方向的示意图。左：无倾角刚性风轮，中：考虑倾角的刚性风轮，右：考虑倾角的柔性风轮。

为了避免倾斜轴向导致的载荷不平衡，Bladed 引入了一种校正措施，以使风场与倾斜轴向对齐。该机制在图2中清晰地说明。可以观察到，当风场与倾斜轴向对齐时，风轮受到的载荷将与偏航角无关。这有效地消除了由于倾斜轴向引起的系统偏航角依赖性。这种校正也可以应用于浮式风力发电机，如图3所示。

图2：针对陆上风力发电机或固定基础海上风力发电机，将风场对准倾斜的轴向。左：方位角相关的系统，右：方位角无关的系统。

图3：针对浮动式风力发电机，将风场对准倾斜的轴向。左：方位角相关的系统，右：方位角无关的系统。

初始条件程序中确定风轮中心轴，此时风轮未加速，模态变形使得弹性载荷与外部载荷平衡。在初始条件程序中，进行迭代计算，并在每次迭代中调整风场，使其与倾斜的风轮中心轴平行。一旦找到合适的条件，风场方向就被固定在初始条件程序中最后确定的倾斜角度。随后，系统使用固定的风场方向进行扰动，用于计算Bladed中可用的任何线性化类型计算。此过程针对每个模拟的工作点进行。

22 KiB Raw Blame History Unescape Escape