基于仿真的结构优化方法(txt+pdf+epub+mobi电子书下载)


发布时间:2020-08-17 17:08:01

点击下载

作者:毛虎平

出版社:电子工业出版社

格式: AZW3, DOCX, EPUB, MOBI, PDF, TXT

基于仿真的结构优化方法

基于仿真的结构优化方法试读:

版权信息书名:基于仿真的结构优化方法作者:毛虎平排版:小暑暑出版社:电子工业出版社出版时间:2019-05-01ISBN:9787121365720本书由电子工业出版社有限公司授权北京当当科文电子商务有限公司制作与发行。— · 版权所有 侵权必究 · —前 言 Foreword

人类行为和自然进化均与最优化紧密相关。工程师通过调整各种参数使得产品性能达到最优;制造商则设计不同的加工过程使生产率最高、生产成本最低或产品性能最好;自然系统的进化是优胜劣汰的过程,可使各种生物在最恶劣的条件下最适宜地生存;物理系统会自然地趋于能量最低的状态。

产品建模与仿真优化的最终目的是实现产品的优化设计。仿真优化是指基于系统仿真的参数优化,它针对仿真模型建立优化问题,采用相关的优化搜索算法进行求解,是基于仿真的目标或/和约束的优化问题。基于仿真的优化根据仿真问题的类型分为静态优化和动态优化。基于仿真的动态优化又可分为基于仿真的动态响应优化、基于仿真的动态特性优化和基于仿真的动态疲劳优化。基于仿真的动态优化不仅是仿真优化问题,而且还是动态优化问题。

机械结构几乎都是在动态载荷环境下工作的,其各种性能都是依赖于时间的函数。为了提升机器的动态性能,动态优化非常必要。然而,当前动态优化设计理论与方法难以适应现代化产品设计的需要。传统的机械结构优化几乎都是静态优化,没有考虑由于动态载荷的作用而产生的动态效应,即在静态力作用下,应用经典优化算法优化,难以实现动态载荷作用下机器性能最佳。另外,即便进行动态优化,也只是直接进行优化。由于动态分析的复杂性和高耗时性,直接进行动态优化收敛速度极慢,甚至发散。鉴于此,中北大学于2007年将作者派往华中科技大学攻读博士学位,师从动态仿真优化专家陈立平教授和吴义忠教授,重点研究结构在动态载荷作用下的优化理论与方法。本书系统总结了作者自2007年以来的研究成果。

本书章节安排如下:第1章绪论,主要介绍本书所做研究的目的、意义,以及国内外研究概况;第2章基于谱元法的结构动态分析方法,主要介绍任意载荷振动问题分析的Chebyshev谱元法、聚集单元谱元法在承受冲击载荷结构动态分析中的应用和非线性振动分析的Chebyshev谱元法;第3章基于时间谱元法的动态响应优化方法,主要介绍机械结构动态响应优化模型、动态响应优化方法、线性单自由度系统的最优化设计、线性两自由度减振器最优化设计和汽车悬架系统动态响应优化设计;第4章基于面向所有节点等效静态载荷的模态叠加法的结构动态响应优化方法,主要介绍模态叠加法、等效静态载荷法和关键时间点集;第5章基于局部特征子结构法的连续结构优化方法,主要介绍连续结构优化问题描述和子结构法;第6章基于子结构平均单元能量的结构动态特性优化方法,主要介绍结构动态特性优化问题描述、结构平均单元能量、子结构划分和基于子结构平均单元能量的结构动态特性优化方法的实施;第7章基于节点里兹势能主自由度的结构动态缩减方法,主要介绍节点里兹势能计算及主自由度选择、构造缩减系统和基于节点里兹势能主自由度的结构动态缩减方法的实施。

本书涉及的研究工作得到了国家自然科学基金项目“基于等效体积应变静态载荷的解空间谱元离散关键点方法及其在结构动态响应优化中的应用”(资助号:51275489)和山西省自然科学基金项目“基于能量评价PDOFs的多子结构精简等效静态载荷和刚体模态分离解析梯度的复杂结构动态响应优化基础研究”(资助号:201701D121082)的资助,在此表示诚挚的谢意。同时,中北大学机电工程学院的郭保全副教授,中北大学能源动力工程学院的张艳岗副教授、王强副教授、王艳华副教授、王军博士、刘勇博士、赵利华博士、王英博士、郑利峰博士、高鹏飞硕士等为本书的研究提出了许多宝贵建议,硕士生田力、刘鑫等在本书的写作修改方面做了大量工作,在此一并表示衷心感谢!

由于结构优化是一个具有挑战性的研究领域,并且仍处于蓬勃发展阶段,加之作者水平有限,本书难免存在疏漏之处,敬请读者批评指正。第1章绪论产品建模与仿真的最终目的是实现产品的优化设计。仿真优化是指基于系统仿真的参数优化,它是针对仿真模型建立的优化问题。采用相关的优化搜索算法进行求解的一整套技术,是基于仿真的目标或和约束的优化问题。其原理如图1.1所示,即基于模型仿真给出的输入关系构造优化模型,通过输出给优化算法得到最佳的输入量。如制造系统、交通系统、电力系统、化工系统[1-3]等的许多工程问题都能归结为基于仿真的优化问题。基于仿真的优化根据仿真问题的类型分为静态优化和动态优化。基于仿真的动态优化又可分为基于仿真的动态响应优化、基于仿真的动态特性优化和基于仿真的动态疲劳优化。基于仿真的动态优化不仅是仿真优化问题,而且是动态优化问题。从狭义上讲,动态响应优化是在动态载荷作用下的结构动态特性的最优化;从广义上讲,动态响应优化是指目标函数优化或约束函数与时间相关、而设计变量与时间无关的优化,当然后者包含前者。图1.1 基于仿真模型优化的原理

基于仿真的优化的特点如下:(1)系统输入输出的关系有两种,一种是缺少结构信息,不存在解析表达式,仅能通过仿真得到;另一种是虽然存在解析表达式,但是获得解析解比较困难,如微分方程或偏微分方程等,或者采用近似方法获得解析表达式,如结构的运动模型。(2)在仿真优化时,一种仿真耗费时间较少,另一种仿真耗费时间则较多,且缺少针对这种情况的仿真优化算法,导致优化过程十分耗时,甚至不可能实现优化。(3)对于大型连续结构,优化难度非常大,主要表现在:大型连续结构有限元分析耗时较长,寻优过程效率较低;连续结构参数化比较困难,但对于杆单元、平面单元、梁单元等简单单元来说,常以单元截面积、单元厚度及梁截面积为设计变量,很容易实现参数化。根据子结构方法的优势与优化过程中各个子功能的特点重新思考优化是一种有效的解决方案。(4)设计变量及其取值范围的确定是对结构优化问题建立合理优化模型的前提,它们的取值直接决定了优化问题的收敛性与高效性。工程师一般结合经验尽量多取设计变量,尽量扩大设计变量的取值范围,这样做是为了避免漏掉关键设计变量及其包含的重要取值范围,但这样确定设计变量的取值范围没有任何理论依据。

鉴于基于仿真的结构优化工程背景和上述特点,它一直是很多领域,尤其是机械制造和航空航天等领域学者和工程师们共同关注的重要课题。随着计算机技术、人工智能技术和数学分析方法的发展,对基于仿真的结构优化研究更加迫在眉睫。1.1 本书研究目的和意义

本书研究是在确定性优化的基础上进行的,其目的如下:(1)在谱元法的基础上,实现对任意载荷振动问题的Chebyshev时间谱元法分析、对承受冲击载荷结构的聚集单元谱元法动态分析和对非线性振动的Chebyshev时间谱元法分析。(2)研究基于时间谱元法的结构动态响应优化,并将谱元法与模态叠加法相结合应用于结构动态响应优化。

本书研究基于时间谱元法的系统动态响应设计,深入探讨时间域内的离散动态响应,将运动微分方程组转化为代数方程组,精确解出瞬态响应,用Guass-Lobbato-Legendre(GLL)点法和关键点法处理时间约束;针对瞬态动力学分析的复杂性和等效静态载荷转化的不确定性,提出基于模态叠加的所有节点等效静态载荷法,并将其应用于动态响应优化。从模态叠加原理出发,分析动态响应与各模态的关系,然后通过详细分析等效静态载荷法的原理,给出基于模态响应的所有节点等效静态载荷的计算表达式,最后提出关键时间点集的所有节点等效静态载荷法,采用谱元离散插值且微分获得了时间关键点,并与邻近的GLL点组成关键时间点集。(3)基于局部特征子结构方法的连续结构优化和基于子结构平均单元能量的结构动态特性优化。

本书从优化过程的各个子功能和连续结构的几何特征出发,将连续结构划分为参数化子结构、超单元和状态变量子结构,以参数化子结构的几何特征为设计变量,建立连续结构评价的目标函数;以最小的连续结构质量为优化目标,将参数化子结构和状态变量子结构所承载的连续结构应力应变作为约束条件,建立连续结构评价的最优化数学模型(其中隐含着整体结构局部几何不变而且不包含所需状态变量的局部构造,即超单元)。本书采用基于梯度的序列二次规划法对模型求解,实现一种基于局部特征子结构的连续结构高效优化方法。从建立优化模型出发,结合结构几何特性,将整体结构划分为多个准设计变量子结构,对于桁架结构来说,将每根杆作为一个子结构;对于连续结构而言,将多个单元作为一个子结构,再结合结构平均单元体积应变能及平均单元动能与结构动态响应贡献的关系,将平均单元体积应变能较大的结构单元作为尺寸应该变大的子结构,将平均单元动能较大的结构单元作为单元尺寸应该变小的子结构,从而确定设计变量的合理范围。本书推导了平均单元能量与结构动态响应的关系,构造了结构动态特性优化模型,调用基于梯度的优化器进行迭代寻优,实现了基于子结构平均单元能量的结构动态特性优化。1.2 国内外研究概况

对基于仿真的结构优化方法的研究涉及结构优化和基于仿真的优化两个方面。1.2.1 结构优化研究概况

结构优化的研究内容涉及子结构法、结构灵敏度计算、等效静态载荷法、谱元法及结构动态响应优化。

1)子结构法

结构优化设计是建立在精确的结构数学模型基础上的,而对于大型复杂结构,其动态响应还是以试验结果为准。对振动试验计算机仿真的研究难度很大,其核心问题是如何得到准确的振动响应。振动响应计算比静力计算、振动特性计算等难度都大得多。当前迫切需要提高大型复杂结构动态响应的计算精度与计算效率,子结构法是最合适的选择。

子结构法充分利用各个子系统的动态特性,以简便的计算过程获得可靠的原系统动力特性参数或动态响应,它是解决大型复杂结构动力分析的有效计算方法。根据采用的边界条件不同,子结构法可分为[4][5][6]四大类:约束子结构法、自由子结构法、混合子结构法[7][8][9]和载荷子结构法。在20世纪60年代初,Hurty首先提出模态综合法的概念,并建立约束模态综合法的基本框架,半个多世纪以来,模态综合法取得了很大的进展,其方法已定性化并已成为结构动[5]态分析的一种常规方法。2000年,Craig在Hurty的基础上进行了改进,形成约束子结构法。固定界面模态综合法是指在子结构的全部界面上附加约束,该方法在子结构的对接界面坐标中,不再区分静定约束坐标和赘余约束坐标,这使得固定界面模态综合法得到了广泛应用。自由界面子结构模态综合法是指以子结构的自由界面保留主模态集为其中一个子集,切断界面上的连接,将整体系统划分为空间中毫无关联的若干个子结构,然后利用相邻部件之间的界面位移协调及界面力的平衡条件,将之前完全释放的连接界面连接成一个整体。由于自由界面子结构模态综合法比较符合当前动态测试水平的要求,当需要通过试验来校核解析模型的可靠性时,这种方法很有吸引力。自由界面子结构模态综合法不仅计算精度优于固定界面子结构模态综合法,而且不受子结构界面坐标数目的限制,但是该方法复杂。混合界面子结构模态综合法是为了克服上述两种方法的缺点提出的。载荷子结构法改善了低阶模态的精度,但是随着模态数量的增加,其精度也将下降。以上述方法为基础还发展了很多其他方法,并得到了广泛应[10,11]用。

子结构法是一种通过对求解模型分块处理来降低仿真计算求解时间的高级有限元分析方法,只要分块合理并建立模型的子结构数据库,便能极大地提升求解效率。[12]

张焰等人将子结构法用于对汽车车架的降噪分析,求解效[13]率提升了94%。高鹏飞等人应用ANSYS前处理模块对活塞内冷油腔进行分离,使用子结构法中超级单元的思路对活塞的机械应力及温度分布进行分析,通过与传统计算方案对比,验证了子结构法在活塞结构强度校核及温度分布方面应用的高效性及准确性。高鹏飞等人[14]为了提升活塞结构优化设计的求解效率,首次将子结构法引入结构的优化设计中,以某增压柴油机活塞为研究对象,研究结果显示,其优化迭代收敛时间减少了74.14%。在实际问题中,常常需要对模型的某些局部特征进行优化设计,为了克服传统意义上建模和优化效[15]率低、参数化信息易丢失等缺点,刘波等人充分运用三维CAD软件和CAE软件的优点,以Pro/ENGINEER软件和ANSYS自身的脚本语言APDL为平台进行模型参数化建立,并提出了三种面向局部特征[16]优化的参数化设计方法。李芸等人将子结构法运用到大底盘双塔连体结构中,通过对其进行结构静力学分析及模态分析,充分展现了子结构法在大型复杂模型仿真计算应用中的优势。柴国栋等人[17]应用ANSYS中的子结构模块对电子设备箱进行模态分析,并寻[18]求提升整体结构刚度的方法。张明明等人将子结构法运用到柴油机曲轴的有限元仿真分析中,通过对曲轴进行子结构划分及子结构数据库建立,组合出柴油机曲轴的子结构有限元模型,对其进行结构静力学分析及模态分析,并与传统有限元算法进行对比,发现两种计算方法的误差不超过工程上的误差允许范围,再一次证明了子结构法[19]在结构有限元仿真计算领域的可行性。丁阳等人将子结构法引入钢框架结构抗倒塌性能的评估中,通过对两个五层钢框架结构的抗[20]倒塌评估可以证明该思路的高效性及准确性。丁晓红等人在对汽车座椅骨架进行拓扑优化设计时引入了子结构的思想,通过逐步逼近的方法得到了座椅骨架的最优结构,缩小了该座椅骨架的体积。张[21]盛等人通过对比多重多级子结构法与模态综合法在结构模态分析计算中的精度,证明多重多级子结构法更加高效准确,在结构的高[22]阶频率计算方面表现良好。张帆等人在对客车车身进行拓扑优化的过程中引入子结构法,将不参与优化的部分设计成一个子结构,并通过适当的方式与待优化部分进行节点连接,缩减了整体计算模型[23]的矩阵阶数,提升了车身拓扑优化设计的效率。李志刚等人将高架铁路浮桥进行子结构划分,重组计算后也证明了子结构法的高效性。笔者对模型采用基于梯度的序列二次规划法进行求解,以某柴油机活塞连续结构为例进行分析优化,并从收敛性和高效性方面与传统[24]优化方法进行比较,证明该方法的优越性。

另外,笔者研究了多领域仿真优化中SQP算法的并行处理与调度策略,提出了基于多领域仿真的SQP算法并行优化问题中的抽象调度模型,即等式约束离散变量优化模型,并对算法理论的可行性做了深入探讨;采用机群系统构建了并行仿真优化环境,在自主研发的多领[25]域统一建模与仿真平台MWorks下实现并行优化模块。

大型复杂结构动力学分析需要对拥有大量自由度的模型进行计算,在高频激励力作用下,要求计算步长非常小,导致计算耗时指数级增加。为了提高计算效率,可在保证一定精度的情况下,用少自由度模型代替多自由度模型,即模型缩减。所谓模型缩减,是通过一定的变换,将对总体结构动力学分析影响较小的次自由度,用对总体结构动力学分析影响较大的少量自由度表示,以达到减小计算规模的目的,其中少量自由度就是PDOFs。然而如何从庞大的自由度中选择PDOFs,目前在结构动力学领域仍属极具挑战的问题。

不过目前学术界提出了一些选择PDOFs的原则,最具代表性的有:①将结构振动方向定为PDOFs;②在质量或转动惯量相对较大而刚度相对较小的位置选择PDOFs;③在施加力或非零位移的位置选择PDOFs,这些原则仅仅是指导思想,具体选择PDOFs时,随意性较大。PDOFs的位置和数目直接影响模态分析的精度。Jeong等人[26]提出了一种阻尼系统基于自由度能量分布比率的PDOFs选择方法,为了估计结构的能量分布,采用双边Lanczos算法得到里兹向量,利用获得的里兹向量计算能量分布矩阵,将DOFs对应的最低的瑞利[27]商作为PDOFs。Kim等人提出了一种用于特征问题缩减的自由度分析选择方法,该方法根据结构系统模态自由度相关的能量进行选择,将能量分布矩阵加权行的值作为选择PDOFs的有效准则。Cho等[28]人提出了一种单元级的能量估计方法,构建了一个小规模的有限元模型,通过里兹向量来计算每个单元的能量,将该能量值排序,并把小的能量值作为PDOFs。

2)结构灵敏度计算

在结构优化中,灵敏度计算需要非常多的资源,因此有很多关于灵敏度计算效率的研究。灵敏度计算有三种方法,即有限差分法、基于离散方程的分析法和基于连续方程的分析法。其中,基于离散方程的分析法又分为分析法和半分析法,而基于连续方程的分析法基本上是完全分析法。分析法又包括直接差分法和伴随变量法。

在有限差分法中,中心差分法是最常用的一种,其中,目标函数和约束函数的灵敏度可表示为

式中,b为设计变量向量,ξ为特征值,z为节点位移向量。

当然,除了中心差分法,还有向前差分法和向后差分法,但中心差分法精度最高。有限差分法处理简单,而且可以采用现有仿真软件,通过将仿真模型看作一个黑箱函数来获得。但这种方法非常耗资源,特别是当需要反复进行有限元分析时。

当仿真代码计算耗费时间很少时,有限差分法是最好的灵敏度计算方法。然而,对于工程问题,仿真一般都采用商业有限元软件,此时分析法更适用。基于离散方程的分析法可以表示为

在式(1.2)中,最难计算的是,尤其是的计算很耗时。计算有两种方法,即直接差分法和伴随变量法。后者引入一个伴随方程,此时有

在优化过程中,伴随变量法需要求解的次数和被激活的约束个数相等。而直接差分法的求解次数和设计变量的个数相等。因此,评判两者的有效性需要具体问题具体分析。

在分析法中,质量矩阵和刚度矩阵的差分计算是最困难的。商业软件中对质量矩阵和刚度矩阵的计算通常采用基于有限元的有限差分法,然而其精度依赖于扰动尺寸的大小,特别是在大型复杂结构形状优化中,更是如此。基于连续方程的分析法是从积分公式开始,表示为

式中,,计算δz至关重要。基于连续方程的分析法同样也分为直接差分法和伴随变量法两种。

由于半解析灵敏度分析(SAM)法结合了分析法的精度和有限差分法的高效性,而且适合在商业软件中应用,因此该方法一直是研[29]究热点。1973年,Zienkiewicz和Campbell提出了半解析灵敏度[30][31]分析法。之后Barthelemy等人和Pauli在研究中发现,半解析灵敏度分析法在某些应用中出现了不精确的现象。为了解决这个问[32]题,Olhoff等人提出中间差分方案来求解刚度矩阵的微分。1993[33]年,Cheng和Olhoff研究发现了半解析灵敏度分析法存在不精确现象的真正原因,即当单元中存在刚体运动时,就会表现出不可靠的精度,其本质是刚体运动与截断误差直接相关。从这个原因出发,[34]Keulen和Boer提出了精细的半解析灵敏度分析(RSAM)法,它基于精确的刚体模态差分来消除由刚体模态导致的灵敏度误差。然而,当扰动尺寸比较大时,RSAM法也不能获得足够的精度。因此,在半解析法中,需要考虑高阶项,在高阶项的扩展中逆矩阵可以通过[35]Neumann级数展开,同时对基于模态分解的RSAM也进行了一定[36]研究,并描述了其在非线性结构分析中的应用。Cho和Kim结合模态分解及Neumann级数展开研究了RSAM法。

3)等效静态载荷法

等效静态载荷法是将在等效静态载荷作用下的结构位移场与在动[37]态载荷作用下某一时刻的位移场等效,研究人员在前期研究中将等效静态载荷的概念进行了扩展:等效静态载荷不仅要代替动态载荷作用下产生的位移场,而且要代替体积应变能,也就是说其代替效果包含位移场和体积应变能。在等效载荷思想的驱动下,演化出两类方法:①基于时间关键点的等效静态载荷的结构动态响应优化;②基于所有结构动态分析时间步或设计者指定时间细分点结构等效静态载荷的结构动态响应优化。后者考虑了所有可能,将结构动态分析的每一步或设计者指定时间细分点的每一点都等效为一组静态载荷。对于过小的载荷步来说,该方法耗时太长,而对于稍大的载荷步来说,结构动态响应分析的相对精确性会受到一定影响。当然,如果有足够好的计算环境,这个问题影响不大。如果采用设计者指定的时间细分点,虽然不存在精度问题,但存在如何细分等问题。无论哪种方式,当考虑所有的时间点时,都存在庞大的约束条件和载荷状态,这给优化算[38]法带来了极大挑战。笔者研究了如何高效识别关键时间点,在关键点时刻,通过体积应变能等效将动态载荷更加合理地转化为静态[39]载荷,然后应用SQP多初始点方法求解,使其较好地收敛。然而,在该研究中没有考虑载荷位置矢量,即没有考虑等效静态载荷作用位置,只是将等效静态载荷作用在动态载荷作用的位置或凭经验应该加载荷的位置,这无疑存在一定的随机性。并且等效静态载荷作用的位置不同,再加上其取值空间的不确定性,所需计算时间与获得的结果自然不同。而且求解等效静态载荷的本质是一个优化问题,这样以等效静态载荷作为设计变量的优化需要确定其取值空间,这又存在一定的随机性,这些随机性会导致结果具有不确定性,而且求解等效静态载荷很耗时。针对结构动态分析的复杂性和等效静态载荷转化的不确定性,笔者提出基于模态叠加的所有节点等效静态载荷法,并将其应用到动态响应优化中,对124杆桁架结构进行动态响应尺寸优化和对18杆桁架结构进行尺寸与形状混合优化设计的结果表明,该方[40]法是可行的和有效的。针对结构动态响应优化中动态分析的复杂性与高耗时问题,笔者提出了基于全局动态应力解空间谱单元插值的关键时间点识别方法,找到了结构动态响应下最危险的时刻。具体来说,首先,利用模态叠加法,获得结构的模态应力分布,并计算全局动态应力解空间;然后,利用谱单元离散动态应力绝对极大值点曲线,采用Lagrange插值并调用区域细分全局优化求解器,找到全局动[41]态应力的极大值与极小值,即得到关键时间点。

4)谱元法

谱元法(Spectral Element Method,SEM)基于弹性力学方程弱形式基础,在有限单元上进行谱展开,该方法具备有限元方法适应任意复杂介质模型的韧性和谱方法的精度,又称为高阶有限元方法或谱方法的域分解。[42]

Patera于1984年提出谱元法,并应用于流体动力学,其将有限元方法的处理边界和结构的灵活性与谱方法的快速收敛性结合起来。在相同精度的情况下,谱元法采用了较少的单元,减小了计算开销。谱元法包括空间谱元法、时间谱元法和空间-时间谱元法。空间谱元法利用区域嵌入技术,将实际问题中复杂的几何区域嵌入规则的矩形区域中,构造适当的谱元空间(相当于有限元方法的有限元空[43]间),解决了谱方法对区域的要求。时间谱元法则在有限元空间的基础上构造谱时间单元,然后在每个单元内进行插值,最后求解线性方程组。空间-时间谱元法是在一个谱单元内,将空间或时间离散为与GLL多项式零点或Chebyshev多项式零点相对应的网格点,在这[44]些点上进行Lagrange插值。从理论上来说,在一定点数上插值,当这些点是对应的正交多项式的零点时,获得的插值精度最高[45]。

关于谱元法的许多工作都取得了一定的进展。谱元法广泛应用于[46][47]可压流体和不可压流体的数值模拟。Pathria采用谱元法解[48]决了非光滑域的椭圆问题。Hesthaven提出了使用开边界条件区[44]域分解的谱方法。M.H.Kurdi将时间谱元法用于常微分方程的整[45]体求解,笔者在M.H.Kurdi工作的基础上,将时间谱元法用于结[49]构动态响应仿真。波的传播为了扩展谱元法的适应性,针对承受冲击载荷的结构动态响应问题,从谱单元离散方案出发并根据冲击载荷的特点,以冲击载荷最大值点为中心将谱单元尺寸按一定比例等比例向两侧扩大,实现单元尺寸与载荷特征相适应。在此基础上,将动力学方程转化为一阶线性微分方程组,通过Bubnov-Galerkin法获得离散线性方程组并采用高斯消元法求解。将其与等距谱元法进行比[50]较,可证明该方法的可行性和有效性。笔者研究了用Chebyshev时间谱元法求解任意载荷作用下的振动问题,从Bubnov-Galerkin法出发,在第二类Chebyshev正交多项式极点处重心Lagrange插值来构造节点基函数并分析其特性,推导了任意载荷作用下振动问题的Galerkin谱元离散方案,利用最小二乘法求解线性方程组;以线性载荷、三角载荷、半正弦波载荷作用下的振动问题及正弦载荷作用下的悬臂梁振动问题为例,验证了该方法的可行性,并与配点法进行比较,[51]进一步说明了该方法的高精度性和可靠性。笔者还研究了用Chebyshev时间谱元法求解非线性振动问题,从Bubnov-Galerkin法出发,在第二类Chebyshev正交多项式极点处重心Lagrange插值来构造节点基函数并分析其特性,推导了非线性振动问题的Galerkin谱元离散方案,利用Newton-Raphson法求解非线性方程组;对于非线性单摆,还需要将二分法和重心Lagrange插值结合来求解角频率;以Duffing型非线性振动和非线性单摆振动问题为例,表明了该方法的可[52]行性和高精度。针对结构动态响应方程自由度很大,而谱元法是大自由度的矩阵与整体时间矩阵张量的乘积,求解起来非常耗内存和时间的问题,笔者提出了逐步时间谱元法。将仿真时间划分为很小的时间段,在每个时间段内划分单元,并在每个单元中采用谱展开近似,这样处理具有有限元处理复杂结构及边界的灵活性、谱方法的高[53]精度及快速收敛、逐步划分仿真时间的高效性等特点。针对传统结构动态响应优化方法的不足,笔者提出了元模型混合自适应优化与时间谱元法相结合的方法,从Bubnov-Galerkin法出发,深入探讨时间域内的离散动态响应,将整体结构动力学方程转化成代数方程组,精确、高效地求解动态响应;依据优化问题的特点,采用自适应策略选择对应的元模型进行优化;在优化过程中利用均匀网格获取有潜力点的数量,并将局部优化与多元模型混合自适应方法融合,使得优化结果更可靠。为了处理与时间相关的约束,笔者提出将关键点及其相邻Gauss-Lobatto-Legendre点组成集合的关键点集方法。但不论是精确性还是效率,元模型混合自适应优化与时间谱元法相结合的方[54]法都更优于关键点集方法。

谱元法的精度既可以通过增加每个单元上谱方法的自由度实现,也可以通过增加单元的数目来实现,最好的情形是每个单元上的自由度可以自由调节而不会相互制约,这样的谱元法才具有足够的灵活性。

一般地,采用时间谱元法求解结构动力学方程是不合适的,原因在于:如果离散的单元过多,则解线性方程组的过程费时并影响其应用;如果离散的单元过少,则动态响应不够准确,难以满足工程要求。另外,由于结构有限元离散后,自由度一般非常大,而时间谱元法求解结构动力学方程需要矩阵求逆运算,因此工程应用困难。解决此问题的方法是采用时间分段求解,并且在每一段中采用逐元技术[55],然而,这样也不能从根本上解决其工程应用问题。在工程结构动力学分析中,结构单元数比较多,逐元技术每一单元的矩阵求逆运算也阻碍了它的工程应用。

5)结构动态响应优化[56]

20世纪60年代,Niordson提出了结构动态特性优化的概念并进行了相应的研究,从此拉开了结构动态优化——结构动态特性优化的序幕。早期的结构动态特性优化方法是分布参数结构优化方法,属于解析法,由于该方法中偏微分方程求解困难,因此只适用于一些简单的结构,对于大型复杂结构无能为力。随后,准则法和数学规划法得以发展,目前这部分内容已经比较成熟。

在结构设计中,精确获得外部载荷非常关键,但是在多数状况下困难重重,因此,一般先设置某个静态载荷进行静态优化设计。从严格意义上来说,结构所受载荷都是动态的。动态因子法可以将动态优化问题转化为静态优化问题,可是这样处理往往会造成结构过设计或[57]欠设计,所以结构动态优化直接采用动态载荷更为合理。[58]

Wang等人应用数学规划方法对动态载荷作用下的平面正交钢框架结构进行了优化设计,在优化时以结构固有频率不小于一定值、最大动位移和动应力不大于一定值为约束条件,以结构总质量为优化目标,但没有考虑结构阻尼。另外,其研究发现,结构参数的可行域在对其进行动态响应优化设计时一般是不连续的。秦健健等人[59]针对某柴油机连杆质量过大的问题,采用基于ANSYS的有限元结构仿真分析方法,利用APDL语言建立了柴油机连杆的有限元模型,在有限元分析的基础上,利用ISIGHT集成优化软件结合多岛遗传算法对连杆杆身进行优化设计,从而使杆身的质量降低了6.02%。Lin等[60]人使用单元重构法和结构形状渐进算法,对同时具有静态约束和动态约束的结构尺寸和节点坐标进行了有效的优化设计。[61]Pantelides等人针对其在优化时初始设计点为非可行点,使用一般优化方法不一定能收敛于全局最优解的缺点,将MISA算法(改进的模拟退火法)应用于同时考虑结构动位移和动应力约束的结构动力优化问题,并将MISA算法与一般优化方法进行综合比较,验证了[62]MISA算法的优势。Min等人利用均匀化和直接积分方法对冲击载荷作用下的薄板结构进行拓扑优化设计,此项工作具有一定的前瞻性[63]和开创性。Du等人与以往考虑固有频率和动响应位移的动力优化不同,他们主要考虑如何降低结构的声辐射强度,在分析中忽略结构与声传播媒介的耦合作用,在简谐激励力的作用下,计算分析其结构参数的灵敏度,在此基础上成功进行了振动结构的拓扑优化。

概括起来,结构动态响应优化设计分为三个研究方向:①与时间相关约束的处理;②灵敏度分析;③近似。笔者所在课题组对①和③两个方向进行研究,提出基于GLL点集的处理与时间相关约束的方法[45,63,64]。该方法采用满足精度要求的较少谱单元求解运动微分方程,在每个单元内,对其高次Lagrange插值函数进行一维搜索找到单元绝对值极值点,将所有单元响应绝对最大值及其相邻的2个GLL点作为约束,将最大值附近的其他一些点包含其中,构成GLL点集约束。构件在动态载荷作用下产生的位移非常小或仅仅考虑某一方向的位移时,其几何形状或尺寸或多或少都会发生变化,这时构件内部每个单元体都会因动态载荷作用引起形状的相对改变,笔者所在课题组将动态载荷变化和体积应变对应起来,提出了等效体积应变静态载荷法[65],推导了以静态载荷作用的体积应变和动态载荷变化的函数关系,实现等效体积应变。1.2.2 基于仿真的优化研究概况

仿真优化迭代过程中需要调用仿真程序来计算目标函数和约束函数的值。响应面方法是提高仿真优化效率的有效途径。响应面12(Response Surface)是指输出响应变量Y与一组输入变量(x, x,n…, x)之间的函数关系。通常,响应面反映的是某个计算密集的复杂原模型(如多领域仿真模型、FEA模型、CFD模型等)的近似模型。因此,响应面又称为代理模型(Surrogate)或元模型(Meta-model),即模型的模型。

响应面方法(Response Surface Method,RSM)是指通过构造原模型的响应面来解决原模型的设计或分析等问题的近似方法。据报[66]道,福特汽车公司进行一次汽车碰撞模型的仿真分析需要36~160小时。要实现该模型两个变量的设计优化,假设平均需要50次迭代寻优,而每次迭代需要进行一次仿真计算,则获得该优化问题的解需要75天至11个月。同样,要实现FEA模型或CFD模型的设计优化可能需要更长的时间。这在实践中几乎是不可接受的。因此,在过去的20年中,响应面方法应运而生,而且得到了快速发展。该方法能够减少优化迭代过程中原模型的仿真次数,而且响应面都是基于采样点数据构造的,而采样点估值计算都是彼此独立的(传统的优化迭代过程是序列估值的),可以方便地通过并行计算来获得,因此,该方法可以大大提高复杂分析模型的设计优化效率。

根据文献[67, 68],RSM的作用包括以下4个方面:(1)模型近似:这是RSM的基本功能。建立一个复杂原模型在其全局定义域内的响应面近似模型,可以利用该近似模型实现新未知设计点的快速估值。(2)设计空间探索:基于建立的响应面模型可以帮助工程师或设计人员进行参数试验、灵敏度分析,以及响应变量与输入参数之间的函数关系可视化,从而帮助工程师更好地理解原模型的特性。(3)优化问题的准确表达:基于对响应面模型的设计空间探索,特别是灵敏度分析,可以帮助设计人员构造更加准确的设计优化问题。例如,可以从设计变量集合中剔除那些非敏感的参数,从而减少设计变量的维数。根据参数试验也可以缩减搜索区间,从而减少采样区间,进而减少优化迭代次数。同样地,通过分析,一个多目标设计优化的问题可能简化成单目标优化问题;而是一个简单的单目标设计优化问题,通过设计空间探索,可能需要建立多目标设计优化问题才能解决。(4)优化方法的支持:这是目前RSM的主要应用领域。利用建立的响应面模型可以辅助完成各种涉及原模型仿真的设计优化问题,如全局优化、多领域仿真优化、多目标优化、多学科设计优化及概率设计优化,包括可靠性优化、稳健优化等。

响应面方法中的一个重要环节是构造响应面模型。响应面模型种类较多,分别适合不同的需求,常用的有多项式回归(Polynomial Regression Surrogate,PRS,通常称为RSM模型)、Kriging插值模型、径向基函数(Radial Basis Functions,RBF)模型、支持向量回归(Support Vector Regression,SVR)模型、神经网络(Neural Network,NN)模型,以及基于RBF和NN混合的RBNN模型;还有基于样条的多元自适应回归样条(Multivariate Adaptive Regression Spline,MARS)模型、BMARS(B样条MARS)模型和NURBs模型,还有归纳学习(Inductive Learning)模型、最小插值多项式(Least Interpolating Polynomial,LIP)模型等。

要构造一个响应面模型,第一步是在设计空间内采样,形成设计点集S;再进行仿真计算(也称“昂贵”计算,Expensive Calculation)获得响应数据集Y;最后根据不同的算法由S和Y构造不同的响应面模型。试验设计方法提供了各种采样策略,主要包括两大类:边缘分布型和全空间分布型(Space Filling)。边缘分布型也称经典采样方法,利用该方法采样,其采样点主要分布在设计域的边界附近,典型的边缘分布型采样方法有:全因子/部分因子试验(Full/Fractional Factorial Design)、中心复合试验(Center Composite Design,CCD)、Box-Behnken等,还有Taguchi、D-Optimal、Plackett-Burman等方法。全空间分布型是指采样点布满整个设计域,该类型的采样方法有:简单网格、拉丁超立方设计(Latin Hypercube Design)、正交表(Orthogonal Array),还有随机采样、一致设计(Uniform Design)、混杂网络(Scrambled Nets)、蒙特卡罗仿真和Hammersley序列设计等。

一般地,通过一次采样构造一个响应面模型是不合适的,原因在于:如果采样点过多,则构造过程费时且可能影响使用;如果采样点过少,则构造的响应面模型不够准确,难以满足应用需求。另外,由于原模型的性态未知,难以确定合适的采样方法。解决此问题的方法是基于序列自适应采样来构造序列响应面模型。序列自适应采样的主要思想是根据近似值与真实值之间的误差来确定采样点的疏密。序列探索试验设计(Sequential Exploratory Experiment Design,SEED)[69]方法是此类试验设计方法的代表,i-Sight软件中使用模拟退火算法来进行自适应采样。

在使用响应面方法解决实际工程问题时,须综合考虑如下5个方面的因素:(1)响应面的精确度。毫无疑问,响应面模型的精确度是近似的基本要求。(2)原模型的仿真估值次数,即构造响应面模型所需的总采样点数目。由于每次估算的计算费用较高,需要限制原模型的仿真次数。在相同精度下,原模型仿真估值次数越少,则该响应面模型构造效率越高。(3)构建和优化响应面的时间。如前所述,响应面模型的构造往往是一个逐步精细的过程,序列自适应采样构造逐步精细的响应面模型是目前响应面方法的一个研究热点,特别是在采样点逐步增多的情况下,如何快速更新响应面模型以实现其增量构造算法是各种响应面方法值得探索的课题。(4)响应面模型占用的存储空间。显然,响应面模型本身所占的内存空间越大,则构造过程越慢,利用其进行估值也越慢。因此,在同等情况下,响应面模型本身所含的信息越少,也即所占的内存空间越小越好。(5)利用响应面模型对给定点估值的速度。建立响应面模型的最终目的是使用其进行估值,而且通常这个估值过程被称为“便宜计算”(Cheap Calculation),因此,在应用如优化迭代过程中被大量执行。可见,如果对给定点的估值速度太慢,就会使理想的“便宜计算”变得不“便宜”。1.3 本书主要研究内容、主要创新点与组织结构1.3.1 主要研究内容

本书针对基于仿真黑箱函数模型结构优化的工程实际问题,研究有效且高效的优化技术,主要研究内容如下:(1)本书研究用Chebyshev时间谱元法求解任意载荷作用下的振动问题,从Bubnov-Galerkin法出发,深入分析在第二类Chebyshev正交多项式极点处重心Lagrange插值构造的节点基函数及其特性,推导了任意载荷作用下振动问题的Galerkin谱元离散方案,利用最小二乘法求解线性方程组。为了扩展谱元法的适应性,针对承受冲击载荷的结构动态问题,从谱单元离散方案出发并根据冲击载荷的特点,以冲击载荷最大值点为中心将谱单元尺寸按一定比例等比向两侧扩大,实现单元尺寸与载荷特征相适应。在此基础上,将动力学方程转化为一阶线性微分方程组,通过Bubnov-Galerkin法获得离散线性方程组,应用高斯消元法求解。研究用Chebyshev时间谱元法求解非线性振动问题,从Bubnov-Galerkin法出发,深入分析在第二类Chebyshev正交多项式极点处重心Lagrange插值构造的节点基函数及其特性,推导了非线性振动问题的Galerkin谱元离散方案,利用Newton-Raphson法求解非线性方程组。(2)本书研究基于时间谱元法的系统动态响应设计。深入探讨时间域内的离散动态响应,将运动微分方程组转化成代数方程组,精确求解瞬态响应,用GLL点法和关键点法处理时间约束。以弹簧减振器设计为例,引入人工设计变量,详细分析两种处理约束方法的优缺点,也说明了此方法的正确性。这些内容可为进一步研究动态响应优化提供参考,如在此基础上研究复杂系统的灵敏度分析,以提高此方法的实用性等。(3)本书针对瞬态动力学分析复杂性和等效静态载荷转化的不确定性,提出基于模态叠加的所有节点等效静态载荷法,并将其应用于动态响应优化。首先从模态叠加的原理出发,分析了动态响应与各模态的关系;然后通过详细分析等效静态载荷法的原理,给出利用模态响应的所有节点等效静态载荷的计算表达式;最后提出关键时间点集的所有节点等效静态载荷法,采用谱元离散插值且微分获得了时间关键点,并与邻近的GLL点组成关键时间点集。(4)本书为了实现连续结构优化的可行性和高效性,本书提出一种基于局部特征子结构的优化方法。从优化过程的各个子功能和连续结构的几何特征分析出发,将连续结构划分为参数化子结构、超单元和状态变量子结构,以参数化子结构的几何特征为设计变量,建立连续结构评价的目标函数;以最小的连续结构质量为优化目标,将参数化子结构和状态变量子结构所承载的连续结构应力应变作为约束条件,建立连续结构评价的最优化数学模型,其中隐含着整体结构局部几何不变且不包含所需状态变量的局部构造。对于模型求解,采用基于梯度的序列二次规划法进行求解。以某柴油机活塞连续结构优化为例进行分析优化,并与传统优化方法从收敛性和高效性方面进行比较,说明了本书方法的合理性和优越性。(5)为了提高结构优化的可行性和高效性,本书提出了一种基于子结构平均单元能量的结构动态特性优化方法。从建立优化模型出发,结合结构几何特性将整体结构划分为多个准设计变量子结构。对于桁架结构来说,将每根杆作为一个子结构;对于连续结构而言,多个单元作为一个子结构,再结合结构平均单元体积应变能及平均单元动能与结构动态响应贡献的关系,将平均单元体积应变能较大的结构单元作为尺寸应该变大的子结构,平均单元动能较大的结构单元作为单元尺寸应该变小的子结构,从而确定设计变量的合理范围。本书推导了平均单元能量与结构动态响应的关系,构造了结构动态特性优化模型,调用基于梯度的优化器进行迭代寻优。(6)为了提高结构动态分析的效率,本书提出了一种基于节点里兹势能主自由度的结构动态缩减方法。阐述了改进缩减系统法的原理,分析了里兹向量提取过程,定义了节点里兹势能,并利用其作为捕捉能精确反映结构动态特性的主自由度的依据,给出了节点里兹势能的计算公式。通过计算分析圆柱曲板和曲轴两个实例,验证了本书方法的可行性和优越性。研究结果表明,在里兹向量空间定义节点里兹势能更容易捕捉高精度的动态特性,加权系数可提高高阶频率的精度,在结构缩减中,主自由度约取总自由度的1/3比较合适。1.3.2 主要创新点

本书研究的主要创新点如下:(1)将时间谱元法应用于机械动态问题仿真中,提出了任意载荷振动问题分析的Chebyshev谱元法、承受冲击载荷结构动态分析的聚集单元谱元法和非线性振动分析的Chebyshev谱元法;将时间谱元法应用于机械系统动态响应优化中。应用谱元法求解机械系统动态响应,改善了传统求解动态响应时误差大的缺点,达到谱收敛精度。这样,动态响应优化就可以在超曲线或超曲面上找到满足所有时间约束变化的目标函数。在处理约束上,采用了GLL点法和关键点法两种方法,并比较了这两种方法的优缺点。(2)将子结构法应用于连续结构优化。将优化三要素分别与连续结构的几何特征结合起来,分别定义了参数化子结构、超单元和状态变量子结构。设计变量对应参数化子结构的几何参数,目标函数和约束函数对应状态变量子结构的响应值,超单元是连续结构中既不包含设计变量也不包含状态变量的部分。这样定义不仅可以提高优化效率,而且能够减小有限元参数化的难度。(3)提出一种基于子结构平均单元能量的结构动态特性优化方法。对于一些特殊结构,如桁架结构,将结构划分为多个准设计变量子结构,分别以平均单元体积应变能和动能确定变大子结构和变小子结构,进而确定设计变量范围,为建立精确的优化模型奠定了基础。1.3.3 组织结构

各章组织结构如图1.2所示。图1.2 各章组织结构第2章基于谱元法的结构动态分析方法[70]振动问题的分析在瞬态分析、机械工程、故障诊断等领域具有不可替代[71]的地位,是机械设计中一种重要的分析。如何获得精确的振动响应一直[72]是研究者的工作重点。振动问题分析的数值方法不断发展,不管何种方法都各有利弊。其中,具有[73]代表性的是摄动传递矩阵法和有限差分法(FDM)。前者是考虑系统的随机性,后者是直接将微分方程(组)转化为代数方程组,其数学概念简单直观,表达式简单,容易编程,时间步长直接决定计算的收敛性和精确性[45][74-77]。对Chebyshev拟谱方法的研究表明,仿真时间范围较小,可以获得很高的精度,否则,获得的解没有意义。利用有限元法变分原理和差分方法的优点,可以进一步获得一定精度的近似解,并且当形函数为正交多项式零点或极点的插值基函数时,将其称为谱元法,正交多项式为Chebyshev正交多项式时,称其为Chebyshev谱元法。[42,78,79]Steven Orszag于1969年提出谱方法,通过大量研究表明其具有高阶数值分析快速收敛的优点,然而不能处理复杂空间域等缺点限制了其发展。考虑低阶有限元方法在非结构化域的灵活性和谱方法的高精度及谱收敛[42]的特点,Patera于1984年提出谱元法,采用在GLL的零点Lagrange插值[80]和p型节点基函数,并应用到流体动力学中。Dimitri Komatitsch将有限元方法的灵活性和谱方法的准确性结合,并将其引入三维地震波计算,利用高阶Lagrange插值对单元上的波场进行离散,然后根据高斯-洛巴-勒让德积分规则对单元进行积分。近年来,谱元法已经被应用于科学和工程的很多领域[44,81]。文献[82]采用勒让德四边形谱元近似Black-Scholes方程,并将其应用于欧洲彩虹和篮子期权的定价。文献[83]采用谱元法对球面几何中的[64]浅水方程进行分析,并将其与其他模型进行了比较。笔者应用Lobatto-Legendre正交多项式的零点将振动微分方程的时间域谱展开,通过Galerkin谱元离散方案获得精确解,并提出GLL点集法处理与时间相关的约束,进一步[84]提出逐步时间谱元法,缩短了CPU时间。P.Z.Bar-Yoseph 等人采用时间[85]谱元法求解非线性混沌动态系统。U.Zrahia和P.Z.Bar-Yoseph 采用空间-时间耦合谱元法求解二阶双曲方程。然而很少有文献使用Chebyshev谱元法分析振动问题。本章提出采用Chebyshev谱元法进行任意载荷作用下的振动问题分析,通过稳定性很好的重心Lagrange插值近似振动解函数,结合有限元方法的灵活性,得出振动问题的h收敛和p收敛两种收敛方式,并与配点方法进行比较。2.1 谱元法[42]

谱元法是Patera在1984年提出来的,应用于流体动力学,其将有限元法处理边界和结构的灵活性与谱方法的快算收敛性结合起来。在要求相同精度的情况下,谱元法能够采用较少的单元减小计算开销。在一个单元内,谱元法将时间离散为与GLL多项式零点相对应的网格点,在这些点上进行Lagrange插值。

在一定点数上插值,从理论上讲,当这些点是对应正交多项式的[86]零点时,获得的插值精度最高。如图2.1所示,均匀点分布和GLL点分布近似Runge函数[见式(2.1)]。在均匀分布的有限元法中,随着插值次数的增加,近似的数值误差极大或近似失败,此时GLL插值有明显的优势。GLL点在标准区间的分布如图2.2所示。图2.1 均匀点和GLL点插值近似Runge函数

微分方程包括动力学方程或方程组,常用来描述自然界的一些物理现象。差分法是直接基于这些微分方程或方程组进行离散化的。但在多数情况下,描述同一物理过程或现象,可用不同的形式。从物理意义上的守恒定律出发,可以推导出变分原理。而变分问题与微分方程定解问题在某种意义下等价。谱元法是基于变分原理的一种离散计算方法。下面讨论一阶线性微分方程的初值问题[见式(2.2)]。高阶线性微分方程可以化成一阶线性微分方程组。v

式中,x是依赖于时间的状态变量,且,N是状态变量的s个数;t是时间;f(x,t)是状态变量x和时间t的函数;A是状态变量的耦合矩阵,与时间无关。图2.2 GLL点在标准区间的分布el注:ratio表示最后一个单元与第一个单元的比值;N表示单元数;p表示单元内插值节点数。2.1.1 瞬态和稳态响应分析

应用谱元法将每个状态变量在给定时间段内离散化,并近似为m次Lagrange多项式:k

式中,为第j个单元的第k个m次Lagrange多项式;ξ为定义在[-1,1]上的GLL点;为第j个单元上未知节点在GLL点的值。Lobatto多项式是Legendre多项式微分定义的正交多项式,见文献[86]第146~151页。其中,ξ∈[-1,1]是经过式(2.4)映射获得的。

式中,分别为第j个单元的第一个节点和第二个节点;ξ为经过域变换后的变量,且ξ∈[-1,1]。这个域如图2.3所示。图2.3 时间域可离散为谱单元,每个单元近似为基于GLL点的m次Lagrange多项式

将式(2.3)代入式(2.2),在每个单元上应用Bubnov-Galerkin[87]法使得插值误差最小,得到(j)

式中,n=0,1,…,m ;h是第j个单元的长度;是m次Lagrange多项式的第n个插值的基函数。对式(2.5)每一部分得到

式(2.6)中的积分是用Gauss-Lobatto求积公式[见式(2.7)]求的。q

式中,I是ξ的一般函数;ω是Gauss-Lobatto积分在q点的权值[87]。每个单元可用矩阵形式表示为

其中,

由于相邻两个单元共享其中一个元素,所以应满足:[86]

通过式(2.8),利用连接矩阵C,将一个状态变量的所有谱单元组装起来,就得到了一个状态变量的Galerkin近似方程:uωu

式中,B、B是全局微分矩阵与全局权矩阵;F (X)是激励力的全局形式,其中,

试读结束[说明:试读内容隐藏了图片]

下载完整电子书


相关推荐

最新文章


© 2020 txtepub下载