浙大化工机械-复杂系统分析与建模-课程作业1

复杂系统分析及建模课程大作业余雏麟 110280021、计算条件1.1 结构示意图容器的结构示意图如图 1 所示：?65012 1220R10120?46R6R1020 120R1516?54 36120R1020R10 R325R6120 20R1516?54?46图 1 容器结构示意图其中球型封头与筒体的过渡采用堆焊球封头的方式。1.2 计算要求（1）运用有限元法预测爆破压力 bP（2）若容器壁温从 200降至室温 20，同时内压 满足如下条件，p求疲劳寿命，并说明热应力利弊。sin,0432,8bpth2、容器爆破过程分析塑性材料制造的压力容器的爆破过程如图 2 所示。在弹性变形阶段（OA 线段） ，器壁应力较小，产生弹性变形，内压与容积变化量成正比，到 A 点时容器内表面开始屈服，与 A 点对应的压力为初始屈服压Py PbAB C变 形 量pO图 2 厚壁圆筒中压力与变形关系力 ；在弹塑性变形阶段（AC 线段） ，随着内压的继续提高，材料从内壁向外壁sp屈服，此时，一方面因塑性变形而使材料强化导致承压能力提高，另一方面因壁厚不断减薄而使承压能力下降，但材料强化作用大于壁厚减薄作用，到 C 点时两种作用已接近，C 点对应的压力是容器所能承受的最大压力，称为塑性垮塌压力；在爆破阶段（CD 线段） ，容积突然急剧增大，使容器继续膨胀所需要的压力也相应减小，压力降落到 D 点，容器爆炸，D 点所对应的压力为爆破压力 。bp对于内压容器，爆破过程中，内压和容积变化量的关系与材料塑性、容器容积和壁厚有关。对于脆性材料，不会出现弹塑性变形阶段。虽然塑性垮塌压力大于爆破压力，但工程上往往把塑性垮塌压力视为爆破压力。本报告也近似认为塑性垮塌压力就是爆破压力。3、爆破压力的经验计算公式Faupel 厚壁圆筒爆破压力计算公式：，2()ln3sbsbpKoiR, ， ，50sMPa530a361o25i6.87bp4、爆破压力数值计算4.1 计算简化根据结构的特点，从偏于危险和节约计算成本的角度，建模时只考虑整体的1/4，几何模型如图 2 所示:图 2 几何模型4.2 材料属性、单元选择及边界条件材料为 Q345R，弹性模量 205000MPa，泊松比 0.3，热膨胀系数 ，实际1.26e应力应变关系曲线如图 3 所示：图 3 Q345 真实应力应变关系曲线采用有限元软件 Abaqus6.9 进行数值模拟计算，单元采用具有沙漏控制的减缩积分单元 C3D8R。有限元模型的边界条件如图 4 所示：图 4 边界条件4.3 网格收敛性验证：对结构采用两种不同大小的网格进行划分，以爆破压力作为评价标准，考察网格的收敛性。本分析采用的网格模型和较细的网格划分模型的网格图片分别如图 5所和图 6 所示。它们相应的爆破压力计算结果如表 1 所示。图 5 本报告采用网格图 6 较细的网格表 1网格模型 单元数 爆破压力本报告采用的网格 19345 49.40较细网格 24611 49.41从表 1 可知，本报告采用的网格划分符合收敛性条件。4.4 数值计算分别采用弧长法和非线性稳定性算法计算容器的爆破压力。4.4.1 弧长法采用弧长法得到的结构在爆破时的 Mises 应力云图如图 7 所示。容器的爆破压力为 49.40MPa。爆破压力与弧长关系曲线如图 8 所示：图 7 容器爆破时 Mises 应力云图图 8 爆破压力与弧长关系曲线4.4.2 非线性稳定性算法采用 ABAQUS 的 dynamic/implicit 模块，运用非线性稳定算法得到的结构的Mises 应力云图如图 9 所示，载荷与弧长关系曲线如图 10 所示：图 9 非线性稳定算法得到的容器 Mises 应力云图图 10 爆破压力曲线图得到容器的爆破压力为 49.33MPa。4.5 有无开孔及开孔大小对爆破压力的影响为研究有无开孔及开孔大小对爆破压力的影响，考虑如表 2 所示的几个模型：表 2模型名称 模型说明模型 1 原始模型模型 2 仅封头开孔模型 3 无开孔模型 4 开孔尺寸变大计算结果如表 3 所示：表 3模型名称 爆破压力（MPa） 爆破位置模型 1 49.40 筒体接管相贯处内壁模型 2 49.41 封头接管相贯处内壁模型 3 49.45 筒体模型 4 48.93 筒体接管相贯处内壁模型 2、模型 3 和模型 4 在爆破时相应的 Mises 应力云图分别如图 11、12 和 13所示。图 11 模型 2 爆破时 Mises 应力云图图 12 模型 3 爆破时 Mises 应力云图图 13 模型 4 爆破时 Mises 应力云图在前面的计算结果可以看出，在本文采用的计算模型的前提下，有无开孔，及开孔大小对爆破压力影响很小，但是对爆破位置却影响很大。4.6 倒角大小对爆破压力的影响在原始模型的基础上，研究筒体与接管相贯处内倒角大小对爆破压力的影响，考虑如下的几个模型：表 4模型名称 模型说明模型 1 原始模型倒角 R=10mm模型 5 倒角 R=15mm模型 6 倒角 R=5mm计算结果如表 5 所示：表 5爆破压力（MPa） 爆破位置49.40 筒体接管相贯处内壁49.399 筒体接管相贯处内壁49.41 筒体接管相贯处内壁模型 5 和模型 6 在爆破时相应的 Mises 应力云图分别如图 14 和 15 所示。图 14 模型 5 爆破时 Mises 应力云图图 15 模型 6 爆破时 Mises 应力云图对于该容器来说，倒角的大小变化对爆破压力影响很小，但对爆破位置影响很大。5、极限载荷数值计算5.1 极限分析方法极限分析的基本概念：在加载过程中，结构中的高应力区首先进入塑性，当载荷继续增加时塑性区便不断扩大，同时还出现应力重分布现象。当载荷增大到某一极限值时，由理想塑性材料制成的结构将变成不稳定的几何可变机构，从而丧失承载能力，出现不可限制的塑性流动，此时载荷不变但应变能无限增加，这种状态称为塑性极限状态，相应的载荷称为极限载荷（全屈服载荷） 。假如只计算极限载荷而不计及极限状态到达以前的变形过程，不考虑加载的历史，这种分析方法就称为极限分析法。极限分析基本假设：材料为理想弹塑性，变形足够小，满足比例加载条件。极限载荷确定的方法：采用清华大学陆明万教授等提出的零曲率法，该方法具有分散性小的特点。5.2 材料属性材料基本参数：材料为 Q345R，弹性模量 205000MPa，泊松比 0.3，初始屈服应力。30sMPa5.3 极限载荷的理论计算值筒体整体屈服时，按 Mises 屈服失效判据，极限载荷理论计算值为：2ln36.9spKa5.4 数值计算分别采用 Newton 法和弧长法计算容器的极限载荷。5.4.1Newton 法在容器内壁所有与介质接触的表面施加压力载荷 p=60MPa，采用 Newton 法计算得到的应力云图如图 16 所示：图 16 极限分析应力云图载荷-最大点等效塑性应变曲线，如图 17 所示：0. 0.10.20.30.40.50.60.0.10.20.30.40.50.6p(MPa) 图 17 LPF 曲线图采用零斜率准则确定极限载荷： 600.610404=36.62 MPa36.24.15MPa允 许 极 限 载 荷5.4.2 弧长法在容器内壁所有与介质接触的表面施加压力载荷 p=60MPa，采用 risk 方法计算得到的应力云图如图 18 所示：图 18 极限分析容器的 Mises 应力云图载荷-最大点等效塑性应变曲线，如图 19 所示：图 19 LPF- 曲线图采用零斜率准则确定极限载荷： 600.592061=35.52MPa35.2.681允 许 极 限 载 荷 MPa6、应力分类法6.1 应力分类的概念一次应力超过材料屈服极限时，将会引起过量的总体塑性变形而造成结构破坏。分析设计一次加二次应力强度的控制值采用了 ，即以结构是否安定为判断2s依据。6.2 材料属性材料基本参数：材料为 Q345R，弹性模量 205000MPa，泊松比 0.3。6.3 名义弹性应力计算方法材料为线弹性，取单位压力 1 ，计算得到的应力云图如图 20 所示：MPa图 20 名义弹性应力分析应力云图及分析路径选取关键路径并对其进行应力评定，此处仅列出最危险路径的应力线性化数据：表 6应力分类 计算值 MPa（ ） 许用值 Pa（ ） 许 用 值允 许 压 力 计 算 值 MPa（ ）lP11.92 238.5 20.01lbQ19.14 477 24.92采用应力分类法计算得到的允许压力为 20.01 。MPa7、疲劳寿命计算7.1 基本概念工作循环：由初始状态进入新状态，随后又回到初始状态开始点的过程。包括启动停止循环、正常工作循环和设计中必须考虑的任何紧急状态或异常情况由起始到恢复的循环。本次计算时仅仅考虑正常工作循环。应力循环：指应力由初始值开始，经过代数最大值和代数最小值，然后又返回初始值的循环。一个工作循环可以引起一个或多个应力循环。根据容器的工作条件，一个工作循环包括两个应力循环。疲劳评定是以结构应力循环引起的应力差波值为基础的。根据工作循环和应力循环的定义，结合本题给定的条件。此工作循环包括 2 个应力循环。疲劳评定方法：当所考虑点的主应力方向在循环中变化时，应按如下步骤确定交变应力强度幅：（1）确定所考虑点在整个应力循环中与时间相对应的包括总体和局部结构不连续以及热效应所引起的六个应力分量；（2）选取循环条件的极端点（代数最大值或最小）所对应的时刻，记下这一时刻的各应力分量；（3）循环中的每一时刻，从每个与该时刻对应的应力分量中减去 i 时刻的相应的应力分量记下波动应力分量；（4）在循环的每一时刻，计算由六个波动应力分量所导出的主应力波动范围。这些波动主应力的方向虽然在循环中变化，但其编号应保持不变；（5）计算在整个应力循环中，相对于时间的波动主应力差；（6）确定各波动主应力差的最大波动范围，计算交变应力强度幅。7.2 内加热对容器疲劳的影响分析对于厚壁圆筒在内压与内加热同时作用下的应力特点来说，如果由内压引起的应力与内加热温差所引起的热应力同时存在，在弹性变形前提下筒壁的总应力则为两种应力的叠加。内加热情况下内壁应力叠加后得到改善，而外壁应力有所恶化。外加热时则相反，内壁应力恶化，而外壁应力得到很大改善。通常内压作用下应力最大点都出现在接管和筒体相贯处的内壁，疲劳评定时也通常选取此点。因此从疲劳的角度分析，内加热时对冲内压的容器的疲劳是有益的。相对应的是对与容器受外压的情况，内加热对容器的疲劳是不利的。7.3 采用 ABAQUS 软件对容器进行的疲劳计算与评定分析方法：考虑热应力影响时，采用间接耦合方法，即先计算温度场再计算应力场的方法。温度场计算时仅仅考虑热传导。考虑几何非线性效应。采用固定时间步长以 0.2小时为单位增量单位。对容器内壁施加随时间变化的温度载荷，外壁保持温度恒定为 20，且仅考虑热传导效应。7.3.1 不考虑温度时的容器疲劳评定不考虑温度时容器在第 1 小时和第 3 小时的应力云图如图 21 所示：图 21 第 1 小时和第 3 小时对应的容器名义弹性应力云图以节点 1389 为疲劳评定点，相关应力数据如表 7 所示：表 7应力分量 任意时刻应力 最大或最小应力 应力差(MPa) (MPa) (MPa)S11 19.1852 -19.1852 38.3704S22 381.043 -381.043 762.086S33 24.1375 -24.1375 48.275S12 30.7763 -30.7763 61.5526S13 7.20228 -7.20228 14.40456S23 -5.99224 5.99224 -11.9845相应三个主应力值为 767.44，57.72，25.56(MPa)主应力波动差值为 767.44-25.56=741.88MPa 对应的交变应力幅为 370.1MPa采用 JB4732 疲劳载荷计算公式：对应的应力循环次数约为 3472 次。采用 ansys 的疲劳模块，得到：评定点在事件中各载荷点的应力值：LIST FATIGUE STRESSES FOR LOCATIONS 1 TO 5 IN STEPS OF 1EVENTS 1 TO 3 - LOADS 1 TO 10* POST1 FATIGUE STRESSES *LOCATION 1 NODE 1389 EVENT 1 (TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 1SX,SY,SZ 19.180 381.04 24.140 SXY,SYZ,SXZ,TEMP 30.780 -5.9900 7.2000 0.0000 (TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 2SX,SY,SZ -19.180 -381.04 -24.140 SXY,SYZ,SXZ,TEMP -30.780 5.9900 -7.2000 0.0000 疲劳评定PERFORM FATIGUE CALCULATION AT LOCATION 1 NODE 0* POST1 FATIGUE CALCULATION *LOCATION 1 NODE 1389 EVENT/LOADS 1 1 AND 1 2 PRODUCE ALTERNATING SI (SALT) = 371.94 WITH TEMP = 0.0000 CYCLES USED/ALLOWED = 0.1000E+05/ 3445. = PARTIAL USAGE = 2.90242CUMULATIVE FATIGUE USAGE = 2.90242即采用 ansys 疲劳模块计算得到的疲劳次数为 3445 次7.3.2 考虑温度时的疲劳计算与评定考虑温度时容器在第 0.2 小时的温度分布云图如图 22 所示。第 1 小时、第 3 小时、第 5 小时、第 7 小时对应的应力云图分别如图 2326 所示。图 22 第 0.2 小时对应的容器温度云图图 23 第 1 小时对应的容器应力云图图 24 第 3 小时对应的容器应力云图图 25 第 5 小时对应的容器应力云图图 26 第 7 小时对应的容器应力云图以节点 1389 为疲劳评定点，相关应力数据如表 8 所示：表 8应力分量 任意时刻应力 (MPa) 最大或最小应力 (MPa) 应力差 (MPa)S11 18.72 -21.53 40.25S22 376.18 -405.36 781.54S33 23.24 -28.61 51.85S12 30.33 -32.99 63.32S13 6.96E+00 -8.41 15.37S23 5.95 6.2 -0.25相应的三个主应力值为 786.91，60.89，25.84(MPa)最大主应力波动差值为 786.91-25.84=761.07MPa采用 JB4732 疲劳载荷计算公式：对应的应力循环次数约为 3190 次。采用 ansys 的疲劳模块，得到：评定点在事件中各载荷点的应力值：LIST FATIGUE STRESSES FOR LOCATIONS 1 TO 5 IN STEPS OF 1EVENTS 1 TO 3 - LOADS 1 TO 10* POST1 FATIGUE STRESSES *LOCATION 1 NODE 1389 EVENT 1 (TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 1SX,SY,SZ 18.720 376.18 23.240 SXY,SYZ,SXZ,TEMP 30.330 5.9500 6.9600 0.0000 (TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 2SX,SY,SZ -21.530 -405.36 -28.610 SXY,SYZ,SXZ,TEMP -32.990 6.2000 -8.4100 0.0000 疲劳评定PERFORM FATIGUE CALCULATION AT LOCATION 1 NODE 0* POST1 FATIGUE CALCULATION *LOCATION 1 NODE 1389 EVENT/LOADS 1 1 AND 1 2 PRODUCE ALTERNATING SI (SALT) = 380.54 WITH TEMP = 0.0000 CYCLES USED/ALLOWED = 0.1000E+05/ 3203. = PARTIAL USAGE = 3.12215CUMULATIVE FATIGUE USAGE = 3.12215即采用 ansys 疲劳模块计算得到的疲劳次数为 3203 次7.4 采用 ansys 软件对容器进行的疲劳计算与评定分析方法：采用 ansys 直接软件 12.0 进行分析。考虑温度时采用直接耦合法。考虑材料非线性、几何非线性效应。采用固定时间步长以 0.2 小时为单位增量单位。对容器内壁施加随时间变化的温度载荷，外壁保持温度恒定为 20，且仅考虑热传导效应。有限元模型：单元为 solid5。分析单元为 6993 个，通过收敛性验证。网格模型如图 27 所示：图 27 容器有限元网格图网格质量检查：SUMMARIZE SHAPE TESTING FOR ALL SELECTED ELEMENTS-| Element count 6993 SOLID5 |-Test Number tested Warning count Error count Warn+Err %- - - - -Jacobian Ratio 6993 0 0 0.00 %Warping Factor 6993 0 0 0.00 %Any 6993 0 0 0.00 %-边界条件：边界条件如图 28 所示：图 28 容器边界条件示意图7.4.1 不考虑温度时的容器疲劳评定不考虑温度时容器在第 1 小时和第 3 小时的应力云图如图 29 所示：图 29 第 1 或第 3 小时对应的容器应力云图评定点在事件中各载荷点的应力值：LIST FATIGUE STRESSES FOR LOCATIONS 1 TO 5 IN STEPS OF 1EVENTS 1 TO 5 - LOADS 1 TO 10* POST1 FATIGUE STRESSES *LOCATION 1 NODE 4792 EVENT 1 (TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 1SX,SY,SZ 0.98587E+07 -0.58100E+07 0.35948E+09SXY,SYZ,SXZ,TEMP 0.21767E+07 0.22379E+06 -0.40174E+08 0.78886E-30(TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 2SX,SY,SZ 0.20798E-04 -0.95407E-06 0.35215E-05SXY,SYZ,SXZ,TEMP 0.31515E-05 -0.42409E-05 0.63062E-05 0.78886E-30(TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 3SX,SY,SZ -0.98585E+07 0.58101E+07 -0.35948E+09SXY,SYZ,SXZ,TEMP -0.21767E+07 -0.22376E+06 0.40174E+08 0.78886E-30(TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 4SX,SY,SZ -0.20616E-04 0.26279E-04 0.22150E-05SXY,SYZ,SXZ,TEMP -0.31923E-05 -0.43002E-05 -0.68344E-05 0.78886E-30(TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 5SX,SY,SZ -0.98585E+07 0.58101E+07 -0.35948E+09SXY,SYZ,SXZ,TEMP -0.21767E+07 -0.22376E+06 0.40174E+08 0.78886E-30(TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 6SX,SY,SZ -0.84352E-05 0.51755E-05 0.52036E-05SXY,SYZ,SXZ,TEMP -0.35773E-05 0.49984E-05 -0.82783E-05 0.78886E-30(TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 7SX,SY,SZ 0.98586E+07 -0.58100E+07 0.35948E+09SXY,SYZ,SXZ,TEMP 0.21766E+07 0.22379E+06 -0.40174E+08 0.78886E-30(TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 8SX,SY,SZ 0.19627E-04 -0.12463E-04 0.13471E-05SXY,SYZ,SXZ,TEMP 0.76218E-05 -0.51002E-07 0.58006E-05 0.78886E-30疲劳评定：PERFORM FATIGUE CALCULATION AT LOCATION 1 NODE 0* POST1 FATIGUE CALCULATION *LOCATION 1 NODE 4792 EVENT/LOADS 1 3 AND 1 7 PRODUCE ALTERNATING SI (SALT) = 0.37026E+09 WITH TEMP = 0.78886E-30CYCLES USED/ALLOWED = 0.1000E+05/ 3495. = PARTIAL USAGE = 2.86082CUMULATIVE FATIGUE USAGE = 2.860827.4.2 考虑温度时的容器疲劳评定考虑温度时容器在第 0.2 小时的温度分布云图如图 30 所示。第 1 小时、第 3 小时、第 5 小时、第 7 小时对应的应力云图分别如图 3134 所示。图 30 第 0.2 小时对应的容器温度云图图 31 第 1 小时对应的容器应力云图图 32 第 3 小时对应的容器应力云图图 33 第 5 小时对应的容器应力云图图 34 第 7 小时对应的容器应力云图评定点在事件中各载荷点的应力值：LOCATION 1 NODE 4792 EVENT 1 (TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 1SX,SY,SZ 0.12759E+08 -0.95293E+07 0.33127E+09SXY,SYZ,SXZ,TEMP -0.23766E+07 0.10794E+06 -0.35301E+08 177.50 (TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 2SX,SY,SZ 0.24866E+07 -0.31887E+07 -0.24179E+08SXY,SYZ,SXZ,TEMP -0.39028E+07 -99587. 0.41768E+07 155.00 (TOTAL STRESS - ITEMS 1-7) LOCATION 1 EVENT 1 LOAD 3SX,SY,SZ -0.77863E+07 0.31533E+07 -0.37963E+09SXY,SYZ,SXZ,TEMP -0.54289E+07 -0.30631E+06 0.43655E+08 132.50