非线性数值分析

非线性数值分析（精选7篇）

非线性数值分析 篇1

1 引言

压力容器已经广泛应用于化工、医疗、航天、航空、核工程等众多领域。目前, 压力容器正在向大型化和多参数化的方向发展, 尤其是在石油化工行业。与此同时, 压力容器也面临着高温、高压等更加复杂的环境。这就意味着对压力容器的设计有了新的要求, 传统的圆形压力容器已经不能满足要求, 对于非圆形压力容器的研究已经广泛展开。

试验应力方法、解析方法和数值方法是目前应力分析的主要方法。随着科技的发展, 出现了新的应力测定法方法, 比如光弹法、电测法、散斑法和全息干涉法等。其中, 数值分析方法已经广泛应用于工程之中。韩敏[2]利用ANSYS有限元软件对压力容器进行了应力分析, 分析结果与真实情况基本一致。马亚娟[3]用有限元的方法对压力容器的封头进行了应力分析, 其数值分析与试验的误差在10%左右。

本文将采取数值分析与试验相结合的方法, 利用有限元软件ANSYS对椭圆形压力容器进行应力分析, 并用实验数据来验证数值分析的正确性。从而大大节约了压力容器的设计成本和设计周期。并且在线性分析的基础上进行非线性数值分析, 数值分析与试验的误差在6%左右, 取得了较好的一致性。

2 非线性数值分析由于制造工艺的限制, 压力容器在制造过程中总会存在一定的缺陷, 而这种缺陷在有限元分析的建模中是难以实现的, 这就会引起数值分析与实际情况的误差。

非线性数值分析是在线性分析的基础上, 引入缺陷度的概念, 对有限元模型施加合适的初始缺陷, 再进行数值分析。缺陷度定义为设计最大误差与线弹性分析节点位移最大值的比值。

非线性数值分析的流程是在线性数值分析的基础上进行的。为了更详细的介绍非线性数值分析的流程, 将其所有的步骤展开如下:

第一步:在UG中建模。

第二步:将建好的模型导入ANSYS中。

第三部:对模型施加边界条件, 包括约束和载荷, 需要注意载荷的方向。

第四步:进行网格划分, 需要设置合适的网格尺寸, 尺寸太大会引起求解精度的下降, 尺寸太小会引起求解速度的下降。

第五步:进行线弹性求解。

第六步:通过后处理, 得到模型节点位移的最大值, 并计算缺陷度。

第七步:对模型施加合适的初始缺陷。

第八步:再次进行线弹性分析, 并观察结果。

3 有限元方法的验证

马亚娟利用ANSYS对球形封头及筒体进行了应力分析, 本文将采用非线性数值分析的方法对文献[3]中的封头及筒体进行分析, 以验证该方法的准确性。

在A N S Y S中建立球形封头及筒体的有限元模型及边界条件。筒体内径154.5mm, 壁厚8mm, 封头内径154.5mm, 壁厚8.1mm, A端施加固定约束, 限制其Y方向的位移, B端施加对称约束, 使其关于X轴对称, 模型内部施加1MPa的均布压力。

在模型的指定位置设置8个硬点, 以观察这些位置的应力。硬点的坐标分别是:0, 158.5, 0;40, 153.4, 0;80, 136.8, 0;110, 114.1, 0;130, 90.7, 0;150, 51.2, 0;156, 28, 0;158.5, -20, 0。经过计算, 通过后处理命令观察指定硬点的周向和径向应力。

为了说明分析的准确性, 将有限元分析的结果与测量值进行对比, 见表1。其中, ‘σθ实’表示实测周向应力, ‘σθ非’表示非线性分析的周向应力, ‘σθ线’表示线性分析的周向应力, ‘σr实’表示实测径向应力, ‘σr非’表示非线性分析的径向应力, ‘σr线’表示线性分析的径向应力, 从上到下依次对比8个硬点的应力。

从表1中可以看出, 除了y=51.2处的硬点, 对于其他的硬点, 有限元的计算结果与测量值吻合性较好, 相对误差在6%左右, 能够准确的进行压力容器的应力计算。而且非线性分析的结果比线性分析的结果偏小, 更加接近真实值, 这都是在分析中考虑了压力容器的初始几何缺陷的缘故。

4 椭圆断面压力容器数值分析

椭圆断面压力容器具有与圆形断面压力容器不同的特点和受力状态, 有必要对其进行研究。本文将采用有限元的方法对椭圆断面压力容器进行应力分析, 对于非圆形压力容器的进一步研究有一定的参考价值。

椭圆断面压力容器并非真正意义上的椭圆, 只是近似的椭圆形。其由半径是R的大圆弧和半径是r的小圆弧组成, 大小圆弧的厚度分别是δ1和δ2, 其值可以相等也可以不相等, 大小圆弧的圆心角分别是2θ和2φ。

本文中设置椭圆断面压力容器的具体参数如下:R=150mm, r=30mm, θ=30°, φ=60°, δ1=δ2=8mm。由于椭圆形压力容器是对称结构, 本文在分析中将建立模型的1/4部分, 然后施加边界条件。对模型施加对称约束, 然后对模型内部施加1MPa的均布压力。

通过计算, 得到了椭圆断面压力容器的周向和径向应力云图, 如图1所示。从图1左半部分可以看出, 压力容器大圆弧中部外侧的周向应力最大, 小圆弧中部外侧的周向应力最小, 表明椭圆形压力容器的周向应力分布是不均匀的。观察图1右半部分可以看出, 径向应力分布方式与周向应力的分布有所不同, 径向应力最大值出现在压力容器大圆弧中部外侧和小圆弧中部内侧, 最小值出现在小圆弧中部外侧, 径向应力的分布同样是不均匀的。

结语

椭圆形压力容器的应力分布与圆形压力容器完全不同, 其沿着周向的应力大小各不相同, 这对于非圆形压力容器的设计有一定的参考价值。在非圆形压力容器的设计中, 根据有限元分析的结果, 将应力较大的部分的壁厚加厚, 将应力较小部分的壁厚适当减薄, 使其应力的分配尽量均布化。如此设计, 既可以保证非圆形压力容器能够承受更大的压力, 又可以保证压力容器不会太笨重。

非线性数值分析考虑了压力容器在制造过程中的几何缺陷, 使得数值分析的结果与真实值更加接近, 有限元分析的可靠性得到了一定的提升。

但是在分析中, 还会出现数值分析与实验结果不一致的地方, 这是由于有限元模型与试样的不完全一致造成的, 这个问题将是本文作者今后研究工作的重点。

摘要：本文采用有限元的方法对于椭圆断面压力容器进行了非线性数值分析。并用一个算例验证了非线性数值分析的可靠性, 结果表明数值分析与实测结果的一致性较好, 有限元分析完全可以满足工程需要。最后, 通过分析得出的结论是椭圆形压力容器的应力分布是不均匀的。这个结论对于非圆形压力容器的设计有较大的参考价值和指导作用。

关键词：椭圆断面,压力容器,非线性,有限元,几何缺陷

参考文献

[1]Liu P, Xu P, Han S, et al.Optimal design of pressure vessel using an improved genetic algorithm[J].Journal of Zhejiang University SCIENCE A, 2008, 9 (09) .

[2]韩敏.利用ANSYS软件对压力容器进行应力分析[J].煤矿机械, 2008, 29 (01) :73-74.

[3]马亚娟.压力容器的球形封头和椭圆形封头的应力测定及分析[D].西安工程大学, 2014.

非线性数值分析 篇2

以一种新型的.Boussinesq型方程为控制方程组,采用五阶Runge-Kutta-England格式离散时间积分,采用七点差分格式离散空间导数,并通过采用恰当的出流边界条件,从而建立了非线性波传播的新型数值模拟模型.通过对均匀水深水域内波浪传播的数值模拟说明,模型能较好地模拟大水深水域和强非线性波的传播.通过设置不同的入射波参数来进行潜堤地形上波浪传播的物理模型实验,并将数值解与物理模型实验结果进行了比较.

作 者：张洪生 冯文静 王亚玲 吴中 缪国平ZHANG Hong-sheng FENG Wen-jing WANG Ya-ling WU Zhong MIAO Guo-ping 作者单位：张洪生,冯文静,缪国平,ZHANG Hong-sheng,FENG Wen-jing,MIAO Guo-ping(上海交通大学,船舶海洋与建筑工程学院,上海,30)

王亚玲,WANG Ya-ling(上海海事大学,交通运输学院,上海,35)

吴中,WU Zhong(河海大学,交通与海洋工程学院,江苏,南京,210024)

非线性数值分析 篇3

HRB500级钢筋作为GB 50010—2010《混凝土结构设计规范》[1]推荐采用的受力主筋,具有强度高,延性好、降低材料消耗等优点,但对于提高构件的变形能力的作用并不明显,可能导致结构出现过大的裂缝从而不能满足结构的适用性要求[2]。而聚丙烯纤维作为一种增强材料,可以有效增强混凝土抗冲击性及柔韧性,提高普通混凝土力学性能[3,4]。可见,在高强钢筋混凝土中掺入适量聚丙烯纤维成为一种提高强钢筋混凝土构件抗裂性能的科研思路。同时,现有文献中虽然也有一些聚丙烯纤维混凝土数值模拟的分析[5,6,7],但是对配置HRB500级钢筋的聚丙烯纤维混凝土构件的受力性能的数值模拟分析尚未发现。因此,本文通过配置高强钢筋的聚丙烯纤维混凝土梁的受弯试验来分析聚丙烯纤维对混凝土抗裂性能的影响,并在试验分析基础上利用有限元软件ANSYS模拟聚丙烯纤维混凝土梁的工作性能。

1 试验概况

对6根不同聚丙烯纤维掺入量的配置HRB500级钢筋的混凝土梁进行了受弯性能试验,构件分为L、LX-A及LX-B三组,每组两根混凝土梁。其中,L系列的混凝土梁未掺入聚丙烯纤维,LX-A、LX-B系列为掺入定量聚丙烯纤维的混凝土梁,其掺入量分别为0.9kg/m3和1.2kg/m3。本试验皆为破坏性试验,采用两点对称集中加载方式。试验梁的截面形式如图1所示,配筋情况如表1所示。聚丙烯纤维的参数如表2所示,混凝土及HRB500级钢筋的相关力学性能如表3、表4所示。

注:fcu,100—混凝土立方体100mm×100mm×100mm抗压强度;fcu—混凝土标准立方体150mm×150mm×150mm抗压强度,fc—混凝土轴心抗压强度;ft—混凝土抗拉强度。

2 试验结果分析

2.1 试验现象

在试验施加荷载作用下,L系列在梁的跨中纯弯段范围内首先出现竖直裂缝,随着试验荷载的不断增加,裂缝数量逐渐增多,裂缝宽度也在增大,但是裂缝间距却不断减小,当试验加载逐渐接近试验梁的极限荷载时,试验梁的裂缝数量趋于稳定,裂缝间距不再变化,但是裂缝的宽度却增长较快,尤其是纯弯段内的主裂缝宽度增加最快,直到试验结束。掺入定量聚丙烯纤维的LX-A、LX-B系列的试件裂缝发展规律与L系列试件基本一致,但是同L系列试件相比,LXA、LXB系列试件的裂缝数量更多,平均间距也明显减小,如图2所示。由此可知,掺入聚丙烯纤维的试验梁裂缝数量更多,通过更多的裂缝来承担试验梁开裂的能量,从而提高了试件的抗裂能力。

2.2 裂缝宽度及限值分析

GB 50010—2010在正常使用状态验算时规定:“三级裂缝控制等级时,钢筋混凝土构件的最大裂缝宽度可按荷载准永久组合并考虑长期作用影响的效应计算,三级裂缝控制等级是允许出现裂缝的构件”[1],即文中试验梁均属于三级裂缝控制构件,其最大裂缝宽度不应超过最大裂缝宽度限值0.3mm。依据试验实测的短期裂缝宽度乘以长期荷载作用下的裂缝增大系数1.5,得到推算的使用阶段裂缝宽度,并与0.3mm的限值进行比较,结果如表5所示。

从试验计算结果表5可以得出,LX1-A裂缝宽度为L1裂缝宽度的87.12%,LX1-B裂缝宽度为L1裂缝宽度的87.43%,LX2-A裂缝宽度为L2裂缝宽度的81.59%,LX2-B裂缝宽度为L2裂缝宽度的81.78%。即掺入0.9kg/m3及1.2kg/m3的聚丙烯纤维可使试验梁的正常使用阶段裂缝宽度减小,且减小数值相差不多。由此可见,聚丙烯纤维的掺入能够有效的减小高强钢筋混凝土梁的裂缝宽度,从而提高混凝土梁的抗裂性能。同时,试验梁在使用阶段裂缝宽度值都小于0.3mm,按GB 50010—2010《混凝土结构设计规范》规定,配置HRB500级钢筋的聚丙烯纤维混凝土受弯梁在使用阶段能够满足裂缝宽度限值的要求。

注:ωexps,max—短期最大裂缝宽度的实测值;ωexpl,max—推算的长期最大裂缝宽度。

3 数值模拟分析

3.1 模型建立

结合ANSYS中钢筋混凝土整体式和分离式的两种建模思路[8],把聚丙烯纤维作为一种类似钢筋的加强材料和混凝土整体式建模,钢筋和混凝土还是以分离式建模[9]。混凝土采用钢筋混凝土实体单元Solid65,钢筋采用Link8单元,聚丙烯纤维按掺量在混凝土实常数设置中输入三个方向的体积率。依据试验情况,为避免应力集中现象,在两个支座处设置刚性垫块,采用Solid 45单元,混凝土破坏准则采用Willianop-warnke 5参数强度准则,其参数设置参照相关文献[8,10],裂缝处理方式为弥散固定裂缝模型。混凝土及钢筋的计算模型如图3所示。

3.2 本构关系

(1)混凝土的本构关系采用GB50010—2010给定的曲线,采用ANSYS中的多线性各向同性强化(MISO)模型。

(2)钢筋的本构关系按理想弹塑性考虑,采用ANSYS中的双线性各向同性强化(BISO)模型[11]。

(3)聚丙烯纤维采用理想线性模型,采用厂家提供的弹性模量和泊松比。

3.3 计算结果对比

(1)试验梁极限荷载有限元数值模拟结果和试验结果对比如表6所示。

由表6可以看出,6根试验梁的极限荷载的模拟结果与试验值相差很小,最大差值仅为5.3%,模拟效果较好。

(2)试验梁有限元数值模拟裂缝开展图如图4所示。

由图4可以看出,其裂缝分布规律与分布范围与试验现象(图2)相符,裂缝首先出现于纯弯段,加载中后期,混凝土梁弯剪区段陆续出现裂缝,可见裂缝分布模拟结果较好。

(3)试验测得跨中纯弯段纵筋应变发展情况,作荷载-钢筋应力曲线对比,如图5所示。

由试验梁LX1-B的纵筋荷载—应力曲线可以看出,模拟值与试验值吻合较好,当荷载大于60k N·m时,模拟值与试验值之间差值逐渐增大,这主要是钢筋建模时采用理想弹塑性模型的结果。

(4)试验梁荷载-挠度曲线对比

试验梁LX1-B荷载-挠度曲线如图6所示,ANSYS计算结果与试验实测值曲线发展规律一致,差值在10%以内,最大差值出现在荷载加载20~40k N·m的范围内。主要原因是在这个荷载范围内混凝土开裂,而在建立有限元模型时,没有考虑各种材料之间的粘结力,导致在刚开裂的这段区域计算值与实测值产生相对较大的差值。

通过对试件的极限荷载、裂缝开展图、跨中纯弯段纵筋的荷载-应力曲线以及荷载-挠度曲线的数值模拟值与试验值的对比,发现模拟值与实测值吻合较好,说明采用本文的有限元模型和建模方法是正确的,能够有效模拟高强钢筋聚丙烯纤维混凝土梁的受力过程。

4 结论

(1)聚丙烯纤维混凝土受弯梁和未掺入聚丙烯纤维的混凝土受弯梁裂缝发展规律相同,但掺入聚丙烯纤维的混凝土梁裂缝细而密,裂缝宽度较同条件下未掺入聚丙烯纤维的混凝土梁小,掺入聚丙烯纤维能够有效提高混凝土梁的抗裂性能。

(2)配置HRB500级钢筋的聚丙烯纤维混凝土梁在正常使用极限状态能够满足GB 50010—2010《混凝土结构设计规范》裂缝宽度限值的要求。

(3)采用ANSYS有限元软件对聚丙烯纤维混凝土进行非线性数值模拟,计算结果与试验值较为接近,差值在10%以内,该方法能够有效的模拟聚丙烯纤维混凝土梁的受力性能。

参考文献

[1]GB50010-2010混凝土结构设计规范[S].北京:中国建筑工业出版社,2010.

[2]李艳艳.配置500MPa钢筋的混凝土梁受力性能的试验研究[D].天津:天津大学,2008.

[3]孙家瑛.聚丙烯纤维对高性能混凝土抗折强度、抗冲击性能影响研究[J].混凝土,1999(6):19-21.

[4]朱江,苏健波.聚丙烯纤维混凝土的力学性能研究[J].广西工学院学报,2000,11(2):60-64.

[5]王金晶,刘志奇.透水性聚丙烯纤维混凝土的耐久性研究[J].混凝土,2009(11):31-33.

[6]隋红军,贾宝新.聚丙烯纤维混凝土锚喷结构内力有限元分析[J].科学技术与工程,2007,7(15):160-165.

[7]吴永详.马洪旺.聚丙烯纤维混凝土抗拉增韧性能的初探及其在地铁各工程中的试用[J].四川建筑科学研究,2009,35(5):167-170.

[8]江见鲸,陆新征.钢筋混凝土结果有限元分析[M].北京:清华大学出版社,1998.

[9]胡飞玲,聚丙烯纤维混凝土梁的试验分析[D].浙江:浙江大学,2004.

[10]A.F.Ashour,C.T.Morley,Three-dimensional nonlinearfinite element modeling of reinforced concrete structures[J].Finite Elements in Analysis and Design,1993,15(1):43-55

覆冰单导线舞动非线性数值模拟 篇4

基于上述原因, 本文采用考虑扭转自由度的两节点索单元模拟覆冰导线舞动。首先推导索单元的参数矩阵, 给出了刚度矩阵的显式表达, 其次考虑导线运动过程中的几何非线性和动荷载与位移及速度的耦合作用, 最终得到了非线性动力学有限元方程。采用对加速度中心差分、对速度向后差分的时间积分法计算了覆冰单导线的舞动响应。

1 覆冰导线舞动的动力有限元方程

导线的覆冰舞动是典型的几何非线性问题, 其伸长变形一般认为在线弹性小变形范围［11，12］。更新拉格朗日格式的有限元方程为

式 ( 1) 中: M为质量矩阵;为t + Δt时刻节点的加速度; C为t时刻结构的阻尼矩阵;为t + Δt时刻节点的速度; Ktt+Δt为t至t + Δt时刻结构的刚度矩阵; qt +Δt为t + Δt时刻的节点位移增量向量;Ft +Δt为t + Δt时刻的荷载向量; Qt为t时刻的内部节点荷载向量。

1. 1 单元分析

由于覆冰舞动分析时需要考虑导线的扭转自由度, 现在构建两节点索单元时考虑了导线的3 个平动自由度u, v, w以及一个扭转自由度 θ。沿导线轴向的自然坐标为s, 如图1 所示。

以t时刻为基准进行分析, 从t时刻到t + Δt时刻的位移为ui, 扭转角为 θ, 单元节点位移向量qe由下式定义

式 ( 2) 中右下标表示单元中局部节点的编号, 从t时刻到t + Δt时刻单元上任意一点的位移和转角可以表示为:

形状函数矩阵N为:

形函数N1= s / L, N2= ( 1 - s) / L, I4为4 ×4 阶的单位矩阵, L是单元的长度。

单元质量矩阵采用一致质量矩阵可表示为

式 ( 4) 中 ρ 为导线的线密度; A为导线的截面面积;Jx为导线的扭转惯量。

考虑轴向扭转自由度, 索单元的应变增量 Δε和应力增量 Δσ 分别为

应变增量可表示为

式 ( 6) 中:

将式 ( 5) 按泰勒公式展开考虑前二阶, 并将应变增量分为线性部分和非线性部分可得:

单元的轴向应变矩阵:

得到单元应变矩阵: B = BL+ BNL。

整理可得用节点坐标表示的单元应变矩阵

H为系数矩阵:

从t时刻到t + Δt时刻单元的刚度矩阵为

式 ( 12) 中: σt +Δt为t + Δt时刻的单元应力向量; σt为t时刻的单元应力向量; D是弹性矩阵。

结构阻尼的确定比较困难, 对于导线可以采用Raleigh阻尼［13］,

式中: cij为结构阻尼矩阵对应位置的元素; mij为质量矩阵对应位置的元素; kij为刚度矩阵中对应位置的元素, ξi、ξj分别为对应于频率 ωi、ωj的阻尼比。

由于覆冰导线的截面随时间发生变化, 因此t时刻的空气动力荷载与t时刻的位移和速度是耦合的。在求解作用在导线上的空气动力荷载时, 应先求解t时刻的位移和速度, 在编写程序时做相应处理。空气动力荷载由下式确定

式中 ρair为空气密度; d为裸线直径; U为水平风速;CL为升力系数; CD为阻力系数; Cθ为扭转系数。

空气动力系数Ci, 其值取决于覆冰的形状, 一般需要风洞试验测定［14］, 通常将由风洞试验得到的气动力系数拟合为攻角 α 的3 次多项式近似:

由此可以计算单元节点外荷载向量:

内部节点荷载向量:

式 ( 16) 中T为导线t时刻的轴向拉力; M为导线t时刻的扭矩。

1. 2 邻档输电导线和绝缘子串的模拟

邻档输电导线和绝缘子串均采用线性弹簧模拟, 邻跨导线等效弹簧刚度为Kc1［15］, 绝缘子串在x和z方向的等效弹簧刚度为Kx和Ky［16］。

式中L0为无张力时导线原长; E为导线的弹性模量; Py为单位长度导线受到的竖直方向上的荷载;L1为档距; T为静止状态下张力的水平分量; L2为绝缘子串的长度; W1为绝缘子串的重量。

2 动力平衡方程求解

2. 1 递推公式推导

考虑空气动力荷载与位移和速度的耦合, 在数学表达上即为动荷载Ft是位移q和速度q·的函数。

对加速度采用中心差分, 对速度采用向后差分:

将式 ( 19) 代入式 ( 1) 得

整理上式可得采用差分法求解方程的递推公式:

Ft是qt和的函数, 式 ( 21) 中qt已知, 可以通过式 ( 20) 求得。qt和为已知后, 可以通过函数关系求得Ft。此时递推公式右部都是已知量, 递推可以进行下去。

在递推求解起步时, 初始节点速度为已知, 可由求出qt －Δt后即可通过式 ( 21) 递推求解。

2. 2 计算程序

利用MATLAB编制计算程序, 流程图如图2 所示。时间间隔 Δt = 0. 005 s, 模拟时间为2 000 s。

3 数值算例及分析

3. 1 与文献的计算结果进行对比

本算例利用文献[6]的导线舞动试验结果验证本方法的正确性。舞动实验中导线的覆冰用D型木块模拟, 档距为125. 88 m, 覆冰导线中点的弧垂为1. 38 m, 转角位移0. 1725 rad, 舞动时风速约为4. 1 m / s。绝缘子串长度为2. 1 m, 自重490 N, 其他物理参数如表1 所示。

根据本文给出的索单元, 通过验证, 采用20 个两节点索单元对本算例中的覆冰导线进行离散。

图3 和图4 分别为计算所得的覆冰导线的舞动轨迹和位移时间历程, 由此可以看出, 开始振动时导线在平衡位置附近做小幅度的摆动, 此时水平振幅较大, 随后水平振幅减小竖直振幅增加, 最后均趋于稳定, 稳定后轨迹近似为椭圆。

将计算结果与文献[6]所得的计算结果进行对比。在舞动稳定状态下, 文献中舞动轨迹中心位置为 ( 0. 04, 1. 38) , 本文中舞动轨迹中心位置为 ( 0. 034, 0. 053) , 这是由于文献算例在计算时将式 ( 14) 中空气动力系数常数项对应的力作为静力施加在模型上, 本文将其作为动力荷载施加在导线上;文献计算结果竖向最大振幅为0. 79 m, 水平最大振幅为0. 02 m, 本文结果分别为0. 793 m和0. 020 6m, 两者的相对差分别为0. 4% 和2. 9% 。结果表明, 对于覆冰单导线, 本文的方法及计算结果是合理可行的。

3. 2 导线在不同风速下的舞动分析

根据前面论述的理论和方法, 选择文献[6]中的算例, 考察风速对舞动的影响。

在初始攻角不变的情况下, 不同风速下单根导线的运动轨迹如图3 及图5 ~ 图8 中所示, 可以看出随着风速增大, 导线的振幅不断增大, 但导线竖向振幅增幅较大, 水平向振幅增加较小; 当风速增大到一定程度, 竖向振幅增加开始变得不明显。舞动是风能逐渐积累的过程, 在风速没有到达某个临界值时, 导线自身的衰减能抵消外界风速输入的能量, 但当风速太大导线自身的衰减已不能抵消外界输入的能量, 导线位移会越来越大。由于导线舞动的本质是非线性振动, 最终其舞动振幅会达到某个极限值处稳定。

4 结论

( 1) 采用两节点索单元模拟覆冰导线舞动, 并考虑了空气动力荷载与位移及速度的耦合作用, 运用对加速度中心差分、对速度向后差分的求解方法, 通过算例验证了本文舞动计算方法及程序的正确性, 且具有较高的计算精度及效率。

( 2) 导线舞动时竖向振幅大于水平振幅, 当风速较小时, 导线振幅越来越小最终趋于稳定; 当风速较大时, 导线的竖向振幅远远大于水平振幅。在初始攻角一定时, 导线舞动的竖向振幅及水平振幅都会随着风速增大而增大。

非线性数值分析 篇5

关键词：煤层气,非线性,拟稳态,COMSOL Multiphysics

煤层岩孔隙属于孔隙-裂缝介质, 尤其对于低渗透储层, 孔隙结构更复杂, 在裂隙界面与煤层气存在较强的吸附作用, 复杂作用的影响使煤层气在低渗透储层中渗流规律并不遵循达西渗流模型。本文通过对低渗透煤层的孔隙结构及气-固作用力方面进行研究, 建立煤层气非线性渗流模型, 并应用COMSOL Multiphysics软件进行数值模拟研究。

1 低渗透煤层孔隙特征及渗流模型引入

煤是一种具有复杂孔隙结构的聚合物质, 煤层气主要吸附在孔隙的内表面。煤层气与固体表面作用力的强弱决定了吸附量, 同时也影响着煤层气的渗流。因此, 煤层气-固界面作用对于致密低渗透煤层介质煤层气渗流的影响, 比中高渗透煤层中煤层气渗流的影响大得多。低渗透煤层气渗流遵循一种非线性渗流规律[3,4]。本文引入新的煤层气渗流模型:

2 低渗透煤层气低速非线性渗流数学模型

2.1 气体的控制方程

气体的物质控制方程为

其中

忽略重力的影响:

则

2.2 解吸-扩散方程

煤层气的扩散视为拟稳态, 遵循Fick第一定律, 单位体积煤体所吸附的气体质量为

2.3 特性参数

本文采用由Palmer和Mansoori提出的孔渗模型, P&M模型考虑基质的应变变化线弹性方程及煤层气解吸后煤基质收缩的影响。

根据煤储层孔渗透率和孔隙度的关系

2.4 煤层气拟稳态数学模型

将方程 (9) 代入方程 (5) , 则得到:

2.5 初始条件及边界条件

考虑煤储层中甲烷的吸附解吸过程, 设置模型初始条件和边界条件。

初始条件:

定压内边界:

封闭边界:

3 模型分析

模型尺寸为100m×100m, 井位于模型中心, 模型为封闭边界, 取模型的1/4进行数值模拟研究, 选用COMSOL Multiphysics建立模型。

考虑三种渗流规律模型计算365d后, 煤层孔隙压力分布如图2所示, 定产量生产时, 非线性渗流模型和启动压力梯度模型需要消耗更多的能量, 启动压力梯度模型由于悲剧地夸大了地层的阻力, 近井地带压力降落更大。达西渗流模型生产后, 整个气藏均动用, 压力降落均匀, 非线性渗流模型和启动压力梯度模型由于存在最小启动压力梯度和启动压力梯度, 存在未启动区, 即图2中压力水平段。

如图3所示了不同非线性系数影响下孔隙压力分布, 非线性渗流越严重, 气藏的压力降落越大, 同等的产能需要消耗的地层能量越大, 同时存在的非启动区越大, 因此, 考虑低渗透煤层气非线性渗流规律对研究低渗透煤层气储层的生产具有重要的意义。

4 结论

(1) 考虑低渗透煤层气孔隙结构及吸附特征, 建立非线性渗流模型, 模型具有普遍适用性

(2) 对比达西渗流模型和启动压力梯度模型, 非线性渗流模型更好的描述低渗透煤层气的渗流规律, 同等产能, 非线性渗流模型和启动压力梯度模型压力降落更严重, 消耗更多的能量。

(3) 研究了不同非线性系数影响下, 储层中压力分布规律及储层非启动区范围。

参考文献

[1]Zhang H B, Liu J and Elsworth D.Howsorption-induced matrix deformation affectsgas flow in coal seams:a new FE model[J].International Journal of Rock Mechanics andMining Sciences, 2008, 45 (8) :1226-1236

[2]邓英尔, 马宝岐, 刘慈群.煤层气—固界面作用与吸附模型[J].煤田地质与勘探, 2003, 32 (3) :24-26

[3]张先敏, 冯其红, 陈东等.阜新盆地煤层气渗流规律实验[J].重庆大学学报, 2011, 34 (4) :58-61

[4]邓英尔, 谢和平, 黄润秋等.低渗透孔隙-裂隙介质气体非线性渗流运动方程, 四川大学学报, 2006, 38 (4) :1-4

非线性数值分析 篇6

实验及分析研究[1]表明,钢筋混凝土双向受弯构件在受荷过程中,随着荷载的增大,挠度的不断增大,截面形心将逐渐偏离荷载作用平面,从而不可避免产生次生扭矩,且次生扭矩的大小随着荷载的增加和构件裂缝产生及破坏的发展不断变化,并对双向受弯构件的受力性能产生一定影响。两点对称加载时,扭转角很小,当发生正截面破坏时,扭转效应可忽略不计。文中通过Matlab编写电算程序,模拟了双向受弯的钢筋混凝土梁非线性破坏过程。

1 基本假定

1)平截面假定:梁受弯以后,截面混凝土、钢筋应变分布符合平截面假定。钢筋混凝土构件开裂前能近似满足这一假定,构件开裂后,裂缝处截面显然不能满足平截面假定。但是国内外的大量试验表明,若采用大标距应变计测得标距范围内的平均应变,则从加载到构件破坏全过程的平均应变值都能很好的满足平截面假定。

2)钢筋和混凝土之间粘结可靠,无相对滑移。

3)忽略剪切变形、忽略次生扭矩对构件变形的影响。

2 计算公式

1)应变、应力以受拉为正、受压为负,力以拉力为正、压力为负。弯矩、曲率均为正。长度单位mm,力单位N。

2)钢筋应力—应变关系:

取:εu=0.01,ε0=0.002;

-ε0≤ε≤ε0,σs=Esε;

εu≥ε≥ε0或-εu≤ε≤-ε0,σs=σsmax;

ε>εu或ε<-εu,σs=0。

3)混凝土应力—应变关系:

取:εcu=-0.003,ε0=-0.002,εcu=ftk/Ec。其中,ftk为混凝土最大拉应力。

ε≤εu,σ=0.3σmax;

ε0>ε>εu,σ=σmax{1-[200(ε-ε0)]2};

0≥ε≥ε0,σ=2σmaxε/(εu+ε0);

εcu2≥ε≥0,σ=Ecε;

ε≥εcu2,σ=0。

当梁角处受拉钢筋屈服后,不再考虑混凝土的受拉承载力。

4)当梁角处的钢筋屈服时,开始考虑塑性铰。但不考虑塑性铰的发展过程。

按下式计算:

lp=2(A+B)。

其中,A为受压作用点到中和轴距离;B为最远端钢筋到中和轴距离;C为受压合力作用点。

5)截面采用网格划分。

对混凝土,截面宽度b划分n,高度h划分m,分块后,计算信息可以用m×n的矩阵存储。对钢筋,计算信息采用p×q的矩阵进行存储。

6)已知M—ϕ关系曲线,将梁沿长度方向分成n等份,给定荷载P,可得到梁的弯矩分布图。根据梁中任一节点i处的弯矩Mi,由截面的M—ϕ关系可确定该节点处的曲率ϕi,于是节点i处的转角$\theta {}_i=\frac{1}{2}{\displaystyle \sum _{j=1}^i(}\text{ϕ}{}_j+\text{ϕ}{}_{j+1})\text{Δ}X$,挠度$\delta {}_i=\frac{1}{2}{\displaystyle \sum _{j=1}^i(}\theta {}_j+\theta {}_{j+1})\text{Δ}X$。

3 程序流程图[2]

3.1M—ϕ关系程序流程图(见图1)

3.2P—Δ关系程序流程图(见图2)

ϕm为梁跨中曲率;Mm为梁跨中弯矩;δ1i为位移曲线1;δ2i为位移曲线2;δ1m为曲线1梁跨中位移;δ2m为曲线2梁跨中位移。

4 对比试验[3]

试验对比选取文献[3]中的两组数据,见表1。

拉力筋为直径25 mm,其余三根均为12 mm。整个梁长2 m。

4.1 试验结果

弯矩曲率关系原文未给出:

PAS-18-L1的破坏荷载为F=7.04 t;PAS-18-L3的破坏荷载为F=8.15 t。

最大位移值可以从文中的图中通过比例尺测得,测得的最大位移值约为25 mm。

4.2 程序得到结果(见图3)

4.3 结果对比分析

1)符合情况与试验结果有一定偏差,文中提到了支托钢筋滑移引起的附加转角对构件变形影响很大,约占极限挠度的20%,而程序并没有考虑到这一影响。这可能是产生偏差的原因。2)文中采取的是不对称配筋,截面的扭转在荷载作用角很小时,很明显,通过程序计算,当输入角度大于10/π时,就必须考虑扭转对荷载角的修正。3)塑性铰的计算按照文中最初的计算方法算得的分别为240 mm,233 mm,试验结果为358 mm,307 mm。

参考文献

[1]吴光宇,项贻强,徐兴.钢筋混凝土双向受弯梁次生扭矩的分析研究[J].工程力学,2008,22(4):17-18.

[2]王沫然.MATLAB与科学计算[M].北京:电子工业出版社,2000.21-22.

[3]陈忠汉,朱伯龙,钮宏.斜向受力钢筋混凝土压弯构件的非线性分析[J].土木工程学报,1984(4):24-25.

非线性数值分析 篇7

1 基本方程

1.1 激发方程

令第次本征模式轴向电场为

可得扰动电子注激励的电场为[5,6,7]:

式中:定义Kn=为n次谐波的耦合阻抗;(r⊥)为电子注横向分布函数ds。所有本征模中只有个别模式与电子注同步,且除了电子流i(z)激发的同步波之外,还有输入的“冷波”,即,则具有外加激励源的同步场为:

式中:K0,Γ0,β0分别为该同步模式的耦合阻抗、传播常数、相位常数。式(1)两边同时对z求导2次得到熟知的激发方程:

1.2 运动方程

相对论下电子的运动方程为:

能量守恒方程为:

式(4)代入式(3)可得一维电子运动方程:

又由:

所以式(5)可写为:

其中:Ez=Ecz(线路场)+Esz(空间电荷场);γ=(1-v2/c2)-1/2为相对论因子,c为真空中的光速,v为电子速度。

1.3 电子流复振幅方程

电子流是时间的非简谐周期函数,含有高次谐波,用傅氏分析。

2 慢变系统中归一化

上述是在实验室坐标系下得到的迅变方程,为了处理问题的方便和计算结果普遍性强,通常将其归一化到电子坐标系内,获得慢变方程。

为了表述方便,先引入迅变坐标系的归一化量:归一化距离为ξ=ω/v0z=βez;归一化时间为φ=ωt=2πt/T,归一化场为f=|e|E/mv0ω。则慢变系统中的归一化:归一化轴向距离为θ=Cξ=Cβez;归一化相位Υe=φ-ξ;归一化场幅值为Fcn(θ)=;归一化电流幅值为。

激发方程变为:

运动方程变为:

电子流复振幅方程为:

式(8)～(10)组成了行波管大信号注-波互作用基本工作方程组。其中=I0Kcn/4V0为n次谐波增益参量,非同步参量,dn=α0n/βeCn为衰减常数,rn=bn-idn。

3 空间电荷场的计算

由文献[8]可知z0处圆盘在z处圆盘平面上各点产生的平均空间电荷场为:

其中:Q为圆盘所带电量;b为电子注半径;a为漂移管半径,如图1所示,μ0n为零阶Bessel函数的第n个根。由此可知场是关于z的函数,可以表示为:

其中:B(|z-z0|)是以|z-z0|为变量的函数,由式(11)可以做出如图1所示曲线。

由图1可知,如果用近似式Esz=QB(|z-z0|)代替式(11),所引起的误差很小,而计算式却大大简化了。采用该式时,所有圆盘产生的空间电后为:

其中:ωp为电子等离子体频率;η为电子荷质比;v0为电子注初速度。仿真中用式(13)代替式(9)中的Fsz。

4 休斯结构耦合腔

休斯型耦合腔的剖面结构为图2(a)所示,图2(b)为其截面图,表1为设计S波段耦合腔行波管的结构参数。

5 仿真结果及讨论

在上述尺寸下,对2.87～3.35 GHz频率范围内的休斯结构耦合腔行波管进行了数值分析,仿真中电子注的注电压U0=21 kV,注电流I0=1 A,波的输入功率Pin=300 W。采用四阶龙格-库塔法[9]求解互作用方程组式(8)～(10),其结果如图3～图6所示。

图给出了在中心频率时电子的相位轨迹。由图可知,电子在归一化轴向距离Z=2.4附近获得了较好的群聚,此即为最佳互作用距离。由图4可以看出该频率的波在此处获得最大增益,此后电子注离开最佳互作用区,效率降低,增益下降。

图5给出了中心频率时电子注效率随轴向距离的变化曲线。图中实线为不考虑空间电荷力的情况,虚线为考虑空间电荷力的情况。由图可知,空间电荷力的作用使得饱和位置推后,增益下降,即空间电荷力起发散作用。图中效率刚开始时为负值,是因为电子刚进入互作用区时要得到部分能量,表现电子效率为负值。

图6给出考虑空间电荷场时,2.87～3.35 GHz内各频点的微波增益。由图可知,在给定电子注注电压,注电流,和波的输入功率等参量的情况下,频带内微波增益均大于18.5 dB。

6 结语

对S波段休斯结构耦合腔行波管非线性注-波互作用工作方程组进行了数值求解,求出考虑和不考虑空间电荷场时中心频率f=3.100 7 GHz处的效率,说明空间电荷场对互作用起散焦作用,与文献[2,10]中结论很吻合。求出工作在2.87～3.35 GHz频率范围内耦合腔行波瞬时带宽。仿真中所用管子已制作完成,实验数据对后期管子的热测试有很好的指导意义。

参考文献

[1]李文君,许州,林郁正,等.S波段宽带大功率连续波耦合腔行波管三维模拟设计[J].强激光与粒子束,2006,18(6):965-968.

[2]徐林,杨中海,莫元龙.行波管三维非线性计算机模拟的改进[J].强激光与粒子束,1999,11(3):342-346.

[3]雷文强,杨中海,廖平.耦合腔慢波特性及其注-波互作用的三维软件模拟[J].强激光与粒子束,2003,15(7):673-676.

[4]Kantrowitx F,Tammaru I.Three-dimensional Simulations ofFrequency-phase Measurements of Arbitrary Coupled-cavityRF Circuit[J].IEEE Trans.on Electron Devices,1988,35(11):2 018-2 026.

[5]Vaninstin L A.Waveguide Excited by the Electron Beam[J].Journal of Technology Physics,1953(4):654-657.

[6]李建清,江丽军,莫元龙.等离子体-腔混合模耦合腔行波管非线性注-波互作用分析[J].强激光与粒子束,2004,16(11):1 429-1 433.

[7]张昭洪,吴鸿适.耦合腔慢波结构特性研究[J].电子学报,1986,14(1):7-8.

[8]刘盛刚.微波电子学导论[M].北京:国防工业出版社,1985.

[9]John H Mathews,Kurtis D Fink.数值方法[M].周璐,译.北京:电子工业出版社,2005:399-413.