有限元法在内燃机上的应用.ppt

上传人:牧羊曲112 文档编号:5755855 上传时间:2023-08-17 格式:PPT 页数:109 大小:6.35MB
返回 下载 相关 举报
有限元法在内燃机上的应用.ppt_第1页
第1页 / 共109页
有限元法在内燃机上的应用.ppt_第2页
第2页 / 共109页
有限元法在内燃机上的应用.ppt_第3页
第3页 / 共109页
有限元法在内燃机上的应用.ppt_第4页
第4页 / 共109页
有限元法在内燃机上的应用.ppt_第5页
第5页 / 共109页
点击查看更多>>
资源描述

《有限元法在内燃机上的应用.ppt》由会员分享,可在线阅读,更多相关《有限元法在内燃机上的应用.ppt(109页珍藏版)》请在三一办公上搜索。

1、有限元法在内燃机上的应用入门,前言有限元法简介,内燃机设计提出的问题,曲轴的扭振(曲轴的固有频率、振型、扭振的振幅)轴承的计算(轴心轨迹、油膜厚度、压力)曲轴强度(强度:应力)活塞(温度场:热应力、机械强度:应力)连杆(强度:应力)机体(刚度:变形量、强度:应力)缸盖(鼻梁区:热应力)配气机构动力学(零部件刚度)工具-有限元法(FEA),而不是 FEM 理论,有限元法简介,有限元法的概念,有限元法(有限单元法、有限元素法)Finite Element Method(FEM)Finite Element Analysis(FEA)无限元法:(发动机辐射噪声场)Infinite element m

2、ethod(IFEM)有限元法是可以求解工程、数学、物理方面具有初值、边值条件的微分方程的一种近似的数值计算方法。(数学),经典有限元法求解思路,数理方程、传热学中内部无热源的平面温度场方程区域A可以得到微分方程的解析解(转为齐次方程)区域B、C难以得到解析解,可得近似数值解。如有限差分法、有限元法、边界元法等,有限元法为最主要的方法之一。轴承用有限差分法。噪声用边界元法。流体用有限容积法,有限元解微分方程方法,整个区域中温度T的分布复杂,难以找到一个函数来描述T的分布,可把复杂的区域划分为较小的规则的区域(单元element),则在每个较小的区域温度T的分布简单,可用假定的简单函数来表示,可

3、以通过有关的数学变换(本课主讲内容)将求解微分方程问题变为求解线性方程组的问题(应用线性代数)。,有限元发展过程,工程方面1956年Boeing公司的Turner Clough分析飞机的结构(采用自然离散的方法)1960年Clough处理平面连续弹性体问题,提出Finite Element Mthod1956年-其它大学、公司,数学方面 1943年Courant研究分片连续和最小势能问题1963-1964年Besseling,Melosh,Jones研究FEM和Ritz法的关系和变分原理,有限元法实现方法,有限元法将函数定义在简单几何形状(标准件)的单元域上,且不考虑整个定义域的复杂边界条件,

4、是有限元法优于其他近似方法的原因之一离散的概念,曲轴的离散,曲轴的离散:二维(质量单元+梁单元)三维:四面体单元(棱锥)or 六面体单元(砖),哪些方面可用 FEA,固体力学(力学和土木)流体力学(可用FEA,但多用有限容积法)声场(可用FEA,但多用边界元法BEM)传热学电磁学,FEA在内燃机方面应用,刚度(变形)和强度(应力):曲轴、连杆、活塞、机体、缸盖、配气等所有机械零件振动和噪声:整机振动、曲轴扭振、整机噪声传热:活塞、缸盖、缸套的热负荷流体:进气、排气、缸内的流场、消声器(可用)车辆的振动、驾驶室噪声(可用),几个内燃机方面实例,潍柴6170连杆、曲轴,几个内燃机方面实例(1),潍

5、柴6170机体玉柴6108机体,几个内燃机方面实例(2),潍柴6160表面辐射噪声,几个内燃机方面实例,上汽通用五菱B15D,其他方面实例,广州五羊125摩托车架汽研中心汽车碰撞试验,其他方面实例(2),天津扎努西冰箱压缩机桥梁安全,其他方面实例(3),2002年3月6日伯克利大学环境工程A.Astanehasl 教授在美国众议院科学委员会作证时做了美国911世贸大厦破坏仿真(包括碰撞分析、热分析)认为:大楼倒塌的原因是因为内部钢梁上的防火物质在高温的情况下被破坏,导致失去支撑力量,其他方面实例(3),其他方面分析实例(4),1,CAD-CAE-CAM,FEA 商业软件,SAP(79年)-Al

6、gor:structural analysis program 最早NASTRAN:NASAstructure analysis(动力分析、航空航天NSYS:国内高校学生其它的软件:MARC、ABAQUS、ADINA详细请见网络论坛。,内燃机商业软件,工作条件:气门间歇性启闭冲击载荷(惯性力+气门弹簧力)凸轮面的接触力磨损,24,LMS.Sysnoise(LMS.vituallab):噪声分析有限元(FEA)、无限元(IFEA)、边界元(BEA)软件,第一章:弹性力学基础第二章:一个有限元法的简单例子第三章:平面问题的有限元法(最简单的三角形单元)第四章:平面问题的高次单元第五章:三维实体单元

7、的简介(选讲)第六章:平面温度应力有限元法简介(选讲),课程内容,其他,课时安排:总共32学时,讲课12学时,上机练习20学时 考试方式:1、上机2、期末随堂堂考试学习本课程后达到的要求 1、掌握有限元法的入门阶段的理论基础知识 2、了解有限元分析的基本步骤,学会用Patran软件中的一些基本建模功,能做一些简单零件的静强度分析和有限元模态分析,参考书,有限元法基础 蒋友凉 国防工业出版社有限元法在ICE上的应用 讲义 天津大学教材科有限元法在动力机械上的应用 机械工业出版社有限元及边界元法基础 北航出版社有限元法基础 蒋孝煜 清华大学出版社工程有限单元法基础 石油大学出版社车用发动机设计 吴

8、兆汉 国防工业出版社,第一章 弹性力学基础,本章目的,说明:推导出弹性力学问题的微分方程可用FEA解微分方程可用FEA求解弹性力学问题,非变形体-刚体,变形体-材料力学和弹性力学,材料力学和弹性力学的区别,研究对象:材料力学:研究对象杆、柱、梁、板,长度远远大于厚度、宽度。弹性力学:任意形状弹性变形体(5个基本假设),可喜的是内燃机大多零件是弹性体。研究方法:材料力学:对整个截面建立三组方程。梁弯曲截面变形假设 弹性力学:对任意点建立三组平衡方程。,弹性体的基本假设和变量,5个基本假设3组变量:位移、应力、应变,各种应力名称,任意截面的正应力、剪应力、主应力(第1、2、3主应力)最大拉应力、最

9、大剪应力、米塞斯应力、最大主应力、弯曲应力、扭转应力、拉压应力x、y、z方向应力切向应力和法向应力平面应力和平面应变,弹性体内一点应力状态,3个正应力+6个剪应力=9个应力考虑到剪应力互等,共有6个应力,全部确定了该点的应力状态,弹性体平衡方程,3个方向合力=0,3个方向合力矩=0,共可列6个方程,3个方向合力矩=0,实际上是剪应力互等,只有3个平衡方程3个方向合力=0的推导:,弹性体平衡方程,轮换xyz可得到其他2个方向方程,弹性体平衡方程,共3个平衡方程,弹性体几何方程,几何方程:弹性体任一点位移和应变的关系以x向应变为例推导,y、z及剪应变参照,弹性体几何方程,X方向的应变共有3个正应变

10、+3个剪应变=6个几何方程,用矩阵和向量表示物理方程,物理方程,广义虎克定律应力和应变之间的关系,只和材料有关(弹性模量和泊松比),弹力三组方程:方程数=变量数,变量3个正应力+3个剪应力+3个位移+6个应变=15个变量目前有3个几何方程+6个几何方程+6个物理方程=15个方程但只能得到一组通解,必须有边界条件才能有唯一特解,两类边界条件,力边界条件位移边界条件 用圣维南原理简化内燃机边界条件:曲轴、连杆、活塞的接触边界,弹性力学方程的求解,力法:位移法:(多用,本文用)教材:1-22到1-25页理论上弹性力学方程可求解实际上简单区域的可得理论解,复杂区域只能得到近似数值解FEA,弹性力学虚功

11、原理,弹性体在外力作用下变形后处于平衡状态,假设物体在此基础上发生一个任意小位移,外力在虚位移上做的功等于变形体接受的虚变形功。平面应力问题虚功方程公式,第二章 一个有限元法的简单例子,有限元解变形、应力、应变问题,弹性力学:位移解法。通过弹力三组方程,把应力、应变转换为位移,求解得到位移,再得到应变,再得到应力。求解步骤:1、将弹性体离散化为标准件(单元),得到单元和节点。单元间只能通过节点传递力,(只考虑节点力和节点位移)2、单元内假设一简单的位移函数,利用3组基本方程建立节点力和节点位移关系,结果是线性代数方程组 3、将边界条件(力和位移约束)转化到节点力和节点位移 4、线代方程组,得节

12、点位移-插值得任一点位移-应变应力,商用有限元软件求解步骤,1、建立几何模型2、离散化(网格划分)3、单元属性的指定(包含材料)3.1单元刚阵的生成3.2总体刚阵的组装4、边界条件的施加(4.1实用总刚阵的获得)5、分析求解6、读入分析结果7、查看分析结果,商用有限元软件单元类型,0维单元:集中质量单元、点弹簧、阻尼单元(接地)1维单元:拉压(扭转)杆、梁单元2维单元:板、壳单元(三角形、四边形)平面问题单元3维单元(4面体、5面体、6面体)用实例说明上述应用1阶单元和2阶单元(略),FEA中考虑的几个问题,1、离散化 单元维数单元形状单元大小、疏密(静态计算、瞬态计算)节点位置选取单元和节点

13、的编号单元质量的检查(长度边比、最小角度、翘曲),简单例子,刚度的概念 弹簧振子:刚度系数 2个弹簧振子刚度矩阵k:矩阵中的元素:刚度系数刚度系数-柔度系数,柔度矩阵,如图弹性体各点作用广义力F1、F2 Fn,各点位移为定义柔度系数 为在 j点作用的单位力时在 i产生的位移,则力 在i点产生的位移为,所有力在i产生的总位移则为则所有力在各点产生的位移为:,柔度矩阵,柔度矩阵C,刚度矩阵,如图弹性体各点作用广义力F1、F2 Fn,各点位移为定义柔度系数 为在 j点产生单位位移,需要在 i施加的力,则力在 j点产生 位移需要在 i施加的力为,所有点都产生位移,需要在 i加的力为则要在各点产生的位移

14、,需要在各点加的力为:,刚度矩阵,刚度矩阵k结论:如果知道刚度矩阵k,在力F作用下各点的位移就可以求得,其他各点位移可以插值得到,这就是有限元法的思想。关键是求得k,也就是,单元刚度矩阵,对单元建立此关系,单元刚度矩阵有限元种单元(标准件),如果能求得此关系。任意复杂件离散为这种标准件(单元),则大大简化运算FEA法的核心在于求得单元刚阵k,标准件之一:弹簧单元简例(1),用有限元法求拉伸杆的应力、应变和位移将杆离散化为2个单元,单元1长l1和单元2长l2,3个节点,1、2、3。假定节点1、2、3处的节点力为F1、F2和F3,位移为u1、u2和u3,单元1刚度系数k1,单元2为k2。取单元1,

15、有,标准件之一:弹簧单元简例(2),例子:建立弹簧单元的刚度矩阵用有限元法求拉伸杆的应力、应变和位移将杆离散化为2个单元,单元1和单元2,3个节点,1、2、3。假定节点1、2、3处的力为F1、F2和F3,位移为u1、u2和u3,单元1刚度系数k1,单元2为k2,取单元1,有至此已经得到单元1的单元刚阵,标准件之一:弹簧单元简例(3),单元1刚阵可写为(注意上标和下标)单元2刚阵可写为(注意上标和下标),标准件之一:弹簧单元简例(4),(1)+(2)为可写为,注意F1 F2 F3是节点力,是单元对节点的力这是整个结构的总体刚度矩阵,标准件之一:弹簧单元简例(5),施加边界条件 1、力:节点1 2

16、 3处分别受节点力合力,F1 F2 F3,还有节点3受外力P,考虑到各节点受力平衡,得出 F1=0 F2=0 F3=-P 2、位移:u1=0方程可写为,k1 k2 P已知,2个方程,2个未知量u1 u2可求位移已知,由几何方程可求应变,由物理方程可求应力可见弹力问题是微分方程问题,但最后解代数方程组,由单元刚阵组装总体刚阵,弹簧单元刚阵可以直接写出其他单元尚未求得,第3章和第4章将 由单元刚阵获得总体刚阵 在每个单元刚阵 基础上,把属于总体结构但不在本单元上的所有节点(编号)对应的行(第n行)、列(第n列)矩阵元素补0后,把所有单元刚阵相加即可得到总体刚阵。边界条件施加 1、力P:节点号对应的

17、力向量中元素为-P,其他为0 2、位移:去掉节点号对应的行、列,总体刚阵的特点,(1)对称性 总体刚阵元素关于对角线对称,可以从力的互等作用定理解释(2)任一行(列)元素代数和=0。从结构发生x向和y向刚体运动(此时节点力=0)解释(3)对角线元素为正(4)总体刚阵奇异,必须施加边界条件后才有唯一解,,第三章 平面问题有限元法,平面问题有限元法,有限元法的概念 具有初值、边值调节的微分方程的近似数值解法。求解步骤:1、将弹性体离散化单元和节点。2、单元内假设一简单的位移函数,利用3组基本方程建立节点力和节点位移关系,结果是线性代数方程组 3、将所有单元的线代方程装成整个结构的线性代数方程组 3

18、、将边界条件(力和位移约束)转化到节点力和节点位移 4、线代方程组,得节点位移-插值得任一点位移-应变应力,单元类型、形状和疏密,0维单元:集中质量单元、点弹簧、阻尼单元(接地)1维单元:拉压(扭转)杆、梁单元2维单元:板、壳单元(三角形、四边形)、平面问题单元3维单元(4面体、5面体、6面体)网格大小、疏密 体单元每个节点3个自由度,3个方程,整个结构计算时间大约和单元尺寸的3次方成反比单元节点数(1阶单元和高阶单元)3节点三角形单元和6节点三角形单元、4节点四面体单元和10节点四面体单元,8节点六面体单元(最常用)和20节点六面体单元,单元位移函数,单元位移函数:概念,一般为多项式,积分和

19、微分方便。单元有几个节点,多项式位移函数则有几个未知系数单元位移函数的确定(试着凑出单元节点力和节点力关系)1、单元的节点坐标是已知的 2、单元内任一点都满足位移函数 3、假定节点位移是已知量,目的建立节点力和节点位移关系 4、用待定系数法求解单元位移函数,3节点平面三角形单元(1),每个节点有x,y2个方向,单元有几个节点,位移函数则有几个未知系数。1、右图节点坐标已知i(xi,yi),j(xj,yj),k(xk,yk)2、假定节点位移已知:i:ui,vi,j:uj,vj,k:uk,vk3、单元节点也满足位移函数,代入单元位移函数4、未知量为a1-a6,6个方程 6个位置量,用待定系数法 可

20、求a1-a6,3节点平面三角形单元(2),求得a1-a6 把求得的a1-a6代入,可得三角形单元的位移函数牢记求解的目的在于获得公式:对本三角形单元来说为:,3节点平面三角形单元(3),把位移函数朝上式形式凑,整理成如下形式:可证明整理结果为:位移函数u,v都是1次多项式(1次函数)所以Ni,Nj,Nk一定是一次函数 位移函数u,v只和节点坐标、节点位移、x、y有关,所以 Ni,Nj,Nk一定是只和节点坐标有关的一次函数,而节点坐标决定单元三角形形状,因此Ni,Nj,Nk称为形函数,3节点平面三角形单元(4),把位移函数进一步整理为:对比本章的求解目的,公式右边已经有所接近,平面问题应变和几何

21、方程(1),平面问题应变 代入3节点平面三角形单元位移函数,得到其应变,平面问题应变和几何矩阵,代入3节点平面三角形单元位移函数,得到其应变上公式简写为 由于形函数Ni,Nj,Nk是一次函数,求导后为常数,所以矩阵B中元素一定为常数,也就是3节点三角形单元应变为常量,因此精度较差。对应几何方程,矩阵B 称为几何矩阵。,平面问题应力和弹性矩阵(3),考虑平面问题物理方程(广义虎克定律)矩阵D称为弹性矩阵,只和弹性模量E和泊松比u有关,也就是取决于材料。矩阵D和B中元素都为常量,所以3节点三角形单元为常应力单元,精度较差。,平面问题虚功方程,平面问题虚功方程(作为定律直接应用),t为平面厚度代入3

22、节点三角形单元的应变和应力得到,3节点平面三角形单元刚阵(1),三角形单元刚阵考虑到虚位移可以是任意量,提到积分号外考虑到节点位移也是任意的,3节点平面三角形单元刚阵(2),考虑到几何矩阵元素为常量,应力也是常量代入物理方程,应力考虑到节点位移也是任意的,3节点平面三角形单元刚阵(3),得到了三角形单元刚阵即建立了单元节点力和节点位移见的关系,对比上下公式,kij已经可求,平面三角形单元有限元法(1),1、如右图平面问题,划分网格,得到节点1、2、3、4和单元1、22、求出单元1和单元2的单元刚阵(6x6矩阵)3、将单元1和2的单元刚阵扩充相加,得总体刚阵。(8X8矩阵)4、施加边界条件 A、

23、节点1、2固定,位移为0,不用再求解,在总体刚阵中消去对应的行列(考虑到节点1、2都有x、y2个方向,所以消去4行4列,得实用总刚阵(4X4矩阵)B、由节点平衡方程,得知节点4x方向的节点力为-P,其他节点力=0,可列出节点力和节点位移的线性代数方程组(4个方程),平面三角形单元有限元法(2),求解线性代数方程组,得到各节点位移(实际上只有节点3、4的x、y向位移,共4个,节点1、2位移为0)节点位移已知后,平板上任一点的位移有位移函数求得(近似)知道位移后,平板上任一点应变有几何方程求得,对位移求导,常应变。平板上任一点应力有物理方程求得所有量都已经求得,平面三角形单元有限元法(3),至此所

24、有关于变形、应变、应力的平面问题都已经可以求解。考虑到三角形单元的常应变、常应力,细分网格提高精度通过其他形状单元提高精度,四边形单元,第四章 平面问题高阶单元,3节点平面三角形单元的局限,常应力、常应变单元,精度较差可以细分网格提高精度,但求解量增加,每个节点2个方程。采用高阶单元 1、6节点三角形单元:精度高,求解量大 2、矩形单元:精度高,求解量不大 3、任意四边形单元,矩形单元位移函数,位移函数:4个节点,4个系数参照三角形单元,用待定系数法求a1-a8得到位移函数u,v,整理成形函数的形式因为位移函数为2次函数,所以Ni,Nj,Nm,Np也为2次函数,形函数。,矩形单元应变,位移函数

25、:4个节点,4个系数矩形单元应变。因为形函数为2次函数,求导后几何矩阵B元素为一次函数。所以矩形单元应变是一次函数,矩形单元应力,矩形单元应力(物理方程)矩形单元应力。因为形函数为2次函数,求导后几何矩阵B元素为一次函数。所以矩形单元应力是一次函数三角形单元应变、应力都是常量,矩形单元应变和应力都是一次函数,精度高于三角形单元,可以采用较粗的网格就可以得到较高的精度,而且求解量小,平面问题虚功方程,平面问题虚功方程(作为定律直接应用),t为平面厚度代入矩形单元的应变得到,矩形单元刚阵(1),矩形单元刚阵考虑到虚位移可以是任意量,提到积分号外考虑到节点位移也是任意的,矩形单元刚阵(2),考虑到几

26、何矩阵元素为常量,应力也是常量对比以下公式,矩形单元刚阵已经就得矩形单元刚阵需要积分,几何矩阵B中元素为一次函数,积分方便,这也是单元位移函数采用多项式的原因,矩形单元应力,矩形单元应力(物理方程)矩形单元应力。因为形函数为2次函数,求导后几何矩阵B元素为一次函数。所以矩形单元应力是一次函数三角形单元应变、应力都是常量,矩形单元应变和应力都是一次函数,精度高于三角形单元,可以采用较粗的网格就可以得到较高的精度,而且求解量小,矩形单元有限元法(1),1、如右图平面问题,划分网格,得到节点1、2、3、4、5、6和单元1、22、求出单元1和单元2的单元刚阵(8x8矩阵)3、将单元1和2的单元刚阵扩充

27、相加,得总体刚阵。(12X12矩阵)4、施加边界条件 A、节点1、2固定,位移为0,不用再求解,在总体刚阵中消去对应的行列(考虑到节点1、2都有x、y2个方向,所以消去4行4列,得实用总刚阵(8X8矩阵)B、由节点平衡方程,得知节点4x方向的节点力为-P,其他节点力=0,可列出节点力和节点位移的线性代数方程组(8个方程),矩形单元有限元法(2),求解线性代数方程组,得到各节点位移(实际上只有节点3、4、5、6的x、y向位移,共8个,节点1、2位移为0)节点位移已知后,平板上任一点的位移有位移函数求得(近似)知道位移后,平板上任一点应变有几何方程求得,对位移求导,常应变。平板上任一点应力有物理方

28、程求得所有量都已经求得,矩形单元有限元法(3),可以只用1个矩形单元,4个节点,和三角形单元节点数相同,但,精度高于三角形单元,矩形单元的局限性,矩形单元不适于划分曲边边界处,需要借用三角形单元任意四边形单元适于划分曲边边界,不需要借用三角形单元,单元的变形协调性要求(1),相邻的两个单元变形后在公共边处既不能开裂,也不能重叠,这称为单元变形的协调性,这样的单元称为协调单元如果单元的位移函数在公共边上退化为一次函数(或低于一次),则单元变形协调1、三角形单元位移函数本身就是一次函数,所以单元变形协调 2、矩形单元位移函数为2次函数,但是在如下坐标系下,最右边 方程为x=a,在该边上位移函数退化

29、为一次函数,单元的变形协调性要求(2),3、任意四边形单元,公共边的位移函数仍然为2次函数,单元变形后,在相邻的公共边处会开裂或重叠,任意四边形变形不协调4、任意四边形采用自然坐标系可以解决单元变形协调性问题,内容略,第四章 体单元简介,三种形状体单元,四面体单元(棱锥)五面体单元(楔形单元),只用于四面体和六面体的过度六面体单元(砖形单元)内容略,4节点四面体单元位移函数,4个节点,4个未知系数位移函数为一次函数,类似3节点三角形单元,可以推知其单元的应变和单元应力都是常量,即4节点四面体单元是常应变、常应力单元,精度较差类似可得到4面体单元刚阵。(略)四面体单元法可以求解所有关于体的应力、

30、应变和位移问题(略),高阶体单元,10节点四面体单元8节点六面体单元20节点六面体单元,四面体与六面体单元的选择,4节点四面体单元常应变、常应力,精度差。细化网格才能提高精度,但会使得求解量增加10节点四面体单元应变、应力不是常量,精度高,但是求解量大。每个节点有3个方程8节点六面体单元不需要细化网格即可获得较高精度,求解量不大。四面体单元的最大优点是FEA软件自动划分网格容易,而六面体单元较难自动划分网格,投入的人力大大增加建议:静态应力计算可用10节点的4面体网格。动态应力计算和迭代计算量大的采用六面体网格,第五章 平面温度应力问题简介,平面温度场的控制方程,控制方程平面温度场温度边界条件

31、 1、第一类边界条件 2、第二类边界条件 3、第三类边界条件平面温度场问题涉及到泛函的极值函数,研究生数学内容(略),补充 弹簧单元,弹簧单元的位移函数和形函数,弹簧单元刚阵已知,k为弹性系数,可直接使用弹簧单元位移函数:2个节点,可确定2个系数,为一次函数假定节点位移u1,u2已知,待定系数法求a1和a2,得位移函数将位移函数整理成如下形式,矩阵中元素则为形函数,是一次函数,弹簧单元的应变和应力,知道位移后,应用几何方程,可得应变和几何矩阵可得弹簧单元应力利用虚功方程可得弹簧单元刚阵,弹簧单元有限元法,单元刚阵扩充相加,可得总体刚阵施加位移边界条件,可得到使用总刚阵施加力边界条件,利用平衡方程,可得节点力向量求解可得节点位移由位移函数可得弹簧单元内任一点位移由几何方程可以单元任一点应变由物理方程可得单元内任一点的应力至此应力、应变、位移都可求,上机练习 有限元法在内燃机上的应用,

展开阅读全文
相关资源
猜你喜欢
相关搜索
资源标签

当前位置:首页 > 生活休闲 > 在线阅读


备案号:宁ICP备20000045号-2

经营许可证:宁B2-20210002

宁公网安备 64010402000987号