| 标题 | 非协调元在温度应力场数值计算中的尝试 |
| 范文 | 何勇 张南南 杨锐 郭晓娜 强晟 摘要：在计算结构的温度应力时，采用普通的8节点协调元会产生较大的误差，其中一个重要原因是应变的精度与温度的精度不匹配。为此，引入非协调元理论，在位移插值函数中补充了不协调位移项，虽然在相邻单元边界上位移函数出现了不连续，但是可以使得单元位移函数中的2次项或3次项趋于完全，从而提高其精度。针对温度应力场计算的特点，采用8节点非协调单元的位移插值函数，并用Fortran 语言编成程序。算例验证结果表明：仿真计算中采用协调元计算温度场，采用非协调元计算应力场，其计算结果更接近理论值，是一种较好的匹配方案。 关键词：温度应力；协调元；非协调元；内部自由度 中图分类号：TV3文献标志码：A文章编号：1672-1683（2015）-002-0012-03 近年来，我国很多大型混凝土结构逐步开工建设，因此与其相对应的计算机仿真规模、精度、速度有了更高的要求。而目前常用的温度场应力计算仍存在一些不足，其中，用协调元计算温度场和应力场时，会产生温度应变和弹性应变的精度不匹配，从而造成较大误差。而非协调元在位移模式中包含完全的二次项，应变沿坐标轴可以线性变化的，因此，可考虑用非协调元计算应力场。 自从1973年Wilson[1]成功推导出Q6单元，关于内参型位移非协调元模型的研究取得了丰富的成果。其中，Strang and Fix[2]，Feng[3]对非协调单元稳定性问题进行了系统的数学分析。考虑收敛性问题，通常采用分片检验法。为了强使单元通过分片试验，Taylor[4]等人提出的单元Qm6。在国内，1982年，陈万吉[5]构造了一种新型的用于分析弹性空间问题的八节点六面体单元；1987年，长春[6]等人推导了非协调函数生成的一般公式；1989年， 鹿晓阳[7]等人提出Qmm6单元；九十年代以来，焦兆平[8]等人发展了许多新型内参型单元，探讨了内参型非协调元的合理位移场[9]：2000年，胡胜荣[10]等修正了分片检验的方法，形成了一种半解析半有限元分析方法；2000 年底，张春生[11]提出了内参型附加非协调位移基本项的概念，并且推导出了二维和三维情况[12]下的具体列式；2002年，张春生[13]等从基本力学出发，提出并证明了关于非协调元的两个定理；2010年，任国彪[14]等人构造了一种新的非协调长方体元，并研究了在三维线弹性问题中的稳定性和收敛性。 1非协调元理论 下面以三维8节点六面体单元为例，介绍非协调元的基本理论。 为了改善三维C0型8节点等参单元的性质，提高其精度，在单元的位移差值函数中附加内部无节点的位移项：α1（1-ξ2）、α2（1-η2）、α3（1-ζ2）。单元位移插值函数表示如下： u=∑9[]i=1Niui+α1（1-ξ2）+α2（1-η2）+α3（1-ζ2） v=∑9[]i=1Nivi+α4（1-ξ2）+α5（1-η2）+α6（1-ζ2） w=∑9[]i=1Niwi+α7（1-ξ2）+α8（1-η2）+α9（1-ζ2）（1） 其中：αi（i=1，…，9）为内部自由度，且 Ni=1[]9（1+ξξi）（1+ηηi）（1+ζζi）（i=1，2，…，8）（2） 将式（1）表示成矩阵形式： U=Nae+Nαe（3） 其中：U=（u v w）T、ae=（u1，v1，w1，…，w9）T、αe=（α1，α2，…，α9）T、N=[IN1，IN2，…，IN9] N=1-ξ21-η21-ζ2000000 0001-ξ21-η21-ζ2000 0000001-ξ21-η21-ζ2 （4） 代入几何关系可得： ε=Bae+Bαe（5） 由势能泛函，按照通常有限元的步骤取泛函变分为0，可得： KeuuKeuα KeαuKeαααe ae=Peu Peα（6） 由上式可解出 αe=Keαu-1Peα-Keαuαe（7） 利用上式消去内部自由度αe，得到凝聚后的单元求解方程为： Keae=Pe（8） 经过凝聚后，单元的自由度仍然是原 8 节点六面体等参单元的自由度，求解出单元节点位移后，由（7）式可以求得单元内部自由度，从而可以推求单元的应力。 2程序实现 在计算温度场时，仍然使用协调元，计算位移场和应力场时，采用非协调元。与普通协调元程序相比，非协调元程序不同之处主要包括：在位移模式中增加了非协调项；增加了单元内部自由度向节点的凝聚；根据节点位移反推内部自由度；求解单元应变时增加了单元内部自由度产生的应变。在一个荷载步中非协调元程序的计算流程图见图1。 图1一个荷载步中非协调元程序的计算流程 3算例 本算例考虑计算温度应力时，协调元和非协调元之间的差异，其中协调元采用C0型8节点等参单元，非协调元是相应的8节点非协调元。 计算模型如图2所示，长×宽×高为 5 m×5 m×1.5 m，单元总数为 2 500，节点总数为 3 380。该混凝土的弹模为 100 GPa，泊松比为 0.3，不考虑自重。导热系数为 1 000（kJ/（m·d·℃） ），导温系数为 1.0（m2/d），线膨胀系数为 1×105。模型初温为20 ℃，模型上表面为恒温0 ℃，下表面为恒温 20 ℃，四个侧面为绝热面。该计算模型无任何约束。求解该模型内部温度达到稳定状态时的应力。 图2计算模型 在温度达到稳定时，此模型内部为从0 ℃到 20 ℃线性变化的温度场，初始温度为20℃，因此温度增量也是线性变化的。根据线性的温度场在自由变形下不产生应力[15]，本算例的理论解是模型内无应力。现设置下面两个工况。 工况1：采用协调元计算温度场和应力场。 工况2：采用协调元计算温度场，采用非协调元计算应力场。 工况1的应力计算结果如图3-图6所示，工况2的应力计算结果如图7-图10所示，如无特殊说明，图中应力单位为 MPa。 图3工况1的σx分布云图 图4工况1的σz分布云图 图5工况1的τxy分布云图 图6工况1的τyz分布云图 图7工况2的σx分布云图 图8工况2的σz分布云图 图9工况2的τxy分布云图 图10工况2的τyz分布云图 由应力云图可见，工况 1的计算结果有较大误差，如图3所示，σx值达到0.6 MPa，而图4中的σz值高达1.5 MPa，但是理论解是无应力，故误差甚大。而工况 2 则与理论解非常接近，具有很高的精度，如图7、图8所示，其应力值均在0 MPa附近。工况2中切应力值也明显优于工况一中的计算结果。在工况 1的计算过程中，温度应变和温度具有同一阶精度，而弹性应变却比位移低一阶精度，这造成温度应变和弹性应变的精度不匹配，从而产生较大误差；在工况 2 中，由于位移采用非协调模式，弹性应变沿主应变方向也可以发生线性变化，从而可以更好地与温度应变相匹配，提高了计算精度。 4结语 本文将非协调元理论引入温度应力场的计算程序，并尝试在结构仿真计算中采用非协调元计算应力场，算例表明其精度明显高于协调元。在大体积混凝土温度场应力场仿真计算中，采用协调元计算温度场，采用非协调元计算应力场，可以使得温度应变的精度与温度的精度相匹配，从而达到较好的应力场计算精度，是一种较好的匹配方案。故可以在实际工程温控仿真计算中尝试推广，根据计算结果提出更合理的温控措施。 [HJ1.9mm]参考文献： [1]Wilson，E.L.Taylor，R.L.Doherty，W.P.and Ghabouss.Incompatible Displacement Methods[J]，Numerical and Computer Methods in Structural Mechanic，1973，43. [2]G.Strang and G.J.Fix.An analysis of the finite element method[M].PrenticeHall，Englewood Cliffs，N.J.，1973. [3] K.Feng.On the theory of discontinuous finite elements.[J].Math.Numer.Sinica，1979，1：378385. [4]Taylor，R.L.，Beresford，P.L.，Wilson，E.L.，A Nonconforming Element for Stress Analysis[J]，IJNME，1976，10（6）. [5]陈万吉.一个高精度八结点六面体单元[J].力学学报，1982（3）：262271. [6]ChangChun Wu，MaoGuang Huang and Theodore H.H.Pian，[JP2]Consistency condition and convergence criteria of incompatible elements： general formulation of incompatible functions and its application[J].Computers & Structures，1987，27（5）：639644. [JP] [7]鹿晓阳，刘玉文，许焕然，等人.Wilson非协调元的研究和改进[J].力学报，1989，21（3）：379384. [8]焦兆平，吴长春，黄茂光.内参型非协调元位移试解完备性的研究[J].中国科学技术大学学报，1992，22（3）：308316.[ZK）] |
| 随便看 |
|
科学优质学术资源、百科知识分享平台，免费提供知识科普、生活经验分享、中外学术论文、各类范文、学术文献、教学资料、学术期刊、会议、报纸、杂志、工具书等各类资源检索、在线阅读和软件app下载服务。