欢迎来到淘文阁 - 分享文档赚钱的网站! | 帮助中心 好文档才是您的得力助手!
淘文阁 - 分享文档赚钱的网站
全部分类
  • 研究报告>
  • 管理文献>
  • 标准材料>
  • 技术资料>
  • 教育专区>
  • 应用文书>
  • 生活休闲>
  • 考试试题>
  • pptx模板>
  • 工商注册>
  • 期刊短文>
  • 图片设计>
  • ImageVerifierCode 换一换

    有限元分析理论基础.doc

    • 资源ID:13008645       资源大小:3.21MB        全文页数:25页
    • 资源格式: DOC        下载积分:9金币
    快捷下载 游客一键下载
    会员登录下载
    微信登录下载
    三方登录下载: 微信开放平台登录   QQ登录  
    二维码
    微信扫一扫登录
    下载资源需要9金币
    邮箱/手机:
    温馨提示:
    快捷下载时,用户名和密码都是您填写的邮箱或者手机号,方便查询和重复下载(系统自动生成)。
    如填写123,账号就是123,密码也是123。
    支付方式: 支付宝    微信支付   
    验证码:   换一换

     
    账号:
    密码:
    验证码:   换一换
      忘记密码?
        
    友情提示
    2、PDF文件下载后,可能会被浏览器默认打开,此种情况可以点击浏览器菜单,保存网页到桌面,就可以正常下载了。
    3、本站不支持迅雷下载,请使用电脑自带的IE浏览器,或者360浏览器、谷歌浏览器下载即可。
    4、本站资源下载后的文档和图纸-无水印,预览文档经过压缩,下载后原文更清晰。
    5、试题试卷类文档,如果标题没有明确说明有答案则都视为没有答案,请知晓。

    有限元分析理论基础.doc

    . .有限元分析概念有限元法:把求解区域看作由许多小的在节点处相互连接的单元(子域)所构成,其模型给出基本方程的分片(子域)近似解,由于单元(子域)可以被分割成各种形状和大小不同的尺寸,所以它能很好地适应复杂的几何形状、复杂的材料特性和复杂的边界条件有限元模型:它是真实系统理想化的数学抽象。由一些简单形状的单元组成,单元之间通过节点连接,并承受一定载荷。有限元分析:是利用数学近似的方法对真实物理系统(几何和载荷工况)进行模拟。并利用简单而又相互作用的元素,即单元,就可以用有限数量的未知量去逼近无限未知量的真实系统。线弹性有限元是以理想弹性体为研究对象的,所考虑的变形建立在小变形假设的基础上。在这类问题中,材料的应力与应变呈线性关系,满足广义胡克定律;应力与应变也是线性关系,线弹性问题可归结为求解线性方程问题,所以只需要较少的计算时间。如果采用高效的代数方程组求解方法,也有助于降低有限元分析的时间。线弹性有限元一般包括线弹性静力学分析与线弹性动力学分析两方面。非线性问题与线弹性问题的区别:1)非线性问题的方程是非线性的,一般需要迭代求解;2)非线性问题不能采用叠加原理;3)非线性问题不总有一致解,有时甚至没有解。有限元求解非线性问题可分为以下三类:1)材料非线性问题材料的应力和应变是非线性的,但应力与应变却很微小,此时应变与位移呈线性关系,这类问题属于材料的非线性问题。由于从理论上还不能提供能普遍接受的本构关系,所以,一般材料的应力与应变之间的非线性关系要基于试验数据,有时非线性材料特性可用数学模型进行模拟,尽管这些模型总有他们的局限性。在工程实际中较为重要的材料非线性问题有:非线性弹性(包括分段线弹性)、弹塑性、粘塑性及蠕变等。2)几何非线性问题几何非线性问题是由于位移之间存在非线性关系引起的。当物体的位移较大时,应变与位移的关系是非线性关系。研究这类问题一般都是假定材料的应力和应变呈线性关系。它包括大位移大应变及大位移小应变问题。如结构的弹性屈曲问题属于大位移小应变问题,橡胶部件形成过程为大应变问题。3)非线性边界问题在加工、密封、撞击等问题中,接触和摩擦的作用不可忽视,接触边界属于高度非线性边界。平时遇到的一些接触问题,如齿轮传动、冲压成型、轧制成型、橡胶减振器、紧配合装配等,当一个结构与另一个结构或外部边界相接触时通常要考虑非线性边界条件。实际的非线性可能同时出现上述两种或三种非线性问题。有限元理论基础有限元方法的基础是变分原理和加权余量法,其基本求解思想是把计算域划分为有限个互不重叠的单元,在每个单元内,选择一些合适的节点作为求解函数的插值点,将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式,借助于变分原理或加权余量法,将微分方程离散求解。采用不同的权函数和插值函数形式,便构成不同的有限元方法。1.加权余量法:是指采用使余量的加权函数为零求得微分方程近似解的方法称为加权余量法。(Weighted residual method WRM)是一种直接从所需求解的微分方程及边界条件出发,寻求边值问题近似解的数学方法。加权余量法是求解微分方程近似解的一种有效的方法。设问题的控制微分方程为: 在V域内 在S边界上式中 : L、B分别为微分方程和边界条件中的微分算子; f、g 为与未知函数u无关的已知函数域值;u为问题待求的未知函数混合法对于试函数的选取最方便,但在相同精度条件下,工作量最大。对内部法和边界法必须使基函数事先满足一定条件,这对复杂结构分析往往有一定困难,但试函数一经建立,其工作量较小。无论采用何种方法,在建立试函数时均应注意以下几点:(1)试函数应由完备函数集的子集构成。已被采用过的试函数有幂级数、三角级数、样条函数、贝赛尔函数、切比雪夫和勒让德多项式等等。 (2)试函数应具有直到比消除余量的加权积分表达式中最高阶导数低一阶的导数连续性。 (3)试函数应与问题的解析解或问题的特解相关联。若计算问题具有对称性,应充分利用它。显然,任何独立的完全函数集都可以作为权函数。按照对权函数的不同选择得到不同的加权余量计算方法,主要有:配点法、子域法、最小二乘法、力矩法和伽辽金法。其中伽辽金法的精度最高。2、虚功原理 平衡方程和几何方程的等效积分“弱”形式虚功原理包含虚位移原理和虚应力原理,是虚位移原理和虚应力原理的总称。他们都可以认为是与某些控制方程相等效的积分“弱”形式。虚功原理:变形体中任意满足平衡的力系在任意满足协调条件的变形状态上作的虚功等于零,即体系外力的虚功与内力的虚功之和等于零。虚位移原理是平衡方程和力的边界条件的等效积分的“弱”形式;虚应力原理是几何方程和位移边界条件的等效积分“弱”形式。虚位移原理的力学意义:如果力系是平衡的,则它们在虚位移和虚应变上所作的功的总和为零。反之,如果力系在虚位移(及虚应变)上所作的功的和等于零,则它们一定满足平衡方程。所以,虚位移原理表述了力系平衡的必要而充分条件。一般而言,虚位移原理不仅可以适用于线弹性问题,而且可以用于非线性弹性及弹塑性等非线性问题。虚应力原理的力学意义:如果位移是协调的,则虚应力和虚边界约束反力在他们上面所作的功的总和为零。反之,如果上述虚力系在他们上面所作的功的和为零,则它们一定是满足协调的。所以,虚应力原理表述了位移协调的必要而充分条件。虚应力原理可以应用于线弹性以及非线性弹性等不同的力学问题。但是必须指出,无论是虚位移原理还是虚应力原理,他们所依赖的几何方程和平衡方程都是基于小变形理论的,他们不能直接应用于基于大变形理论的力学问题。3、最小总势能法应变能:作用在物体上的外载荷会引起物体变形,变形期间外力所做的功以弹性能的形式储存在物体中,即为应变能。由n个单元和m个节点组成的物体的总势能为总应变能和外力所做功的差:最小势能原理:对于一个稳定的系统,相对于平衡位置发生的位移总会使系统的总势能最小,即:,i=1,2,3,n有限元法的收敛性有限元法是一种数值分析方法,因此应考虑收敛性问题。有限元法的收敛性是指:当网格逐渐加密时,有限元解答的序列收敛到精确解;或者当单元尺寸固定时,每个单元的自由度数越多,有限元的解答就越趋近于精确解。有限元的收敛条件包括如下四个方面:1)单元内,位移函数必须连续。多项式是单值连续函数,因此选择多项式作为位移函数,在单元内的连续性能够保证。2)在单元内,位移函数必须包括常应变项。每个单元的应变状态总可以分解为不依赖于单元内各点位置的常应变和由各点位置决定的变量应变。当单元的尺寸足够小时,单元中各点的应变趋于相等,单元的变形比较均匀,因而常应变就成为应变的主要部分。为反映单元的应变状态,单元位移函数必须包括常应变项。3)在单元内,位移函数必须包括刚体位移项。一般情况下,单元内任一点的位移包括形变位移和刚体位移两部分。形变位移与物体形状及体积的改变相联系,因而产生应变;刚体位移只改变物体位置,不改变物体的形状和体积,即刚体位移是不产生变形的位移。空间一个物体包括三个平动位移和三个转动位移,共有六个刚体位移分量。由于一个单元牵连在另一些单元上,其他单元发生变形时必将带动单元做刚体位移,由此可见,为模拟一个单元的真实位移,假定的单元位移函数必须包括刚体位移项。4)位移函数在相邻单元的公共边界上必须协调。对一般单元而言,协调性是指相邻单元在公共节点处有相同的位移,而且沿单元边界也有相同的位移,也就是说,要保证不发生单元的相互脱离开裂和相互侵入重叠。要做到这一点,就要求函数在公共边界上能由公共节点的函数值唯一确定。对一般单元,协调性保证了相邻单元边界位移的连续性。但是,在板壳的相邻单元之间,还要求位移的一阶导数连续,只有这样,才能保证结构的应变能是有界量。总的说来,协调性是指在相邻单元的公共边界上满足连续性条件。前三条又叫完备性条件,满足完备条件的单元叫完备单元;第四条是协调性要求,满足协调性的单元叫协调单元;否则称为非协调单元。完备性要求是收敛的必要条件,四条全部满足,构成收敛的充分必要条件。在实际应用中,要使选择的位移函数全部满足完备性和协调性要求是比较困难的,在某些情况下可以放松对协调性的要求。需要指出的是,有时非协调单元比与它对应的协调单元还要好,其原因在于近似解的性质。假定位移函数就相当于给单元施加了约束条件,使单元变形服从所加约束,这样的替代结构比真实结构更刚一些。但是,这种近似结构由于允许单元分离、重叠,使单元的刚度变软了,或者形成了(例如板单元在单元之间的绕度连续,而转角不连续时,刚节点变为铰接点)对于非协调单元,上述两种影响有误差相消的可能,因此利用非协调单元有时也会得到很好的结果。在工程实践中,非协调元必须通过“小片试验后”才能使用。应力的单元平均或节点平均处理方法最简单的处理应力结果的方法是取相邻单元或围绕节点各单元应力的平均值。 1.取相邻单元应力的平均值这种方法最常用于3节点三角形单元中。这种最简单而又相当实用的单元得到的应力解在单元内是常数。可以将其看作是单元内应力的平均值,或是单元形心处的应力。由于应力近似解总是在精确解上下振荡,可以取相邻单元应力的平均值作为此两个单元合成的较大四边形单元形心处的应力。如2单元的情况下,取平均应力可以采用算术平均,即平均应力=(单元1的应力+单元2的应力)/2。也可以采用精确一些的面积加权平均,即平均应力=单元1应力×单元1的面积+单元2应力×单元2面积/(单元1面积+单元2面积)当相邻两单元面积相差不大时,两者的结果基本相同。在单元划分时应避免相邻两单元的面积相差太多,从而使求解的误差相近。一般而言,3节点三角形单元的最佳应力点是单元的中心点,此点的应力具有1阶的精度。 2.取围绕节点各单元应力的平均值首先计算围绕该节点(i)周围的相关单元在该节点出的应力值,然后以他们的平均值作为该节点的最后应力值,即其中,1m是围绕在i节点周围的全部单元。取平均值时也可进行面积加权。有限元法求解问题的基本步骤1.结构离散化对整个结构进行离散化,将其分割成若干个单元,单元间彼此通过节点相连;2.求出各单元的刚度矩阵K(e)K(e)是由单元节点位移量(e)求单元节点力向量F(e)的转移矩阵,其关系式为:F(e)= K(e) (e)3.集成总体刚度矩阵K并写出总体平衡方程:总体刚度矩阵K是由整体节点位移向量求整体节点力向量 的转移矩阵,其关系式为F= K ,此即为总体平衡方程。4.引入支撑条件,求出各节点的位移节点的支撑条件有两种:一种是节点n沿某个方向的位移为零,另一种是节点n沿某个方向的位移为一给定值。5.求出各单元内的应力和应变。 对于有限元方法,其基本思路和解题步骤可归纳为:(1)建立积分方程,根据变分原理或方程余量与权函数正交化原理,建立与微分方程初边值问题等价的积分表达式,这是有限元法的出发点。(2)区域单元剖分,根据求解区域的形状及实际问题的物理特点,将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作,这部分工作量比较大,除了给计算单元和节点进行编号和确定相互之间的关系之外,还要表示节点的位置坐标,同时还需要列出自然边界和本质边界的节点序号和相应的边界值。 (3)确定单元基函数,根据单元中节点数目及对近似解精度的要求,选择满足一定插值条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的,由于各单元 具有规则的几何形状,在选取基函数时可遵循一定的法则。 (4)单元分析:将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近;再将 近似函数代入积分方程,并对单元区域进行积分,可获得含有待定系数(即单元中各节点 的参数值)的代数方程组,称为单元有限元方程。 (5)总体合成:在得出单元有限元方程之后,将区域中所有单元有限元方程按一定法则进行累加,形成总体有限元方程。 (6)边界条件的处理:一般边界条件有三种形式,分为本质边界条件(狄里克雷边界条件 )、自然边界条件(黎曼边界条件)、混合边界条件(柯西边界条件)。对于自然边界条件, 一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件,需按一定法 则对总体有限元方程进行修正满足。 (7)解有限元方程:根据边界条件修正的总体有限元方程组,是含所有待定未知量的封闭 方程组,采用适当的数值计算方法求解,可求得各节点的函数值。 单元刚度矩阵的特性单元刚度矩阵无论在局部坐标系中还是在整体坐标系中都具有相同的三个特性:1)对称性由材料力学中的位移互等定理可知,对一个构件,作用在点j的力引起点i的绕度等于有同样大小而作用于点i的力引起的点j的绕度,即kij(e)= kji(e),表明单元刚度矩阵是一个对称矩阵。2)  奇异性无逆阵的矩阵就叫做奇异矩阵,其行列式的值为0,即|k(e)|=0,这一点可以从例题直接得到验证。其物理意义是引入支撑条件之前,单元可平移。3)  分块性有前面所讲的内容可以看出,矩阵k(e)可以用虚线分成四块,因此可写成如下的分块形式,式中kmn(e)局部坐标系中单元(e)按局部码标记的节点m、n之间的刚度子矩阵刚架结构中非节点载荷的处理的方法在刚架结构以及其他较复杂的结构上,他们所受的载荷可以直接作用在节点上,又可以不直接作用在节点上而作用于单元节点间的其他位置上。后一种情况下的载荷称为非节点载荷。有限元分析时,总体刚度方程中所用到的力向量 是节点力向量。因此在进行整体分析前应当进行载荷的移植,将作用于单元上的力移植到节点上。移植时按静力等效的原则进行。处理非节点载荷一般可直接在整体坐标系内进行,其过程为:1)将各杆单元看成一根两端固定的梁,分别求出两个固定端的约束反力。其结果可直接利用材料力学的公式求得;2)将各固定端的约束反力变号,按节点进行集成,获得各节点的等效载荷总体刚度矩阵的集成法使用刚度矩阵获得的方法获得总体刚度矩阵。在此将其扩展到由整体坐标系中的单元刚度矩阵的子矩阵集成总体刚度矩阵。步骤如下:1)对一个有n个节点的结构,将总体刚度矩阵K划分为n×n各子区间,然后按节点总码的顺序进行编号;2)将整体坐标系中单元刚度矩阵的各子矩阵根据其下标的两个总码对号入座,写在总体刚度矩阵相应的子区间;3)同一子区间内的子矩阵相加,成为总体刚度矩阵中的相应的子矩阵。总体刚度矩阵的特性1)对称性:因为由此特性,在计算机中只需存储其上三角部分;2)奇异性:物理意义仍为在无约束的情况下,整个结构可做刚体运动;3)稀疏性:K中有许多零子矩阵,而且在非零子矩阵中还有大量的零元素,这种矩阵称为稀疏矩阵。大型结构的总体刚度矩阵一般都是稀疏矩阵;4)分块性:平面问题离散化时的规定1)单元之间只在节点处相连;2)所有的节点都为铰接点;3)单元之间的力通过节点传递;4)外载荷都要移植到节点上;5)在节点位移或某一分量可以不计之处,就必须在该节点安置一个铰支座或相应的连杆支座。通过以上的规定来建立平面有限元分析模型。结构对称性的利用规律一般来说,作用在对称结构上的载荷系统分为对称的、反对称的和一般的三种情况。1.结构对称,载荷对称或反对称这种情况下,对称面上的边界条件可按以下规则确定:A.在不同的对称面上,将位移分量区分为对称分量和反对称分量;B.将载荷也按不同的对称面分别区分为对称分量和反对称分量;C.对于同一个对称面,如载荷是对称的,则对称面上位移的反对称分量为零,如载荷是反对称的,则对称面上位移的对称分量为零。如果所分析的结构对称,但载荷是不对称的,也不是反对称的,这时可以将这种结构系统简化成载荷为对称和/或反对称情况的组合,仍可以简化分析过程,提高分析的综合效率。如图a所示,结构对称,载荷一般,可将其载荷分解为图b和图c的组合。图b为对称结构,载荷对x、y轴均为对称,图c为结构对称,载荷对x轴反对称、对y轴对称,此时可取相同的四分之一进行研究,分别施加对称面上节点的边界条件,进行两次分析计算,并将计算结果迭加起来,即可得到原结构四分之一的解答,进而得出整个结构的解答。利用结构的对称性取某一部分建立有限元模型时,往往会产生约束不足现象。例如,若取上例中图c的四分之一建立有限元时,根据上述分析,在两对称面上应加水平放置的滚动铰支座,因此模型在垂直方向存在刚体位移。对这种约束不足问题,利用有限元分析时,必须增加附加约束,以消除模型的刚体位移。在本例中,垂直方向可以用刚度很小的杆单元或边界弹簧单元连接到模型某节点上,使得既消除了模型的刚体位移,又不致于因附加的杆单元或边界弹簧单元刚度太大而影响结构原有的变形状态。单元形态的选择原则单元形态包括单元形状、边中节点的位置、细长比等,在结构离散化过程中必须合理选择。一般来说,为了保证有限元分析的精度,必须是单元的形态尽可能的规则。对于三角形单元,三条边长尽量接近,不应出现大的钝角、大的边长。这是因为根据误差分析,应力和位移的误差都和单元的最小内角的正弦成反比。因而,等边三角形单元的形态最好,它与等腰直角三角形单元的误差之比为sin45°:sin60°=1:1.23。但是为了适应弹性体边界,以及单元由小到大逐渐过渡,不可能是所有的三角形单元都接近等边三角形。实际上,常常使用等腰直角三角形。对于矩形单元来说,细长比不宜过大。细长比是指单元最大尺寸和最小尺寸之比。最优细长比在很大程度上取决于不同方向上位移梯度的差别。梯度较大的方向,单元尺寸要小些,梯度小的方向,单元尺寸可以大一些;如果各方向上位移梯度大致相同,则细长比越接近1,精度越高。有文献推荐,一般情况下,为了得到较好的位移结果,细长比不应超过7;为了获得较好的应力结果,细长比不应超过3。一般情况下,正方形单元的形态最好。对于一般的四边形单元应避免过大的边长比,过大的边长比会导致病态的方程组。边界条件的确定确定边界条件是建立有限元模型的重要一环,合理确定有限元模型的边界条件是成功地进行结构有限元分析的基本要求。一般情况下,建模对象的边界条件是明确的。根据力学模型的边界条件可以很容易确定其有限元模型的边界条件。例如电线杆插入地基的一端为固定端,桥梁一端为固定铰支座,另一端为滚动较支座。但是,在机械工程中,建模对象往往是整个结构中的一部分,在建立有限元模型,确定其边界条件时,必须考虑其余部分的影响。这方面主要考虑如下两类问题。1.边界位置的确定在建立连续弹性体局部区域的有限元模型时,往往取该局部区域为隔离体,取其隔离边界条件为零位移约束,并通过试探校正确定零位移边界条件的位置。例如,进行齿轮齿有限元分析时,取一个轮齿的局部区域为隔离体,如图所示,设定PQRS的边界条件为零位移约束,通过改变边界深度PQ和边界宽度PS研究边界位置对齿根最大拉应力的影响,最后确定合理的边界条件。2.边界条件的确定有些分析对象的边界位置是零部件的连接部位。在建立有限元模型时,必须研究如何给定边界位置上的边界条件,以反映相连接结构的影响。确定这种问题的边界条件是用简单支撑连杆替代相连接结构的作用,使替代后结构的系统刚度等价于原结构的系统刚度。如分析机床主轴和传动轴时,可以利用等刚度的杆单元替代轴承和支座的作用,使轴的分析中包含有轴承和支座的影响。单元和节点编号规则当利用整体刚度矩阵的带状特征进行存贮和求解方程组时,单元节点编号直接影响系统整体刚度矩阵的半带宽,也就是影响在计算机中存贮信息的多少、计算时间和计算费用。因而,要求合理的节点编号使带宽极小化。半带宽的计算公式:半带宽d=(单元节点号的最大差值+1)×节点自由度由此,进行网格节点编号时应使网格中单元节点号的最大差值最小,这样才能保证半带宽最小。试比较下图。图所示网格的四种编号方案中,单元节点标号的最大差值分别为5,3,5,9。显然,图2方案要合理。由此得出结论:沿着短边方向按列-列-列-列地顺序编号比沿着长度方向按行-行-行-行地顺序要合理(半带宽小)平面问题中非节点载荷转换为等效节点载荷由于三角形单元复杂的力学性质,不能像分析刚架时那样简单地利用材料力学公式来求解,而要用虚功方程将加在结构上的非节点载荷转换为等效节点载荷。掌握以下两种常见的非节点载荷的移植结果。1)作用在单元一条侧边上的集中力设Q平行于x方向,如图4-14所示,则等效节点载荷为 若Q平行于y方向,结果与此相仿。2)作用在单元一条侧边上呈三角形 分布的载荷设载荷平行于x方向,如图4-15所示,则等效节点载荷为若分布载荷为集度是q的均布载荷,则其余分量为零。求解时模型是否准备就绪?在求解初始化前,应进行分析数据检查,包括下面内容:1.统一的单位;2.单元类型和选项;3.材料性质参数:考虑惯性时应输入材料密度;热应力分析时应输入材料的热膨胀系数;4.实常数 (单元特性);5.单元实常数和材料类型的设置;6.实体模型的质量特性 (Preprocessor > Operate > Calc Geom Items);7.模型中不应存在的缝隙;8.壳单元的法向;9.节点坐标系;10.集中、体积载荷面力方向;11.温度场的分布和X围;12.热膨胀分析的参考温度。. .word.

    注意事项

    本文(有限元分析理论基础.doc)为本站会员(知****量)主动上传,淘文阁 - 分享文档赚钱的网站仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对上载内容本身不做任何修改或编辑。 若此文所含内容侵犯了您的版权或隐私,请立即通知淘文阁 - 分享文档赚钱的网站(点击联系客服),我们立即给予删除!

    温馨提示:如果因为网速或其他原因下载失败请重新下载,重复下载不扣分。




    关于淘文阁 - 版权申诉 - 用户使用规则 - 积分规则 - 联系我们

    本站为文档C TO C交易模式,本站只提供存储空间、用户上传的文档直接被用户下载,本站只是中间服务平台,本站所有文档下载所得的收益归上传人(含作者)所有。本站仅对用户上传内容的表现方式做保护处理,对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私,请立即通知淘文阁网,我们立即给予删除!客服QQ:136780468 微信:18945177775 电话:18904686070

    工信部备案号:黑ICP备15003705号 © 2020-2023 www.taowenge.com 淘文阁 

    收起
    展开