一种计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法

未命名 08-14 阅读:114 评论:0


1.本发明属于电力系统建模技术领域,涉及一种计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法。


背景技术:

2.近年来,为了建设新型电力系统,以风电、光伏为代表的新能源在我国获得了迅猛发展。风电一般通过变流器并网,采用锁相环实现与电网的同步,由控制算法决定其输出功率,不再具有传统同步机组的惯量响应特性。电力电子器件具有控制灵活、快捷的优点,但风电机组也存在一些不足,其主要表现在:1)抗扰动能力较弱。对频率、电压偏差的耐受能力差,在频率或电压大幅扰动下更易脱网,对电力系统稳定不利。2)具有强非线性和离散性。为了防止变流器过载,风电机组设置了限幅、饱和环节。为保证故障期间的不脱网运行,设计了高/低压穿越控制策略,故风电机组表现出比传统同步机组更为复杂的非线性特性。
3.为了分析新型电力系统的动态特性,需要建立风电场动态模型。由于电力电子器件响应时间在毫秒级,因此,机电暂态模型已无法满足分析精度的要求,需要考虑其电磁暂态特性。由于电磁暂态仿真步长一般为几十微秒,且风电机组模型阶数很高,因此,对于大型风电场,直接采用详细模型进行电磁暂态仿真计算,其计算量、计算速度及数据存储量均难以承受。为此,需要对风电场进行电磁暂态等值。现有的风电场电磁暂态等值方法主要有两类:1)根据短路电流波形的相似性进行分群聚合;2)基于故障前风电机组转子转速进行分群聚合。但上述方法未计及风电机组低电压穿越特性,也未考虑风电场内部电压的差异,因此难以准确地反映风电机组的电磁暂态特性。


技术实现要素:

4.本发明的目的在于克服现有双馈风电场等值方法对风电场内电压差异以及低电压穿越控制考虑较少的不足,提供一种计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法,在保持风电场电磁暂态特性的前提下,对双馈风电场模型进行简化等值,为含双馈风电场系统的电磁暂态分析提供模型。
5.为了解决上述技术问题,本发明采用以下技术方案。
6.一种计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法,包括以下步骤:
7.(1)读取含双馈风电场系统数据并计算故障前稳态潮流:获取系统拓扑结构及导纳矩阵,各风电机组风速,风电场内部集电网络参数,计算故障前稳态时风电场各节点电压;
8.(2)计算考虑低电压穿越控制的故障稳态系统潮流:根据风电机组低电压穿越控制策略修改系统雅克比矩阵,计算短路稳态时风电机组机端电压以及输出的风电功率;
9.(3)基于机端电压跌落及故障前风电功率进行分群:根据风电机组电压跌落幅度
及故障前风电功率,估算风电机组功率恢复时间及风电场功率恢复曲线,据此对风电机组分群;
10.(4)求取等值风电机组及集电网络参数:根据风电功率守恒原则求取等值风电机组风速,根据注入电流不变的原则求取集电网络的等效参数;并对等值机功率恢复速度进行修正;
11.(5)校验双馈风电场等值模型的有效性:对含双馈风电场的算例系统进行等值计算,对等值前后风电场并网点的电压、电流、有功及无功功率进行比较,验证等值模型的有效性。
12.具体地,在步骤(1)中,读取含双馈风电场系统数据并计算故障前稳态潮流,所述的含双馈风电场系统数据包括:各同步发电机输出功率p
g,i
、机端电压v
g,i
,i=1,

,ng,ng为同步发电机总数;风电场中各风电机组风速v
w,j
,j=1,

,nw,nw为风电机组总数;线路阻抗、变压器变比及短路阻抗;根据上述数据构建m
×
m阶系统导纳矩阵y,m为节点总数,y
lk
=g
lk
+jb
lk
,l=1,

,m,k=1,

,m,g
lk
、b
lk
分别为节点l和节点k之间的电导和电纳;
13.将同步发电机所在节点视为pv节点;根据风速-功率曲线查表得出风电机组wj的稳态输出功率将风电机组所在节点视为有功为无功为0的pq节点;对m个节点进行编号,节点1为平衡节点,pq节点的编号为2,

,m-ng,pv节点编号为m-ng+1,m-ng+2,

,m;v、θ分别为各节点的电压及相角向量,潮流计算中相角及电压修正量分别记为δθ=[δθ2,δθ3,

,δθm]
t
、节点有功及无功不平衡量分别为δp=[δp2,δp2,

,δpm]
t
、设故障前系统雅克比矩阵为其中,根据功率修正方程计算故障前系统稳态潮流,得出正常运行时风电机组的机端电压及相角
[0014]
具体地,在步骤(2)中,计算考虑低电压穿越控制的故障稳态系统潮流;根据风电机组低电压穿越控制策略修改系统雅克比矩阵,计算短路稳态时风电机组机端电压以及输出的风电功率;其实现过程包括:
[0015]
读取含双馈风电场系统数据并计算故障前稳态潮流,所述的含双馈风电场系统数据包括:各同步发电机输出功率p
g,i
、机端电压v
g,i
,i=1,

,ng,ng为同步发电机总数;风电场中各风电机组风速v
w,j
,j=1,

,nw,nw为风电机组总数;线路阻抗、变压器变比及短路阻抗;根据上述数据构建m
×
m阶系统导纳矩阵y,m为节点总数,y
lk
=g
lk
+jb
lk
,l=1,

,m,k=1,

,m,g
lk
、b
lk
分别为节点l和节点k之间的电导和电纳;
[0016]
将同步发电机所在节点视为pv节点;根据风速-功率曲线查表得出风电机组wj的稳态输出功率将风电机组所在节点视为有功为无功为0的pq节点;对m个节点进行编号,节点1为平衡节点,pq节点的编号为2,

,m-ng,pv节点编号为m-ng+1,m-ng+2,

,m;v、θ分别为各节点的电压及相角向量,潮流计算中相角及电压修正量分别记为δθ=[δθ2,δθ3,

,δθm]
t
、节点有功及无功不平衡量分别为
δp=[δp2,δp2,

,δpm]
t
、设故障前系统雅克比矩阵为其中,根据功率修正方程计算故障前系统稳态潮流,得出正常运行时风电机组的机端电压及相角
[0017]
设t1时刻系统发生三相对地故障,故障稳态时风电机组的机端电压为则根据风电机组低电压穿越控制策略,风电机组输出有功功率及无功功率分别为:
[0018][0019]
其中,kd为无功功率系数;设故障状态下系统雅克比矩阵则h
(1)
=h
(0)
,j
(1)
=j
(0)
,n
(1)
=n
(0)
+δn,l
(1)
=l
(0)
+δl;其中,δn、δl按如下公式计算:
[0020][0021][0022]
式中,s=2,

,m;t=2,3,

,m-ng;e=2,3,

,m-ng;
[0023]
取各节点电压初始值为1,相角为0,代入功率修正方程求取δv、δθ,对节点电压v和相角θ进行修正,如此迭代直至max{||δv||,||δθ||}≤ε,||δv||和||δθ||分别为向量δv和δθ的最大分量的绝对值,ε为设定的小数,则认为潮流计算收敛,输出故障稳态各节点的电压v
(1)
及相角θ
(1)

[0024]
具体地,在步骤(3)中,基于机端电压跌落及故障前风电功率进行分群,根据风电机组电压跌落幅度及故障前风电功率,估算风电机组功率恢复时间及风电场功率恢复曲线,据此对风电机组分群;其实现过程包括:
[0025]
读取含双馈风电场系统数据并计算故障前稳态潮流,所述的含双馈风电场系统数据包括:各同步发电机输出功率p
g,i
、机端电压v
g,i
,i=1,

,ng,ng为同步发电机总数;风电场中各风电机组风速v
w,j
,j=1,

,nw,nw为风电机组总数;线路阻抗、变压器变比及短路阻抗;根据上述数据构建m
×
m阶系统导纳矩阵y,m为节点总数,y
lk
=g
lk
+jb
lk
,l=1,

,m,k=1,

,m,g
lk
、b
lk
分别为节点l和节点k之间的电导和电纳;
[0026]
将同步发电机所在节点视为pv节点;根据风速-功率曲线查表得出风电机组wj的稳态输出功率将风电机组所在节点视为有功为无功为0的pq节点;
对m个节点进行编号,节点1为平衡节点,pq节点的编号为2,

,m-ng,pv节点编号为m-ng+1,m-ng+2,

,m;v、θ分别为各节点的电压及相角向量,潮流计算中相角及电压修正量分别记为δθ=[δθ2,δθ3,

,δθm]
t
、节点有功及无功不平衡量分别为δp=[δp2,δp2,

,δpm]
t
、设故障前系统雅克比矩阵为其中,根据功率修正方程计算故障前系统稳态潮流,得出正常运行时风电机组的机端电压及相角
[0027]
设t1时刻系统发生三相对地故障,故障稳态时风电机组的机端电压为则根据风电机组低电压穿越控制策略,风电机组输出的有功功率及无功功率分别为:
[0028][0029]
其中,kd为无功功率系数;设故障状态下系统雅克比矩阵则h
(1)
=h
(0)
,j
(1)
=j
(0)
,n
(1)
=n
(0)
+δn,l
(1)
=l
(0)
+δl;其中,δn、δl按如下公式计算:
[0030][0031][0032]
式中,s=2,

,m;t=2,3,

,m-ng;e=2,3,

,m-ng;取各节点电压初始值为1,相角为0,代入功率修正方程求取δv、δθ,对节点电压v和相角θ进行修正,如此迭代直至max{||δv||,||δθ||}≤ε,||δv||和||δθ||分别为向量δv和δθ的最大分量的绝对值,ε为设定的小数,则认为潮流计算收敛,输出故障稳态各节点的电压v
(1)
及相角θ
(1)

[0033]
设t2时刻故障清除,故障后瞬间电压恢复正常值,风电机组有功电流参考值将以k
re
的速率逐渐恢复,而无功电流参考值则恢复到故障前的稳态值;故障后风电机组输出有功功率为:
[0034][0035]
其中,k
re
为风电机组有功功率恢复速率,t
re
,j为风电机组wj的有功恢复时间,
将nw台风电机组的有功恢复时间从小到大排序,设所得数组为t
re

,记t
re
=[0t
re

];则t
re
(k)时刻风电场并网点处的输出功率为:
[0036][0037]
式中,k=1,2,

,(nw+1),sw为双馈风电机组的额定容量;
[0038]
将平面上点记为ak;设a1a2与a1ak之间夹角为β
1,k
,k=2,

,(nw+1),则:
[0039][0040]
cosβ
1,k
反映了[t
re
(1),t
re
(k)]时段内风电场功率恢复曲线的线性度,cosβ
1,k
越接近1,该时间段内功率恢复特性的线性度越好,将该时间段内完成有功功率恢复的风电机组划分为同一群的偏差越小;将cosβ
1,k
≥0.999的风电机组划分为c1群,设c1群包含风电机组数量为n
c1
;类似地,对余下风电机组根据[t
re
(n
c1
+1),t
re
(k)]时段内风电场功率恢复曲线线性度分群,设与之间的夹角为β
2,k
,k=(n
c1
+2),

,(nw+1),则:
[0041][0042]
将满足cosβ
2,k
≥0.999的风电机组划分为c2群,设c2群包含风电机组数量为n
c2
;将余下风电机组归为c3群,设c3群包含风电机组数量为n
c3
,则n
c1
+n
c2
+n
c3
=nw。
[0043]
具体地,在步骤(4)中,求取等值风电机组及集电网络参数,根据风电功率守恒原则求取等值风电机组风速,根据注入电流不变的原则求取集电网络的等效参数;并对等值机功率恢复速度进行修正;其实现过程包括:
[0044]
读取含双馈风电场系统数据并计算故障前稳态潮流,所述的含双馈风电场系统数据包括:各同步发电机输出功率p
g,i
、机端电压v
g,i
,i=1,

,ng,ng为同步发电机总数;风电场中各风电机组风速v
w,j
,j=1,

,nw,nw为风电机组总数;线路阻抗、变压器变比及短路阻抗;根据上述数据构建m
×
m阶系统导纳矩阵y,m为节点总数,y
lk
=g
lk
+jb
lk
,l=1,

,m,k=1,

,m,g
lk
、b
lk
分别为节点l和节点k之间的电导和电纳;
[0045]
将同步发电机所在节点视为pv节点;根据风速-功率曲线查表得出风电机组wj的稳态输出功率将风电机组所在节点视为有功为无功为0的pq节点;对m个节点进行编号,节点1为平衡节点,pq节点的编号为2,

,m-ng,pv节点编号为m-ng+1,m-ng+2,

,m;v、θ分别为各节点的电压及相角向量,潮流计算中相角及电压修正量分别记为δθ=[δθ2,δθ3,

,δθm]
t
、节点有功及无功不平衡量分别为δp=[δp2,δp2,

,δpm]
t
、设故障前系统雅克比矩阵为其中,根据功率修正方程
计算故障前系统稳态潮流,得出正常运行时风电机组的机端电压及相角
[0046]
设t1时刻系统发生三相对地故障,故障稳态时风电机组的机端电压为则根据风电机组低电压穿越控制策略,风电机组输出的有功功率及无功功率分别为:
[0047][0048]
其中,kd为无功功率系数;设故障状态下系统雅克比矩阵则h
(1)
=h
(0)
,j
(1)
=j
(0)
,n
(1)
=n
(0)
+δn,l
(1)
=l
(0)
+δl;其中,δn、δl按如下公式计算:
[0049][0050][0051]
式中,s=2,

,m;t=2,3,

,m-ng;e=2,3,

,m-ng;取各节点电压初始值为1,相角为0,代入功率修正方程求取δv、δθ,对节点电压v和相角θ进行修正,如此迭代直至max{||δv||,||δθ||}≤ε,||δv||和||δθ||分别为向量δv和δθ的最大分量的绝对值,ε为设定的小数,则认为潮流计算收敛,输出故障稳态各节点的电压v
(1)
及相角θ
(1)

[0052]
设t2时刻故障清除,故障后瞬间电压恢复正常值,风电机组有功电流参考值将以k
re
的速率逐渐恢复,而无功电流参考值则恢复到故障前的稳态值;故障后风电机组输出有功功率为:
[0053][0054]
其中,k
re
为风电机组有功功率恢复速率,t
re,j
为风电机组wj的有功恢复时间,将nw台风电机组的有功恢复时间从小到大排序,设所得数组为t
re

,记t
re
=[0t
re

];则t
re
(k)时刻风电场并网点处的输出功率为:
[0055][0056]
式中,k=1,2,

,(nw+1),sw为双馈风电机组的额定容量;
[0057]
将平面上点记为ak;设a1a2与a1ak之间夹角为β
1,k
,k=2,

,(nw+1),则:
[0058][0059]
cosβ
1,k
反映了[t
re
(1),t
re
(k)]时段内风电场功率恢复曲线的线性度,cosβ
1,k
越接近1,该时间段内功率恢复特性的线性度越好,将该时间段内完成有功功率恢复的风电机组划分为同一群的偏差越小;将cosβ
1,k
≥0.999的风电机组划分为c1群,设c1群包含风电机组数量为n
c1
;类似地,对余下风电机组根据[t
re
(n
c1
+1),t
re
(k)]时段内风电场功率恢复曲线线性度分群,设与之间的夹角为β
2,k
,k=(n
c1
+2),

,(nw+1),则:
[0060][0061]
将满足cosβ
2,k
≥0.999的风电机组划分为c2群,设c2群包含风电机组数量为n
c2
;将余下风电机组归为c3群,设c3群包含风电机组数量为n
c3
,则n
c1
+n
c2
+n
c3
=nw;
[0062]
假设c1、c2及c3群对应的等值机分别为w
equ1
、w
equ2
及w
equ3
;以c1群为例,设w
equ1
的稳态功率为则根据风电功率守恒原则,根据风电机组风速-功率曲线应用查表法求出w
equ1
对应的风速v
equ1
;根据注入电流不变原则,等值机w
equ1
注入电流为conj(
·
)为共轭运算;w
equ1
所在节点电压为设风电场并网点电压为则等值机w
equ1
与风电场并网点的联接阻抗类似地,可以求出w
equ2
对应的风速v
equ2
及联接阻抗z
equ2
,以及w
equ3
对应的风速v
equ3
及联接阻抗z
equ3

[0063]
由于同一群内风电机组功率恢复时间有差异,故其风电功率曲线实际上为折线,斜率不断变化,而将该群等值1台等值机时,其功率恢复曲线的斜率为常数;为减小等值误差,根据风电场有功功率恢复曲线对等值机有功恢复速率逐段进行修正:
[0064][0065]
在故障后的初始阶段,w
equ1
、w
equ2
及w
equ3
的有功功率以k
re,1
的恢复速率上升;当w
equ1
恢复到故障前稳态功率时,将w
equ2
及w
equ3
的有功功率恢复速率设为k
re,2
;当w
equ2
恢复到故障前的稳态功率时,再将w
equ3
的有功功率恢复速率设为k
re,3

[0066]
具体地,在步骤(5)中,校验双馈风电场等值模型的有效性,对含双馈风电场的算例系统进行等值计算,对等值前后风电场并网点的电压、电流、有功及无功功率进行比较,验证等值模型的有效性;其实现过程包括:
[0067]
构建含双馈风电场的算例系统,分别计算故障前及故障中的系统稳态潮流,根据
故障期间风电机组的电压跌落及故障前风电功率,估算各风电机组功率恢复时间及风电场功率恢复曲线,根据风电场功率恢复曲线的线性度将风电机组分为3群,分别计算各等值机的等效风速v
equ
、与并网点的联接阻抗z
equ
,构建等值风电场;在含双馈风电场系统中设置三相对地短路故障,应用详细模型及等值模型分别进行仿真,记录风电场并网点注入的有功功率和无功功率,以及并网点的电压及电流,若等值前后的仿真曲线偏差较小,则表明等值模型具有有效性;若等值偏差较大,则增加分群数量,转步骤(3)重新进行风电场等值计算。
[0068]
与现有技术相比,本发明具有以下优点和有益效果:
[0069]
1.本发明提出的计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法,考虑风电机组的低电压穿越控制策略,对故障期间系统雅克比矩阵进行修改,在此基础上计算故障稳态的机端电压,并对故障后风电场有功功率恢复曲线进行估算,避免了繁复的时域仿真计算,极大地减小了计算量和分析速度;
[0070]
2.本发明提出了基于风电场有功功率恢复曲线线性度的风电机组分群方法,该方法兼顾了风电机组稳态功率及电压跌落幅度对风电机组电磁暂态特性的影响,与现有基于风速或短路电流波形相似性的分群方法相比,更全面地考虑了双馈风电场不同时间尺度的电磁暂态特性;同时根据功率守恒及注入电流不变原则计算等值风电机组的风速及联接阻抗,并根据风电场有功功率恢复曲线对等值风电机组的有功恢复速率进行逐段修正,减小了等值前后注入功率曲线偏差,提高了等值精度。
附图说明
[0071]
图1是本发明一个实施例的计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法流程图。
[0072]
图2是本发明一个实施例的含双馈风电场系统结构图。
[0073]
图3是本发明一个实施例实测的故障后风电场有功功率恢复曲线与理论估算曲线对比图,其中实线为实测风电场有功功率曲线,虚线为理论估算的风电场有功功率曲线。
[0074]
图4是本发明一个实施例的等值系统结构图。
[0075]
图5是本发明一个实施例的等值前后风电场输出有功及无功功率曲线对比图,其中细实线为详细模型有功功率曲线,细虚线为等值模型有功功率曲线;粗实线为详细模型无功功率曲线,粗虚线为等值模型无功功率曲线。
[0076]
图6是本发明一个实施例的等值前后风电并网点电压曲线对比图,其中实线为详细模型电压曲线,虚线为等值模型电压曲线。
[0077]
图7是本发明一个实施例的等值前后风电并网点a相短路电流对比图,其中实线为详细模型电流曲线,虚线为等值模型电流曲线。
具体实施方式
[0078]
下面结合附图和实施例对本发明做进一步详细说明。
[0079]
本发明的一种计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法,如图1所示,包括以下步骤:
[0080]
步骤一、读取含双馈风电场系统数据并计算故障前稳态潮流。所述的含双馈风电场系统数据包括:各同步发电机输出功率p
g,i
、机端电压v
g,i
,i=1,

,ng,ng为同步发电机总
数;风电场中各风电机组风速v
w,j
,j=1,

,nw,nw为风电机组总数;线路阻抗、变压器变比及短路阻抗;根据上述数据构建m
×
m阶系统导纳矩阵y,m为节点总数,y
lk
=g
lk
+jb
lk
,l=1,

,m,k=1,

,m,g
lk
、b
lk
分别为节点l和节点k之间的电导和电纳;
[0081]
将同步发电机所在节点视为pv节点;根据风速-功率曲线查表得出风电机组wj的稳态输出功率将风电机组所在节点视为有功为无功为0的pq节点;对m个节点进行编号,节点1为平衡节点,pq节点的编号为2,

,m-ng,pv节点编号为m-ng+1,m-ng+2,

,m;v、θ分别为各节点的电压及相角向量,潮流计算中相角及电压修正量分别记为δθ=[δθ2,δθ3,

,δθm]
t
、节点有功及无功不平衡量分别为δp=[δp2,δp2,

,δpm]
t
、设故障前系统雅克比矩阵为其中,根据功率修正方程计算故障前系统稳态潮流,得出正常运行时风电机组的机端电压及相角
[0082]
步骤二、计算考虑低电压穿越控制的故障稳态系统潮流:根据风电机组低电压穿越控制策略,修改系统雅克比矩阵,计算故障稳态时的系统潮流,求取风电机组机端电压以及输出的风电功率。
[0083]
设t1时刻系统发生三相对地故障,故障稳态时风电机组的机端电压为则根据风电机组低电压穿越控制策略,风电机组输出的有功功率及无功功率分别为:
[0084][0085]
其中,kd为无功功率系数;设故障状态下系统雅克比矩阵则h
(1)
=h
(0)
,j
(1)
=j
(0)
,n
(1)
=n
(0)
+δn,l
(1)
=l
(0)
+δl;其中,δn、δl按如下公式计算:
[0086][0087][0088]
式中,s=2,

,m;t=2,3,

,m-ng;e=2,3,

,m-ng;
[0089]
取各节点电压初始值为1,相角为0,代入功率修正方程求取δv、δθ,对节点电压v和相角θ进行修正,如此迭代直至max{||δv||,||δθ||}≤ε,||δv||和||δθ||分别为向量δv和δθ的最大分量的绝对值,ε为设定的小数,则认为潮流计
算收敛,输出故障稳态各节点的电压v
(1)
及相角θ
(1)

[0090]
步骤三、基于机端电压跌落及故障前风电功率进行分群:根据风电机组电压跌落幅度及故障前风电功率,估算风电机组功率恢复时间及风电场功率恢复曲线,根据风电场功率恢复曲线的线性度对风电机组分群。
[0091]
设t2时刻故障清除,故障后瞬间电压恢复正常值,风电机组有功电流参考值将以k
re
的速率逐渐恢复,而无功电流参考值则恢复到故障前的稳态值;故障后风电机组输出有功功率为:
[0092][0093]
其中,k
re
为风电机组有功功率恢复速率,t
re,j
为风电机组wj的有功恢复时间,将nw台风电机组的有功恢复时间从小到大排序,设所得数组为t
re

,记t
re
=[0t
re

];则t
re
(k)时刻风电场并网点处的输出功率为:
[0094][0095]
式中,k=1,2,

,(nw+1),sw为双馈风电机组的额定容量;
[0096]
将平面上点记为ak;设a1a2与a1ak之间夹角为β
1,k
,k=2,

,(nw+1),则:
[0097][0098]
cosβ
1,k
反映了[t
re
(1),t
re
(k)]时段内风电场功率恢复曲线的线性度,cosβ
1,k
越接近1,该时间段内功率恢复特性的线性度越好,将该时间段内完成有功功率恢复的风电机组划分为同一群的偏差越小;将cosβ
1,k
≥0.999的风电机组划分为c1群,设c1群包含风电机组数量为n
c1
;类似地,对余下风电机组根据[t
re
(n
c1
+1),t
re
(k)]时段内风电场功率恢复曲线线性度分群,设与之间的夹角为β
2,k
,k=(n
c1
+2),

,(nw+1),则:
[0099][0100]
将满足cosβ
2,k
≥0.999的风电机组划分为c2群,设c2群包含风电机组数量为n
c2
;将余下风电机组归为c3群,设c3群包含风电机组数量为n
c3
,则n
c1
+n
c2
+n
c3
=nw。
[0101]
步骤四、求取等值风电机组及集电网络参数:根据风电功率守恒原则求取等值风电机组风速,根据注入电流不变的原则求取集电网络的等效参数,根据估算的风电场有功功率恢复曲线对等值机功率恢复速度进行逐段修正。
[0102]
假设c1、c2及c3群对应的等值机分别为w
equ1
、w
equ2
及w
equ3
;以c1群为例,设w
equ1
的稳态功率为则根据风电功率守恒原则,根据风电机组风速-功率曲线应用查表法求出w
equ1
对应的风速v
equ1
;根据注入电流不变原则,等值机w
equ1
注入电流为
conj(
·
)为共轭运算;w
equ1
所在节点电压为设风电场并网点电压为则等值机w
equ1
与风电场并网点的联接阻抗类似地,可以求出w
equ2
对应的风速v
equ2
及联接阻抗z
equ2
,以及w
equ3
对应的风速v
equ3
及联接阻抗z
equ3

[0103]
由于同一群内风电机组功率恢复时间有差异,故其风电功率曲线实际上为折线,斜率不断变化,而将该群等值1台等值机时,其功率恢复曲线的斜率为常数;为减小等值误差,根据风电场有功功率恢复曲线对等值机有功恢复速率逐段进行修正:
[0104][0105]
在故障后的初始阶段,w
equ1
、w
equ2
及w
equ3
的有功功率以k
re,1
的恢复速率上升;当w
equ1
恢复到故障前稳态功率时,将w
equ2
及w
equ3
的有功功率恢复速率设为k
re,2
;当w
equ2
恢复到故障前的稳态功率时,再将w
equ3
的有功功率恢复速率设为k
re,3

[0106]
步骤五、校验双馈风电场等值模型的有效性:对含双馈风电场的算例系统进行等值计算,对等值前后风电场并网点的电压、电流、有功及无功功率进行比较,验证等值模型的有效性。
[0107]
构建含双馈风电场的算例系统,分别计算故障前及故障中的系统稳态潮流,根据故障期间风电机组的电压降落及故障前风电功率,估算各风电机组的功率恢复时间及风电场功率恢复曲线,根据风电场功率恢复特性的线性度将风电机组分为3群,分别计算各等值机的等效风速v
equ
、与并网点的联接阻抗z
equ
,构建等值风电场;在外部系统中设置三相对地短路故障,分别针对详细模型及等值模型进行仿真,记录风电场并网点注入的有功功率和无功功率,以及并网点的电压及电流,若等值前后的仿真曲线偏差较小,则可验证等值模型的有效性;若等值偏差较大,则增加分群数量,转步骤(3)重新进行风电场等值计算。
[0108]
试验验证:计及双馈风电场电压分布特性及低电压穿越控制的电磁暂态等值方法在含双馈风电场算例系统中的应用。
[0109]
本验证例将上述计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法应用于含双馈风电场算例系统,对其进行了电磁暂态等值计算。其步骤如下:
[0110]
步骤1、读取含双馈风电场系统数据并计算故障前稳态潮流。
[0111]
本实施例的含双馈风电场电力系统的结构如图2所示,双馈风电场含有18台双馈风力发电机,每台风力发电机容量为9mw,所接变压器t1~t
18
的容量为10.5mva,变比为25kv/575v;风电场内风电机通过3条支路联接到风电并网点,在节点4并入电网,t
pcc
的容量为100mva,变比为220kv/25kv。故障前稳态时各风机的风速及功率如表1所示,表1中风电功率标幺值以双馈风力发电机额定容量为基准;算例系统各线路阻抗如表2所示。系统中节点总数m=40,节点1为无穷大母线,系统中无同步发电机,故无pv节点,节点2-40均为pq节点。
[0112]
表1.故障前稳态时各风电机组风速、功率及机端电压
[0113][0114]
表2.线路阻抗
[0115][0116]
在此基础上形成导纳矩阵y,以及雅克比矩阵v、θ分别为各节点的电压及相角,潮流计算中相角及电压修正量分别记为δθ=[δθ2,δθ3,

,δθ
40
]
t
、δv=[δv2,δv3,

,δv
40
]
t
,节点有功及无功不平衡量分别为δp=[δp2,δp2,

,δp
40
]
t
、δq=[δq2,δq3,

,δq
40
]
t
;设故障前系统雅克比矩阵为其中,根据功率修正方程计算稳态时系统潮流,所得风电机组所在节点电压如表1所示。
[0117]
步骤2、计算考虑低电压穿越特性的故障稳态系统潮流。
[0118]
根据风电机组低电压穿越控制策略,修改系统雅克比矩阵,设t1=20.5s,此时系统发生三相对地故障,故障稳态时风电机组的机端电压为则根据风电机组低电压穿越特性,风电机组输出的有功功率及无功功率分别为:
[0119][0120]
取kd=1.5;设故障状态下系统雅克比矩阵则h
(1)
=h
(0)
,j
(1)
=j
(0)
,n
(1)
=n
(0)
+δn,l
(1)
=l
(0)
+δl;其中,δn、δl按如下公式计算:
[0121][0122]
[0123]
s=2,

,40;t=2,3,

,40;e=2,3,

,40;
[0124]
取各节点电压初始值为1,相角为0,代入功率修正方程求取δv、δθ,对节点电压v和相角θ进行修正,如此迭代直至max{||δv||,||δθ||}≤ε,ε取10-6
,迭代7次后,潮流计算收敛,输出故障稳态各节点的电压v
(1)
及相角θ
(1)
,结果如表3所示。
[0125]
表3.故障前后稳态时各风电机组机端电压对比
[0126][0127]
步骤3、基于机端电压跌落及故障前风电功率进行分群。
[0128]
根据风电机组电压跌落幅值及故障前风电功率,估算风电机组功率恢复时间及风电场功率恢复曲线。设t2=20.833s时故障清除,故障后瞬间电压恢复正常值,风电机组有功电流参考值将以k
re
=0.2pu/s的速率逐渐恢复,而无功电流参考值则恢复到故障前的稳态值;故障后风电机组wj输出的有功功率为:
[0129][0130]
其中,t
re,j
为风电机组wj的有功恢复时间,将nw台风电机组的有功恢复时间从小到大排序,设所得数组为t
re

,记t
re
=[0t
re

];则t
re
(k)时刻风电场并网点处的输出功率为:
[0131][0132]
式中,sw=9mva为双馈风电机组额定容量;将转换为以sb=100mva为基准的标幺值;在平面上绘出点记为ak,k=1,2,

,19,如图3所示;
[0133]
设a1a2与a1ak之间的夹角为β
1,k
,k=2,3,

,19,则:
[0134][0135]
cosβ
1,k
反映了从t
re
(1)到t
re
(k)时段风电场功率恢复曲线的线性度,cosβ
1,k
越接近1,该时间段内功率恢复特性的线性度越好,将该时间段内完成有功恢复的风电机组划分为同一群的偏差越小;将cosβ
1,k
≥0.999的风电机组划分为c1群,c1群包含风电机组数量为n
c1
,在本算例中,c1={w
17
,w5,w9,w
15
,w
10
,w3,w
18
},n
c1
=7;类似地,对余下风电机组可计算a8a9与a8ak之间的夹角β
2,k
,k=9,11,

,19,则:
[0136][0137]
根据[t
re
(8),t
re
(k)]时段内风电场功率恢复曲线线性度分群,将cosβ
2,k
≥0.999的风电机组划分为c2群,c2群包含风电机组数量为n
c2
,本算例中c2={w
13
,w4,w7,w
16
,w2,w6},n
c2
=6;余下的风电机组为c3群,c3={w8,w
12
,w
11
,w
14
,w1},n
c3
=5。
[0138]
步骤4、求取等值风电机组及集电网络参数。
[0139]
假设c1、c2及c3群对应的等值机分别为w
equ1
、w
equ2
及w
equ3
;根据风电功率守恒原则求取等值风电机组风速,根据注入电流不变的原则求取集电网络的等效参数;以c1群为例,设w
equ1
的稳态功率为则根据风电功率守恒原则,根据风电机组风速-功率曲线应用查表法求出w
equ1
对应的风速v
equ1
;根据注入电流不变原则,等值机w
equ1
注入电流为conj(
·
)为共轭运算;w
equ1
所在节点电压为设风电场并网点电压为则等值机w
equ1
与风电场并网点的联接阻抗类似地,可以求出w
equ2
对应的风速v
equ2
及联接阻抗z
equ2
,以及w
equ3
对应的风速v
equ3
及联接阻抗z
equ3
;在本算例中,等值风电机组的风速及联接阻抗如表4所示,所得等值系统结构图如图4所示。
[0140]
表4.等值风电机组的容量、风速及联接阻抗
[0141][0142]
由于同一群内风电机组功率恢复时间有差异,故其风电功率曲线实际上为折线,斜率不断变化,而将该群等值1台等值机时,其功率恢复曲线的斜率为常数;为减小等值误差,根据风电场有功功率恢复曲线对等值机有功恢复速率逐段进行修正:
[0143][0144]
计算表明,在本算例中,k
re,1
、k
re,2
、k
re,3
分别为0.1828pu/s、0.1741pu/s及0.1122pu/s,在故障后的初始阶段,w
equ1
、w
equ2
及w
equ3
的有功功率以0.1828pu/s的恢复速率上升;当w
equ1
恢复到故障前稳态功率时,w
equ2
及w
equ3
的有功功率恢复速率下降为0.1741pu/s;当w
equ2
恢复到故障前的稳态功率时,再将w
equ3
的有功功率恢复速率调整为0.1122pu/s。
[0145]
步骤5、校验双馈风电场等值模型的有效性。
[0146]
在matlab/simulink仿真软件平台下构建含双馈风电场的算例系统及其等值系统,在节点2设置三相对地短路故障,故障20.5秒开始,20个周期后故障清除,分别应用详细
模型及等值模型对算例系统的进行仿真计算,等值前后风电场输出有功及无功功率曲线对比图如图5所示;等值前后风电并网点电压曲线对比图如图6所示;等值前后风电并网点a相短路电流对比图如图7所示;由图可见,等值模型仿真结果与详细模型计算结果基本一致,较好地保留了双馈风电场的短路电流特性以及故障期间的有功及无功特性。
[0147]
综上所述,本发明的计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法,弥补了现有双馈风电场电磁暂态等值技术较少计及风电机组低电压穿越控制策略,以及对机端电压对双馈风力发电机电磁暂态特性的影响考虑不足的缺陷。在计及低电压穿越控制的条件下计算故障前及故障中系统稳态潮流,根据风电机组机端电压跌落幅值,并结合故障前、故障中及故障后风电机组的动态特性,快速估算风电机组的有功功率恢复曲线;然后根据风电场功率恢复曲线的线性度对风电机组进行分群,并基于功率守恒的原则求取等值风电机组的风速,根据注入电流不变的原则求取等值机组联接阻抗,根据风电场功率恢复曲线对等值风电机组的有功恢复速率进行逐段调整,有效地提高了等值精度。上述等值方法通过改进潮流计算的方法进行分群和等值,避免了繁复的电磁暂态仿真计算,减小了等值计算的运算量。该方法既考虑了故障瞬间电压突变引起的短路电流响应特性,又计及了电压跌落及低电压穿越控制对风电功率动态的影响,与现有的风电场电磁暂态等值方法相比,较好地保持了不同时间尺度的风电场电磁暂态特性。

技术特征:
1.一种计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法,其特征在于,包括以下步骤:步骤1、读取含双馈风电场系统数据并计算故障前稳态潮流;步骤2、计算考虑低电压穿越控制的故障稳态系统潮流:根据风电机组低电压穿越控制策略修改系统雅克比矩阵,计算短路稳态时风电机组机端电压以及输出的风电功率;步骤3、基于机端电压跌落及故障前风电功率进行分群:根据风电机组电压跌落幅度及故障前风电功率,估算风电机组功率恢复时间及风电场功率恢复曲线,根据功率恢复曲线的线性度对风电机组分群;步骤4、求取等值风电机组及集电网络参数:根据风电功率守恒原则求取等值风电机组风速,根据注入电流不变的原则求取集电网络的等效参数,根据功率恢复曲线对等值机功率恢复速度进行修正;步骤5、校验双馈风电场等值模型的有效性:对算例系统进行等值计算,对等值前后风电场并网点的电压、电流、有功及无功功率进行比较,验证等值模型的有效性;所述的含双馈风电场系统,其数据包括:各同步发电机输出功率p
g,i
、机端电压v
g,i
,i=1,

,n
g
,n
g
为同步发电机总数;风电场中各风电机组风速v
w,j
,j=1,

,n
w
,n
w
为风电机组总数;线路阻抗、变压器变比及短路阻抗;根据上述数据构建m
×
m阶系统导纳矩阵y,m为节点总数,y
lk
=g
lk
+jb
lk
,l=1,

,m,k=1,

,m,g
lk
、b
lk
分别为节点l和节点k之间的电导和电纳;将同步发电机所在节点视为pv节点;根据风速-功率曲线查表得出风电机组w
j
的稳态输出功率将风电机组所在节点视为有功为无功为0的pq节点;对m个节点进行编号,节点1为平衡节点,pq节点的编号为2,

,m-n
g
,pv节点编号为m-n
g
+1,m-n
g
+2,

,m;v、θ分别为各节点的电压及相角向量,节点有功及无功不平衡量分别为δp=[δp2,δp2,

,δp
m
]
t
、设故障前系统雅克比矩阵为其中,根据功率修正方程计算故障前系统稳态潮流,得出正常运行时风电机组的机端电压及相角2.根据权利要求1所述的一种计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法,其特征在于,在步骤2中,所述的计算考虑低电压穿越控制的故障稳态系统潮流中,根据风电机组低电压穿越控制策略修改系统雅克比矩阵,计算短路稳态时风电机组机端电压以及输出的风电功率,其具体过程包括:设t1时刻系统发生三相对地故障,故障稳态时风电机组的机端电压为则根据风电机组低电压穿越控制策略,风电机组输出有功功率及无功功率分别为:
其中,k
d
为无功功率系数;设故障状态下系统雅克比矩阵则h
(1)
=h
(0)
,j
(1)
=j
(0)
,n
(1)
=n
(0)
+δn,l
(1)
=l
(0)
+δl;其中,δn、δl按如下公式计算:+δl;其中,δn、δl按如下公式计算:式中,s=2,

,m;t=2,3,

,m-n
g
;e=2,3,

,m-n
g
;取各节点电压初始值为1,相角为0,代入功率修正方程求取δv、δθ,对节点电压v和相角θ进行修正,如此迭代直至max{||δv||,||δθ||}≤ε,||δv||和||δθ||分别为向量δv和δθ的最大分量的绝对值,ε为设定的小数,则认为潮流计算收敛,输出故障稳态各节点的电压v
(1)
及相角θ
(1)
。3.根据权利要求1所述的一种计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法,其特征在于,在步骤3中,所述的基于机端电压跌落及故障前风电功率进行分群中,根据风电机组电压跌落幅度及故障前风电功率,估算风电机组功率恢复时间及风电场功率恢复曲线,据此对风电机组分群,其具体过程包括:设t1时刻系统发生三相对地故障,故障稳态时风电机组的机端电压为则根据风电机组低电压穿越控制策略,风电机组输出有功功率及无功功率分别为:其中,k
d
为无功功率系数;设故障状态下系统雅克比矩阵则h
(1)
=h
(0)
,j
(1)
=j
(0)
,n
(1)
=n
(0)
+δn,l
(1)
=l
(0)
+δl;其中,δn、δl按如下公式计算:+δl;其中,δn、δl按如下公式计算:
式中,s=2,

,m;t=2,3,

,m-n
g
;e=2,3,

,m-n
g
;取各节点电压初始值为1,相角为0,代入功率修正方程求取δv、δθ,对节点电压v和相角θ进行修正,如此迭代直至max{||δv||,||δθ||}≤ε,||δv||和||δθ||分别为向量δv和δθ的最大分量的绝对值,ε为设定的小数,则认为潮流计算收敛,输出故障稳态各节点的电压v
(1)
及相角θ
(1)
;设t2时刻故障清除,故障后瞬间电压恢复正常值,风电机组有功电流参考值将以k
re
的速率逐渐恢复,而无功电流参考值则恢复到故障前的稳态值;故障后风电机组输出有功功率为:其中,k
re
为风电机组有功功率恢复速率,t
re,j
为风电机组w
j
的有功恢复时间,将n
w
台风电机组的有功恢复时间从小到大排序,设所得数组为t
re

,记t
re
=[0t
re

];则t
re
(k)时刻风电场并网点处的输出功率为:式中,k=1,2,

,(n
w
+1),s
w
为双馈风电机组的额定容量;将平面上点记为a
k
;设a1a2与a1a
k
之间夹角为β
1,k
,k=2,

,(n
w
+1),则:cosβ
1,k
反映了[t
re
(1),t
re
(k)]时段内风电场功率恢复曲线的线性度,cosβ
1,k
越接近1,该时间段内功率恢复特性的线性度越好,将该时间段内完成有功功率恢复的风电机组划分为同一群的偏差越小;将cosβ
1,k
≥0.999的风电机组划分为c1群,设c1群包含风电机组数量为n
c1
;类似地,对余下风电机组根据[t
re
(n
c1
+1),t
re
(k)]时段内风电场功率恢复曲线线性度分群,设与之间的夹角为β
2,k
,k=(n
c1
+2),

,(n
w
+1),则:将满足cosβ
2,k
≥0.999的风电机组划分为c2群,设c2群包含风电机组数量为n
c2
;将余下风电机组归为c3群,设c3群包含风电机组数量为n
c3
,则n
c1
+n
c2
+n
c3
=n
w
。4.根据权利要求1所述的一种计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法,其特征在于,在步骤4中,所述的求取等值风电机组及集电网络参数中,根据风电功率守恒原则求取等值风电机组风速,根据注入电流不变的原则求取集电网络的等效参数,根据功率恢复曲线对等值机功率恢复速度进行修正,其具体过程包括:设t1时刻系统发生三相对地故障,故障稳态时风电机组的机端电压为
则根据风电机组低电压穿越控制策略,风电机组输出有功功率及无功功率分别为:其中,k
d
为无功功率系数;设故障状态下系统雅克比矩阵则h
(1)
=h
(0)
,j
(1)
=j
(0)
,n
(1)
=n
(0)
+δn,l
(1)
=l
(0)
+δl;其中,δn、δl按如下公式计算:+δl;其中,δn、δl按如下公式计算:式中,s=2,

,m;t=2,3,

,m-n
g
;e=2,3,

,m-n
g
;取各节点电压初始值为1,相角为0,代入功率修正方程求取δv、δθ,对节点电压v和相角θ进行修正,如此迭代直至max{||δv||,||δθ||}≤ε,||δv||和||δθ||分别为向量δv和δθ的最大分量的绝对值,ε为设定的小数,则认为潮流计算收敛,输出故障稳态各节点的电压v
(1)
及相角θ
(1)
;设t2时刻故障清除,故障后瞬间电压恢复正常值,风电机组有功电流参考值将以k
re
的速率逐渐恢复,而无功电流参考值则恢复到故障前的稳态值;故障后风电机组输出有功功率为:其中,k
re
为风电机组有功功率恢复速率,t
re,j
为风电机组w
j
的有功恢复时间,将n
w
台风电机组的有功恢复时间从小到大排序,设所得数组为t
re

,记t
re
=[0t
re

];则t
re
(k)时刻风电场并网点处的输出功率为:式中,k=1,2,

,(n
w
+1),s
w
为双馈风电机组的额定容量;将平面上点记为a
k
;设a1a2与a1a
k
之间夹角为β
1,k
,k=2,

,(n
w
+1),则:
cosβ
1,k
反映了[t
re
(1),t
re
(k)]时段内风电场功率恢复曲线的线性度,cosβ
1,k
越接近1,该时间段内功率恢复特性的线性度越好,将该时间段内完成有功功率恢复的风电机组划分为同一群的偏差越小;将cosβ
1,k
≥0.999的风电机组划分为c1群,设c1群包含风电机组数量为n
c1
;类似地,对余下风电机组根据[t
re
(n
c1
+1),t
re
(k)]时段内风电场功率恢复曲线线性度分群,设与之间的夹角为β
2,k
,k=(n
c1
+2),

,(n
w
+1),则:将满足cosβ
2,k
≥0.999的风电机组划分为c2群,设c2群包含风电机组数量为n
c2
;将余下风电机组归为c3群,设c3群包含风电机组数量为n
c3
,则n
c1
+n
c2
+n
c3
=n
w
;假设c1、c2及c3群对应的等值机分别为w
equ1
、w
equ2
及w
equ3
;以c1群为例,设w
equ1
的稳态功率为则根据风电功率守恒原则,根据风电机组风速-功率曲线应用查表法求出w
equ1
对应的风速v
equ1
;根据注入电流不变原则,等值机w
equ1
注入电流为conj(
·
)为共轭运算;w
equ1
所在节点电压为设风电场并网点电压为则等值机w
equ1
与风电场并网点的联接阻抗类似地,可以求出w
equ2
对应的风速v
equ2
及联接阻抗z
equ2
,以及w
equ3
对应的风速v
equ3
及联接阻抗z
equ3
;由于同一群内风电机组功率恢复时间有差异,故其风电功率曲线实际上为折线,斜率不断变化,而将该群等值1台等值机时,其功率恢复曲线的斜率为常数;为减小等值误差,根据风电场有功功率恢复曲线对等值机有功恢复速率逐段进行修正:在故障后的初始阶段,w
equ1
、w
equ2
及w
equ3
的有功功率以k
re,1
的恢复速率上升;当w
equ1
恢复到故障前稳态功率时,将w
equ2
及w
equ3
的有功功率恢复速率设为k
re
,2;当w
equ2
恢复到故障前的稳态功率时,再将w
equ3
的有功功率恢复速率设为k
re,3
。5.根据权利要求1所述的一种计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法,其特征在于,在步骤5中,所述的校验双馈风电场等值模型的有效性中,对算例系统进行等值计算,对等值前后风电场并网点的电压、电流、有功及无功功率进行比较,验证等值模型的有效性,具体实现过程包括:构建含双馈风电场的算例系统,分别计算故障前及故障中的系统稳态潮流,根据故障期间风电机组的电压跌落及故障前风电功率,估算各风电机组的功率恢复时间及风电场功
率恢复曲线,根据风电场功率恢复特性的线性度将风电机组分为3群,分别计算各等值机的等效风速v
equ
、与并网点的联接阻抗z
equ
,构建等值风电场;在外部系统中设置三相对地短路故障,分别针对详细模型及等值模型进行仿真,记录风电场并网点注入的有功功率和无功功率,以及并网点的电压及电流,若等值前后的仿真曲线偏差较小,则可验证等值模型的有效性;若等值偏差较大,则增加分群数量,重新进行风电场等值计算。

技术总结
本发明公开了一种计及双馈风电场电压分布特性和低电压穿越控制的电磁暂态等值方法,包括:读取含双馈风电场系统数据并计算故障前稳态潮流;计算考虑低电压穿越特性的故障稳态系统潮流;基于机端电压跌落及故障前风电功率进行分群;求取等值风电机组及集电网络参数;校验双馈风电场等值模型的有效性。本发明在计及低电压穿越特性的条件下计算故障稳态期间的系统潮流,快速估算风电场的有功功率恢复曲线,在此基础上进行风电机组分群及风电场等值。与现有风电场电磁暂态等值方法相比,本发明较好地保持了不同时间尺度的风电场电磁暂态特性,提高了等值模型的精度,且避免了繁复的电磁暂态仿真计算,减小了等值计算的运算量。量。量。


技术研发人员:周海强 崔晓丹 杨心刚 曹博源 许剑冰 吴家龙 高超 王爽
受保护的技术使用者:河海大学 南瑞集团有限公司
技术研发日:2023.05.25
技术公布日:2023/8/13
版权声明

本文仅代表作者观点,不代表航空之家立场。
本文系作者授权航家号发表,未经原创作者书面授权,任何单位或个人不得引用、复制、转载、摘编、链接或以其他任何方式复制发表。任何单位或个人在获得书面授权使用航空之家内容时,须注明作者及来源 “航空之家”。如非法使用航空之家的部分或全部内容的,航空之家将依法追究其法律责任。(航空之家官方QQ:2926969996)

飞行汽车 https://www.autovtol.com/

分享:

扫一扫在手机阅读、分享本文

相关推荐