一种变循环发动机加速过程最优控制方法

文档序号:25296435发布日期:2021-06-04 11:31阅读:155来源:国知局
一种变循环发动机加速过程最优控制方法

1.本发明涉及变循环发动机控制技术领域,尤其涉及一种变循环发动机加速过程最优控制方法。


背景技术:

2.现代战争要求先进战斗机具备长航程亚声速巡航的能力,同时在作战时又要具备快速反应能力,未来变循环发动机将向长巡航里程、高推重比、宽工作范围三个方向不断发展。通过研究常规发动机速度特性,研究者发现超声速状态下涡喷发动机具有较高的单位推力和较低的单位燃油消耗率而亚声速状态下大涵道比涡扇发动机具有较低的单位燃油消耗率。考虑现代战争对战斗机推进系统的性能要求,涡扇发动机更加适合亚声速飞行,而涡喷发动机更适合超声速飞行。因此,便有了性能更好的变循环发动机。在发动机不同的工作状态下,通过采用调节特征部件的几何形状、物理位置或尺寸大小等不同的技术手段,将涡扇和涡喷两种不同的变循环发动机的性能优势集中一体,从而保证变循环发动机在亚声速巡航状态下以涡扇发动机类似构型工作,从而获得较高的经济性,在超声速作战状态下以涡喷发动机类似构型工作,从而获得持续可靠的高单位推力,达到了将涡扇、涡喷发动机的性能优势融为一体的目的,使变循环发动机在发动机工作全过程中均具有优良的性能。
3.变循环发动机是飞机的心脏,是衡量一个国家航空事业发展水平的重要指标之一,因此对强化动力系统的研究对提升国家航空技术整体水平具有重要意义。由于变循环发动机的工作过程复杂多变,且具有强非线性、多控制变量、时变、复杂的结构特点,因此,对变循环发动机控制问题的研究比一般控制系统更为困难。
4.现代战机对飞机的机动性要求非常高,良好的机动性就要求变循环发动机具有良好的加速性能。加速过程控制是变循环发动机过渡态控制的一种,相较于变循环发动机起动、接通/切断加力、减速控制,加速过程控制对变循环发动机以及飞机性能的影响更为明显。变循环发动机的加速过程直接影响战斗机的重要飞行指标(如:战斗机加速、爬升和紧急着陆复飞等等),因此,研究变循环发动机加速过程的最优控制,改善变循环发动机加速性能具有重要意义。
5.国内外在变循环发动机加速过程的最优控制研究中虽然取得了一定成果,但也存在许多尚未解决的技术难题或待改进之处。比如,单纯形法具有超线性收敛速度,迭代次数少,但是基本单纯形法对初值敏感,易陷入局部最优解,不适宜应用于复杂的变循环发动机加速过程寻优控制中。


技术实现要素:

6.为解决现有技术存在的问题,本发明提出一种变循环发动机加速过程最优控制方法,对单纯形法进行改进,并将改进的单纯形法应用于变循环发动机加速过程寻优控制,实现变循环发动机加速过程的最优控制,提高变循环发动机的加速过程性能,提高飞机的机动性。
7.本发明的技术方案为:
8.首先建立变循环发动机的非线性数学模型,然后以改进单纯形法来进行变循环发动机加速过程寻优,以实现某型航空涡扇变循环发动机加速过程最优。
9.所述一种变循环发动机加速过程最优控制方法,其特征在于:第一步建立变循环发动机的非线性数学模型;第二步根据变循环发动机加速过程确定相应的目标函数和约束函数;第三步以改进单纯形法优化计算;第四步输出最优控制变量给变循环发动机。
10.所述一种变循环发动机加速过程最优控制方法,其特征在于:所述改进单纯形法是在基本的单纯形法上进行改进,主要对反射中心进行改进,并且添加了“顶点平移”策略。其基本思想是首先对n+1个顶点的目标函数值进行最优搜索,确定平移方向;然后将单纯形中心点向目标函数值最好顶点方向适度平移,在迭代的末端过程,n+1个顶点与中心点近乎重合,依靠顶点自身的迭代已经可以很好的逼近最优解,此时如果继续进行顶点平移反而会添加扰动,增加迭代次数。因此,当迭代误差小于进行平移操作的误差阈值时,则放弃顶点平移操作。
11.所述变循环发动机的非线性数学模型为
12.y=f(x)
13.其中为控制输入向量,包括模式选择活门msv打开程度msv,调节主燃油流量w
f
、尾喷管面积a9、风扇导叶角度dvgl和压气机导叶角度dvgh,为输出向量,包括燃油消耗率sfc和变循环发动机推力f,f(
·
)为产生系统输出的非线性向量函数。
14.所述加速过程考虑的约束条件有:涡轮前温度不超温、高压压气机不喘振、高压转子不超转、风扇不超转、燃烧室不富油熄火、主燃烧室供油量不超过其最大供油量等等。优化问题的数学描述如下:
[0015][0016]
其中控制变量x=[msv,w
f
,a9,dvgl,dvgh]
t
,以上各个变量均在相应的变化范围之内取初值。
[0017]
采用线性加权法将多目标函数转化为单目标函数,来确定寻优目标函数。即
[0018][0019]
对上式进行离散化和归一化处理。这样处理的目的是为了消除目标函数中各参数
量纲和量值变化范围的不同对优化结果的影响。最终的寻优目标函数可以写成以下形式:
[0020][0021]
上式中,ω
a
和ω
b
为相应目标函数的权重系数,满足ω
a
≥0,ω
b
≥0,其大小反映相应的寻优目标函数在多目标优化问题中的重要程度。
[0022]
参照目标函数的形式,对变循环发动机约束条件也进行离散化和归一化处理:
[0023][0024][0025][0026][0027]
以上g
i
(x)(i=1,2,...,11)构成约束函数矩阵g(x),考虑约束条件后,目标函数可化为:
[0028][0029]
其中ω=[ω1,ω2,ω3,ω4,ω5,ω6,ω7,ω8,ω9,ω
10

11
]为约束函数的权重调整系数矩阵,其中ω1,ω2,ω3,ω4,ω5,ω6,ω7,ω8,ω9,ω
10

11
为对应约束条件可调整权重系数,ω
·
g(x)的设计用于满足变循环发动机的约束条件。
[0030]
所述改进单纯形法的算法流程为
[0031]
(1)初始化。对n维非线性模型,给定初始顶点x0,其余顶点按下式计算,可构造边长相等的正规单纯形,且k=0。
[0032]
x
(i)
=x
(1)
+a
×
[q,

,q,p
(i)
,q

,q]
t
(i=2,

,n+1)
[0033]
其中,p
(i)
表示第i个元素为p,
[0034][0035][0036]
a是单纯形边长;
[0037]
(2)计算各顶点的目标函数值f(x
(i)
),确定最优点和最差点满足的要求
[0038][0039][0040]
并计算反射中心点和收敛误差err;
[0041][0042][0043]
其中
[0044]
(3)如果收敛误差大于平移操作误差阈值ε
k
,即err>ε
k
,按下式向将单纯形所有顶点向最优点进行平移;否则,执行步骤4);
[0045][0046]
其中,λ∈(0,0.2)是平移系数。
[0047]
(4)进行单纯形反射、收缩、扩张、减小棱长操作计算;
[0048]
(5)如果收敛误差err大于迭代精度ε
e
,k=k+1,返回(2);否则,满足精度要求,迭代计算结束。
[0049]
进一步的,所述控制变量为模式选择活门msv打开程度msv,调节主燃油流量w
f
、尾喷管面积a9、风扇导叶角度dvgl和压气机导叶角度dvgh。
[0050]
有益效果
[0051]
与现有技术相比较,本发明的一种变循环发动机加速过程最优控制方法,对单纯形法进行改进,对反射中心进行了改进,并添加了“顶点平移”策略。顶点平移操作可加快搜索速度,减少迭代次数,快速收敛到最优解。并将改进的单纯形法应用于变循环发动机加速过程寻优控制,实现变循环发动机加速过程的最优控制,在保证变循环发动机安全工作前提下,缩短变循环发动机加速时间,有效改善变循环发动机加速性能,提高飞机的机动性。
[0052]
本发明的附加方面和优点将在下面的描述中部分给出,部分将从下面的描述中变得明显,或通过本发明的实践了解到。
附图说明
[0053]
本发明的上述和/或附加的方面和优点结合下面附图对实施例的描述中将变得明显和容易理解,其中:
[0054]
图1是本发明变循环发动机结构示意图;
[0055]
图2是本发明变循环发动机可调参数示意图;
[0056]
图3是本发明变循环发动机可调部件示意图;
[0057]
图4是本发明变循环发动机双外涵模式气流分布图;
[0058]
图5是本发明变循环发动机单外涵模式气流分布图;
[0059]
图6是本发明变循环发动机加速过程寻优控制流程图;
[0060]
图7是本发明单纯形法的寻优操作示意图;
[0061]
图8是本发明改进的单纯形法流程图。
具体实施方式
[0062]
本发明解决的问题是变循环发动机的加速过程寻优控制。变循环发动机寻优问题
就是为了使变循环发动机的加速过程达到最优,选取最优控制方法寻找一组最优控制量(模式选择活门msv打开程度msv,主燃油流量w
f
、尾喷管面积a9、风扇导流叶片角度dvgl、压气机导流叶片角度dvgh)。
[0063]
以某型航空涡扇变循环发动机非线性数学模型为研究对象,建立加速过程的相应目标函数,利用优化算法对变循环发动机进行优化计算,即可得到加速过程的满足最优性能指标的最优控制变量,在保证变循环发动机安全工作前提下,缩短变循环发动机加速时间,有效改善变循环发动机加速性能。
[0064]
本发明在总结前人成果的基础上,根据变循环发动机的特点,对单纯形法进行改进,并应用于变循环发动机寻优控制中。
[0065]
1、变循环发动机工作原理
[0066]
本发明以带有核心驱动风扇级(cdfs)的双外涵变循环发动机为主要研究对象,其主要结构如图1所示,包含的主要部件有进气道、风扇、核心驱动风扇级、高压压气机、燃烧室、高压涡轮、低压涡轮、混合室、加力燃烧室、尾喷管。相比普通双轴涡扇变循环发动机,其显著的结构特点是在风扇和高压压气机之间增加了cdfs,同时在风扇和cdfs后分别设置副外涵和主外涵。在变循环发动机的不同工作状态下,通过改变cdfs的导叶角度,可以大幅度调节变循环发动机外涵道和核心机的空气流量,进而调节变循环发动机内外涵空气流量、涵道比和增压比等循环参数,使变循环发动机的热力循环调节更加灵活。
[0067]
相比普通双轴涡扇变循环发动机,变循环发动机具有更多的可调部件。带有cdfs部件的变循环发动机主要有8个可调部件,具体如图2所示,可调部件示意图如图3所示。
[0068]
与传统变循环发动机相比,变循环发动机的性能优势主要体现在由于其可调部件增加,通过改变可调部件参数,调节变循环发动机在工作过程中的气动热力循环,在保证推力基本不变时显著降低单位燃油消耗率,大大提高变循环发动机的经济效益,同时可调部件的增加,使控制系统的调节过程更加灵活,风扇、压气机等部件的稳定裕度将大大提高。
[0069]
变循环发动机具有单/双外涵两种典型的工作模式,它们通过模式选择活门msv、fvabi、rvabi等可变阀门实现切换。当msv完全打开时,气流经风扇后被分为两部分,一股气流流入副外涵,这部分气流最终在主外涵出口截面与主外涵气流有效掺混,流入总外涵。另一股气流流入cdfs,这股气流部分经rvabi被引入总外涵,其余气流将流入核心机。由于尾端涵道和rvabi的存在,总外涵气流在出口处会分为两部分,一股气流直接通过尾端涵道流入尾喷管,另一股气流则会进入混合室,与通过核心机的气流进行掺混后经加力燃烧室燃烧后,流入尾喷管,具体气流分布如图4所示。上述工作过程中,主外涵和副外涵均有气流通过,故被命名为双外涵模式。
[0070]
当模式选择活门msv完全关闭时,流经风扇的气流全部流入cdfs,风扇以压气机模式工作,副外涵不再有气流通过,这一过程被命名为单外涵工作模式,其具体气流分布如图5所示。
[0071]
变循环发动机在不同的工作模式下进行切换时,内部热力循环状态会随之发生改变。为保证变循环发动机能够持续保持稳定可靠地工作,平稳实现单双外涵模式转换,在模式切换过程中应该满足以下基本条件:
[0072]
(1)风扇进口流量基本保持不变;
[0073]
(2)风扇增压比基本保持不变;
[0074]
(3)核心驱动风扇级的增压比随切换过程平稳变化;
[0075]
(4)涵道比随msv位移的改变而平稳变化;
[0076]
(5)保证回流裕度始终大于0,即不存在气流绕cdfs的倒流;
[0077]
(6)避免发生持续的超温、超转现象,避免喘振现象。
[0078]
为了满足上述条件,在调节msv位移时,应该配合调节其他可调部件参数,模式选择活门msv打开程度可以表征变循环发动机的工作模式。目前已证实切实可行的模式切换调节策略是:在单外涵至双外涵的模式切换过程中,通过调节msv位移,增大副外涵进口截面面积,为避免风扇压比的大幅度降低,需要配合减小cdfs进口导流叶片角度α
i
,同时减小可调涡轮导向器角度α
t
。由双外涵至单外涵的模式切换过程,调节策略相反。变循环发动机在不同工作模式工作时,为获得理想的涵道比同时保证气流不发生喘振或其他非正常工作状态,需要调节cdfs导流叶片角度α
i
以改变内涵空气流量,使之与变循环发动机工作状态匹配。
[0079]
2、由于变循环发动机加速过程寻优控制需要依据变循环发动机当前工作状态参数做出控制决策,因此,进行加速过程最优控制方法研究时,通常以变循环发动机数学模型取代真实的变循环发动机。由于变循环发动机的建模技术已经非常成熟,这里不再赘述,直接给出建立的变循环发动机非线性模型
[0080]
y=f(x)
[0081]
其中为控制输入向量,包括模式选择活门msv打开程度msv,调节主燃油流量w
f
、尾喷管面积a9、风扇导叶角度dvgl和压气机导叶角度dvgh,为输出向量,包括燃油消耗率sfc和变循环发动机推力f,f(
·
)为产生系统输出的非线性向量函数。
[0082]
3、改进单纯形法的设计
[0083]
变循环发动机的动态性能寻优控制中的最短响应时间控制模式是指在保证变循环发动机安全工作前提下,缩短变循环发动机加速时间。最短响应时间控制模式通常用于变循环发动机加速过程,有效改善变循环发动机的加速性能。变循环发动机加速过程寻优控制流程如图6所示,其基本思想是:首先以所建立的涡扇变循环发动机非线性数学模型为基础,在确保变循环发动机安全运行的前提下,以缩短变循环发动机加速时间为优化目标,然后寻求最优的控制计划,充分挖掘变循环发动机的性能潜力以达到优化的目的。由于变循环发动机具有强非线性性、高复杂性等特点,运用传统的优化方法难以同时提高优化精度和速度,所以必须采用更有效的优化算法来解决此问题。
[0084]
单纯形法不需要计算出目标函数的梯度,能加速计算收敛速度,缩减计算时间,但是常规单纯形法缺点在于太过依赖于初始值,迭代过程过于繁琐,不能直接用在变循环发动机加速过程的寻优控制中。因此,对单纯形法进行改进,并将改进的单纯形法用于变循环发动机加速过程的寻优控制中。
[0085]
单纯形是指n维空间具有n+1个顶点的凸多面体。单纯形法是直接法中计算比较简单,几何概念比较清晰的一种算法。单纯形法是对n维空间的(n+1)个点上的目标函数值进行比较,去掉其中最坏的点,代之以新的顶点、新的点与前面余下的点又构成一个新的单纯形。每次把坏的点去掉,把好的点留下来,逐步地剔除目标函数最优点不可能存在的空间,
直到将最优点包容在单纯形内,再将单纯形的几何尺寸缩小到小于收敛准则,就完成了最优点的搜索。
[0086]
本发明设计一种改进的单纯形法,主要对反射中心进行改进,并且添加了“顶点平移”策略。其基本思想是首先对n+1个顶点的目标函数值进行最优搜索,确定平移方向;然后将单纯形中心点向目标函数值最好顶点方向适度平移,在迭代的末端过程,n+1个顶点与中心点近乎重合,依靠顶点自身的迭代已经可以很好的逼近最优解,此时如果继续进行顶点平移反而会添加扰动,增加迭代次数。因此,当迭代误差小于进行平移操作的误差阈值时,则放弃顶点平移操作。
[0087]
单纯形的不同形成方法,就形成了各种单纯形法。通常有正规单纯形法、特殊单纯形法、long系数表法、利用均匀设计构造初始单纯形法等等。本发明选用正规单纯形法构造初始单纯形。正规单纯形指的是n+1个顶点间的距离都相等的单纯形。构造正规单纯形的方法如下:
[0088]
设n维单纯形的一个顶点为
[0089]
x
(0)
=[a
1 a2ꢀ…ꢀ
a
n
]
t
[0090]
其余n个顶点分别取为
[0091]
x
(1)
=[a1+p a2+q
ꢀ…ꢀ
a
n
+q]
t
[0092]
x
(2)
=[a1+q a2+p
ꢀ…ꢀ
a
n
+q]
t
[0093]
………………………………
[0094]
x
(n)
=[a1+q a2+q
ꢀ…ꢀ
a
n
+p]
t
[0095]
式中:
[0096][0097]
a为单纯形的边长,根据具体情况确定。
[0098]
改进单纯形法的搜寻方法包含4种操作,即反射、延伸、收缩和缩小棱长,如图7所示。
[0099]
本发明选择正规单纯形作为初始单纯形。单纯形法的迭代过程:
[0100]

[0101][0102]
是搜索的第k阶段(k=0,1,

)上e
n
中的第i个顶点,在点处的目标函数值为此外,需要在单纯形内标记给出f(x)的极大值和极小值的x向量,定义当前最差点的
[0103][0104]
和当前最好点的
[0105][0106]
以及当前次差点的
[0107][0108]
e
n
中寻找使目标函数f(x)有更好值的一个顶点的过程中包含四步操作:
[0109]
反射:
[0110]
求在和的联线上关于反射点,反射长度是和之间长度的α倍,即
[0111][0112]
式中,α为给定的反射系数,一般取α=1,此时称为标准反射点;称为反射中心;是在第k阶段上使目标函数f(x)的(n+1)个值中最大的f(x)时对应的顶点。
[0113]
本发明对反射中心进行了改进,新算法的反射中心点是剔除最差点以后剩下点的加权平均值,可由下式给出:
[0114][0115]
其中
[0116]
扩张:
[0117]
如果则在和的联线上按下式
[0118][0119]
计算扩张点。其中γ>1是扩张系数,一般取2.0。
[0120]
如果则用代替且以k=k+1由第一步(即反射)下同,继续进行。否则用代替且以k=k+1由第一步继续进行。
[0121]
收缩:
[0122]
如果对所有的i≠h有按下式计算:
[0123][0124]
其中0<β<1是收缩系数。
[0125]
如果则
[0126][0127]
计算若可用代替且返回第一步,在第k+1阶段上继续搜索。否则转下步。
[0128]
缩小:
[0129]
如果反射失败,则按下式计算:
[0130][0131]
再返回第一步,进行第k+1阶段搜索。
[0132]
每一次得到一个新的单纯形时,是否已得到满意的结果,都要进行检验。设预先给
定精度ε,则按下列收敛准则
[0133][0134]
检验之。若上式成立,停止迭代,输出及否则使k=k+1,返回第一步。以上四个操作如图7所示。
[0135]
本发明在常规单纯形法的基础上,对反射中心进行了改进,并添加了“顶点平移”策略。顶点平移操作可加快搜索速度,减少迭代次数,快速收敛到最优解。改进单纯形法算法流程如图8所示:
[0136]
(1)初始化。对n维非线性模型,给定初始顶点x0,其余顶点按下式计算,可构造边长相等的正规单纯形,且k=0。
[0137]
x
(i)
=x
(1)
+a
×
[q,

,q,p
(i)
,q

,q]
t
(i=2,

,n+1)
[0138]
其中,p
(i)
表示第i个元素为p,
[0139][0140][0141]
a是单纯形边长;
[0142]
(2)计算各顶点的目标函数值f(x
(i)
),确定最优点和最差点满足的要求
[0143][0144][0145]
并计算反射中心点和收敛误差err;
[0146][0147][0148]
其中
[0149]
(3)如果收敛误差大于平移操作误差阈值ε
k
,即err>ε
k
,按下式向将单纯形所有顶点向最优点进行平移;否则,执行步骤4);
[0150][0151]
其中,λ∈(0,0.2)是平移系数。
[0152]
(4)进行单纯形反射、收缩、扩张、减小棱长操作计算;
[0153]
(5)如果收敛误差err大于迭代精度ε
e
,k=k+1,返回(2);否则,满足精度要求,迭代计算结束。
[0154]
如上所述,在迭代过程中连续不断地向最优点移动单纯形,而且单纯形彼此相似,保证了单纯形不退化、不畸形,增强了该算法的收敛性。
[0155]
3、基于改进的单纯形法的加速过程寻优控制
[0156]
在保证变循环发动机安全工作的前提下,采用改进的单纯形法对某型涡扇变循环发动机进行加速过程寻优控制,在保证变循环发动机安全工作的前提下,改进的单纯形法可以有效地缩短加速时间,达到寻优的目的。
[0157]
变循环发动机的加速时间定义为
[0158][0159]
式中:i为转子的转动惯量;n
max
为加速过程结束时的转速;n
idle
为慢车时的转速;δn
ac
为加速过程中涡轮的剩余功率。
[0160]
从上式可以看出:模式选择活门msv打开程度msv选定后,决定加速时间的因素主要是加速过程中的涡轮剩余功率δn
ac
。而涡轮的剩余功率主要决定于高压转子转速n
h
和高压涡轮前总温t
t4
。为缩短加速时间,就必须增大涡轮的剩余功率,也就是必须增大变循环发动机的高压转子转速和提高燃烧室后的温度。因此,本发明选取高压转子转速n
h
和高压涡轮前总温t
t4
作为加速过程寻优控制的目标函数。目标函数的数学表述如下:
[0161][0162]
上式中,n
hd
为高压转子的目标转速,n
h
为高压转子的实际转速。t
t4d
为高压涡轮前的目标总温,t
t4
为高压涡轮前实际总温。
[0163]
要保证变循环发动机在加速过程中稳定工作,本发明考虑的约束条件有:涡轮前温度不超温、高压压气机不喘振、高压转子不超转、风扇不超转、燃烧室不富油熄火、主燃烧室供油量不超过其最大供油量等等。
[0164]
考虑到目标函数、约束条件以及控制变量的影响后,需要寻找一组合适的msv,w
f
,a9,dvgl,dvgh,使变循环发动机加速时间最短,即需要求解如下非线性约束问题:
[0165][0166]
其中控制变量x=[msv,w
f
,a9,dvgl,dvgh]
t
,以上各个变量均在相应的变化范围之内取初值。
[0167]
变循环发动机的加速过程是一个动态过程,需要得到的优化结果是随时间变化的控制变量的轨迹曲线,但改进的单纯形法只适用于静态问题,要达到求解动态问题的目的,必须对目标函数、控制变量及约束条件进行适当的处理。由上式可知,本发明采用的是多目标最优控制方法,采用线性加权法将多目标函数转化为单目标函数,来确定寻优目标函数。

[0168][0169]
对上式进行离散化和归一化处理。这样处理的目的是为了消除目标函数中各参数量纲和量值变化范围的不同对优化结果的影响。最终的寻优目标函数可以写成以下形式:
[0170][0171]
上式中,ω
a
和ω
b
为相应目标函数的权重系数,满足ω
a
≥0,ω
b
≥0,其大小反映相应的寻优目标函数在多目标优化问题中的重要程度。
[0172]
参照目标函数的形式,对变循环发动机约束条件也进行离散化和归一化处理:
[0173][0174][0175][0176][0177]
以上g
i
(x)(i=1,2,...,11)构成约束函数矩阵g(x),考虑约束条件后,目标函数可化为:
[0178][0179]
其中ω=[ω1,ω2,ω3,ω4,ω5,ω6,ω7,ω8,ω9,ω
10

11
]为约束函数的权重调整系数矩阵,其中ω1,ω2,ω3,ω4,ω5,ω6,ω7,ω8,ω9,ω
10

11
为对应约束条件可调整权重系数,ω
·
g(x)的设计用于满足变循环发动机的约束条件。
[0180]
尽管上面已经示出和描述了本发明的实施例,可以理解的是,上述实施例是示例性的,不能理解为对本发明的限制,本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。
当前第1页1 2 3 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1