基于MCFT理论的割线刚度法介绍.doc

基于MCFT理论的割线刚度法介绍F. J. Vecchio于1990年在ACI发表的论文中采用了基于MCFT理论的割线刚度法来进行钢筋混凝土的有限元分析。通过比较割线刚度法与基于切线刚度的Newton-Raphson法，可以得出以下几点结论：1. 在结构的应力-应变曲线进入到下降段后，对于切线刚度法，刚度会出现负值，从而给计算处理带来难度。而割线刚度法则没有这种负刚度问题，但是割线刚度法会在卸载曲线穿越坐标轴时出现数值问题，这是由于此时割线刚度极小造成的。2. 以往采用切线刚度法的程序，在更新应变时，对于一个单元采用了4个高斯积分点，每次迭代得到了一个新的结构位移，就需要在4个积分点上分别求出它们的应变，再以此计算出4个新的割线刚度，用以判断收敛与否。而F. J. Vecchio采用的割线刚度法只涉及了1个积分点，即中心积分点。根据附表1可知，单元应变和位移的关系为：其中：令，则可以得到割线刚度法中心积分点的应变计算式为： 每次循环迭代出的新位移，仅需要换算成中心积分点上的应变，通过一次计算便可以得出割线刚度用于收敛判断。相比于切线刚度法，割线刚度法简化了计算过程。3. 切线刚度法的过程，需要通过外荷载节点力和结构恢复力求解出一个残余力，再以此计算出增量位移来判断是否满足收敛条件，若不收敛，便更新结构位移并重新计算结构的切线刚度矩阵。对于割线刚度法，并不需要求解残余力，而是直接利用外荷载节点力和结构割线刚度来求解出一个全量位移，之后通过计算出的应变来重新计算结构的割线刚度矩阵。二者相比起来，割线刚度法不需要启用一个新的数组来储存残余力值，降低了对程序的存储空间要求并提高了程序的运算速度。从上面三点可以看出，切线刚度法和割线刚度法各有特点和不足，在某些方面，割线刚度法还体现出了一定的优势：1. 割线刚度法相比于切线刚度法的收敛准则更为简单。2. 割线刚度法相比于切线刚度法的有限元模型更为简单。3. 割线刚度法的稳定性较好。所以对割线刚度法进行研究佐证是十分有必要的。F. J. Vecchio在其论文中所述的割线刚度法的具体过程如下：第1步：输入结构和材料的属性数据第2步：计算外荷载节点力第3步：计算材料的刚度矩阵：1) 确定初始应变：，； 2) 计算主方向应变，以及主拉应变角；3) 根据混凝土的本构关系，确定混凝土主拉和主压应力，4) 根据钢筋的本构关系，确定钢筋的主拉和主压应力，；5) 计算混凝土和钢筋的割线刚度元素值： (3-1) 第4步：计算混凝土和钢筋的材料割线刚度矩阵1) 计算混凝土主方向上的材料割线刚度矩阵： (3-2)2) 计算转换矩阵3) 计算整体坐标系（X-Y坐标）下的混凝土割线刚度矩阵 (3-3)4) 计算整体坐标系(X-Y坐标系)下的钢筋割线刚度矩阵： (3-4)第5步：计算组合材料刚度： (3-5)第6步：计算单元刚度矩阵： (3-6)第7步：计算单元预应变产生的位移，第8步：计算单元预应变产生的节点力，第9步：计算单元节点力：第10步：计算结构总刚矩阵： (3-7) 图 3.5 单元尺寸第11步：计算结构总刚矩阵的逆矩阵：第12步：计算节点位移 Fig3.5 Element dimension (3-8) 第13步：计算单元应变（取单元形心处的应变为平均应变）： (3-9) (3-10) (3-11)注：，为相应节点的水平位移和竖向位移。第14步：检查单元应力：第15步：重新计算单元割线刚度值：1) 根据求解出的位移，重新确定主方向上的应变值，；3) 根据混凝土本构模型，计算新的混凝土应力，；4) 根据钢筋本构模型，计算新的钢筋应力，；5) 重新计算钢筋和混凝土的割线刚度元素：，,。第16步 收敛判断详细的算法流程图如下所示：2.1 修正斜压场理论（MCFT）介绍修正斜压场理论（The Modified Compression-Field Theory）是一种用于预测钢筋混凝土膜构件在剪力和正应力作用下的变形性能的分析模型，它可以确定混凝土和钢筋的平均应力应变和局部的应力应变，以及裂缝的宽度和方向，从而确定构件的破坏模式。MCFT是由研究钢筋混凝土受扭和受剪的斜压场理论（Compression-Field Theory）发展而来。斜压场理论忽略了开裂混凝土的拉力，而修正斜压场理论则把裂缝间的拉应力考虑了进来。MCFT把开裂混凝土模拟为一种正交各向异性的材料，是一种完全的转角裂缝模型（Smeared Rotating Crack Model）。与固体的不连续性相悖，开裂混凝土被视作一种裂缝均匀分布的固体连续体。裂缝可以自由转动，并保持和混凝土的主压应力场方向同向。修正斜压场理论包含了三种关系：用平均应变表示的混凝土和钢筋的变形协调关系；用平均应力表示的混凝土和钢筋的应力平衡关系；用平均应力和平均应变表示的混凝土和钢筋的本构关系。开裂混凝土的本构关系是根据多伦多大学所做的钢筋混凝土剪切板实验的结果得到的，所以MCFT采用的是实验观测到的真实的本构模型。同时，由于裂缝均匀分布以及上述三种关系用平均应力和平均应变的概念表述，所以MCFT的一个重点是要考虑裂缝处的局部应力应变条件。图2.1 开裂钢筋混凝土的修正压场理论23Fig. 2.1 Modified Compression Field Theory for cracked reinforced concrete23修正斜压场理论作了以下一些假设：1. 钢筋弥散分布于单元之中；2. 裂缝均匀分布并能自由旋转；3. 单元中，剪应力和正应力均匀分布；4. 应力状态和应变状态一一对应，不考虑加载历史的影响；5. 跨越多条裂缝间距，应力应变可视作是平均的；6. 主应力角和主应变角的方向一致；7. 钢筋和混凝土之间无粘结滑移；8. 混凝土和钢筋有独立的本构关系；9. 忽略钢筋中的剪应力。2.1.1 变形协调条件在混凝土和钢筋构件中，与平均应变相关的变形协调关系如图 所示： 图2.2 混凝土平均应变28Fig.2.2 Average concrete strains28因为假设钢筋和混凝土之间没有粘结滑移，所以对于非预应力钢筋而言，钢筋和混凝土中的平均应变是相等的。变形协调条件可以用一下公式表达：若剪应变值为，则根据莫尔应变圆确定的混凝土平均主拉应变和平均主压应变为：通过莫尔应变圆也可以确定平均主拉应变和应力与轴的夹角：2.1.2 应力平衡关系从构件中取出一个自由体作为研究对象，如图 所示：和方向的力平衡条件要求混凝土平均应力和以及钢筋平均应力和的合力要与施加的正应力和的合力相平衡。弯矩平衡条件则要求混凝土中的平均剪应力要能完全抵抗所施加的剪应力（假设钢筋没有销栓作用）。平均应力的平衡关系可以表述如下：和分别是和方向的配筋率。因为开裂混凝土是主应力方向上的正交各向异性材料，莫尔应力圆可以用于建立混凝土平均应力和与混凝土平均主拉应力之间的关系式：2.1.3 本构关系本构模型需要把变形协调条件中的应变和平衡条件中的应力联系起来。1986年，多伦多大学的Vecchio和Collins做了30个板实验，这些板的尺寸都是，并承受平面内的应力。通过对实验结果的分析，他们得出了开裂混凝土受拉和受压的本构模型。对于受压混凝土，本构关系是主压应力和主压应变的函数。上述板的实验结果表明，在主拉应变同时存在的情况下，抗压强度和刚度会随着的增大而减小。这一现象被称为受压软化，具体表现在混凝土单轴受压时的应力-应变反应的弱化。建议的考虑受压软化的公式如下：其中，是单轴受压混凝土的Hognestad抛物线函数，常用于一般强度的混凝土。是对应于峰值抗压应力的混凝土圆柱体应变（取负值），由混凝土圆柱体单轴受压实验确定。反映的是主拉应变的软化效应。对于受拉混凝土，本构关系是主拉应力和主拉应变的函数。首先，必须要确定开裂强度以及相应的开裂应变。如果信息不足，可以作如下估计： 其中是混凝土的初始刚度，计算式为： 在开裂以前，混凝土在受拉状态下是线弹性的： 开裂以后，由于混凝土和钢筋的粘结作用，拉应力将继续存在于钢筋混凝土裂缝间的混凝土中。为了模拟这种受拉硬化现象，混凝土的拉应力要随着混凝土主拉应变的增大从抗拉强度逐步减小。MCFT建议的关系公式如下：对于处在受拉和受压状态的钢筋而言，MCFT采用了平均应力和平均应变的两折线模型：一个线弹性的上升段和一个屈服平台。公式如下：其中，是钢筋的弹性模量，和分别是钢筋在和方向的屈服应力。2.1.4 考虑局部平衡条件有了相容的应变条件，根据本构关系可以确定混凝土和钢筋中的平均应力，以及它们所平衡的施加剪应力和正应力。然而，如果忽略了构件性能由裂缝处的局部钢筋屈服或沿裂缝的剪切滑移破坏所控制的可能性，那么分析结果将会是不保守的。为了把这些可能性考虑进来，MCFT限制了裂缝处的局部应力以及混凝土的平均拉应力水平。钢筋混凝土的应力场从裂缝间的平均概念到裂缝处的局部概念是不断变化的。图 a所示的是裂缝间垂直于主拉应力方向的截面上的平均应力情况，而图 b所示的是裂缝处截面上的局部应力情况。在裂缝处的截面上，混凝土平均拉应力减小到0。为了跨过裂缝传递平均拉应力，裂缝处的钢筋应力和应变就必须局部提高。在裂缝表面的法线方向根据力平衡条件，用平均应力和局部应力建立的数值等价关系式可以表述为：其中，和是裂缝处的局部钢筋应力，和是钢筋和裂缝截面法线方向的夹角。通过上面的公式可以明显看出，混凝土的平均拉应力受到了裂缝处钢筋屈服的限制。如果用钢筋的屈服强度取代钢筋局部应力，那么公式两个括号中的项就表示了钢筋的富余强度，也就限制了后开裂混凝土的受拉应力不超过：如图 a所示，主平面上没有剪应力。但如果钢筋与裂缝斜交，那么就会在裂缝表面产生局部剪应力。在裂缝切线方向上根据力的平衡关系，可以得到如下公式：不考虑上述方程，局部剪应力可以一直增大，直到发生剪切滑移破坏。剪应力的大小会受到骨料咬合作用的限制，这种骨料咬合作用随着裂缝宽度的增大和最大骨料尺寸的减小而减小。基于Walraven对骨料咬合作用的分析，MCFT对裂缝处剪应力大小的限制为：平均裂缝宽度是混凝土主拉应变和平均裂缝间距的函数：和方向的平均裂缝间距和可以根据粘结属性和钢筋分布情况求出。如果超过了允许的混凝土最大拉应力或者裂缝处的局部剪应力，就需要修正应变状态以求得到较小的混凝土平均拉应力值。2.2 混凝土滞回本构模型图 a所示为混凝土受压的滞回模型，其骨架曲线表示的是结构在单调加载时的反应。根据MCFT可知，混凝土单调本构曲线是以Hognestad抛物线为基础，通过考虑受压软化效应，对其进行修正得到的。设混凝土塑性应变为，当再加载时，混凝土的压应力可以通过下式计算：其中，是之前加载得到的最大压应变，是与相应的应力，是通过骨架曲线计算得到的与相应的应力。如果结构的反应落在骨架曲线上，如，那么就将和更新为和。在每一个荷载步，即时的塑性应变用下式计算：其中，是与骨架曲线上峰值应力相对应的应变，如果即时塑性应变大于塑性补偿应变，那么就将后者更新。任意时刻的卸载产生的应力用如下关系计算：其中卸载模量定义为：这一模型隐含的假设是，当混凝土由拉变为压时，压应力保持为0，直到裂缝完全闭合，即直到。实际上，实验证明这个假设是不正确的。由于之后的应变总是会大于0，且受到裂缝剪切滑移的影响，此时压应力不会等于0。混凝土受拉的滞回本构如图 b所示。骨架曲线由两部分组成：开裂前的曲线和开裂后的受拉硬化曲线。设塑性应变为，再加载时的混凝土拉应力用下式计算：其中，是之前加载得到的最大拉应变，是与相对应的应力，是通过骨架曲线计算得到的与相对应的应力。对于第一个受拉加卸载，塑性补偿应变为常数，其值等于从受压区穿轴确定的塑性应变大小。接下来用于骨架曲线计算的是净塑性应变。在之后的循环中，塑性补偿应变用下式计算：由于目前缺乏合适的模型，从上式可以看出，没有正值的补偿应变。这就使得在加载和卸载的过程中，受拉滞回曲线都会经过原点。当混凝土反应落在骨架曲线上时，最大拉应变和相应的应力都分别更新到和。卸载时的应力可以通过下式计算：其中，卸载模量定义为：2.3 钢筋的滞回本构模型本文的钢筋的应力-应变关系是采用的经Filippou等人改进后考虑了各向同性应变强化的Menegotto-Pinto非线性模型。Menegotto-Pinto模型的表达式为： （2.64）其中， （2.65） （2.66）方程(2.64)表达的是斜率分别为和的两条渐近线的过渡段，如图 所示的线（a）和线（b）。点B是两条渐近线的交点，其应力和应变分别为和。 A点是荷载的反向点，其应力和应变分别为和。b为和的比值，叫做应变强化率。是影响过渡曲线形状的一个参数，它反映了钢筋的包兴格（Bauschinger）效应。图2.13 Menegotto-Pinto 钢筋模型Fig 2.13 Menegotto-Pinto steel model的值取决于当前渐进线交点B的应变和上一次荷载的反向点A的最大或最小应变间的应变差。Menegotto-Pinto模型中的表达式为： （2.67）其中， 是第一次加载时参数的初始值，和跟一样，其值由钢筋试验所确定，将在每次应变反向后更新其值。图2.14 Menegotto-Pinto 钢筋模型中曲率参数的定义34Fig 2.14 Definition of curvature parameter in Menegotto-Pinto steel model34上述Menegotto-Pinto模型主要缺点是，它无法考虑等向应变硬化的影响。针对这个问题，Fillipou (1983)建议把线性屈服渐近线上的应力变换看作是最大塑性应变的函数： （2.68）其中，和分别是钢筋屈服时的应变和应力，是当应变反向时绝对值最大的应变， 与的值由试验确定。考虑了等向应变硬化影响修正后的钢筋滞回曲线如图 所示。图2.15 钢筋滞回本构模型34Fig 2.15 Hysteresis Model for reinforcrment342.4 基于切线刚度的Newton-Raphson迭代法邓兴龙和蒋晓华等编制的非线性有限元程序是采用的基于切线刚度的Newton-Raphson迭代法理论得以实现的。下面对其所涉及的高斯积分以及非线性方程组解法中的Newton-Raphson迭代算法、弧长法和收敛准则做一些简单的介绍。2.4.1 高斯积分基于切线刚度的Newton-Raphson迭代法采用的是积分点的高斯积分，其单元刚度矩阵和单元等效结点恢复力的高斯积分公式介绍如下。至于具体的高斯求积公式的推导过程，可以参见文献。用高斯求积公式表述的单元刚度矩阵计算式为：用高斯求积公式表述的单元等效结点恢复力的计算式为：其中，为应变矩阵，为雅克比行列式。和的值与高斯积分点相对应。针对不同的单元，一般通过经验和试算来确定最佳的积分点数。对于四边形单元而言，的高斯积分点就可以较好地满足精度要求。表 中所示的是各个高斯积分点的坐标和加权系数，图 所示的是高斯积分点的分布情况。表2.3 高斯积分中的和 Table 2.3 and in Gauss Integration点数n积分点坐标加权系数点数n积分点坐标加权系数1024-0.8611363116-0.33998104360.34785484510.65214515492-0.57735026921130.774596669200.5555555560.8888888895-0.9061798459-0.538469310100.23692688510.47862867050.5688888889图2.9 高斯点的位置Fig 2.9 The location of Gauss Points2.4.2 Newton-Raphson迭代法基于切线刚度的Newton-Rapson迭代法如图3.6。0000000000P图3.6 基于切线刚度的Newton-Raphson算法Fig3.6 Newton-Raphson algorithm based on tangent stiffnessNewton-Raphson是一种变刚度的迭代法。首先，取单元的初始刚度为，根据所施加的外荷载等效结点力可以求出位移的第一次近似值为： (3.1)由得到的初始位移可以计算单元的应变值，再结合本构模型求到单元应力，进而得到相应的结点恢复力。然后，用外荷载等效结点力和结构恢复力求出残余力，并连同相应于初始位移的即时切线刚度，求出位移增量，即 (3.2) (3.3)得到位移的第二次近似值为： (3.4)以此类推，第k步的计算过程为： (3.5) (3.6) (3.7)当与的误差达到所规定的限值，或者当足够小的时候，迭代结束。 2.4.3 弧长控制法弧长控制法的基本原理是，通过在结构分析的一般方程中添入一个新的变量，来实现对整个计算过程的自动控制。加入变量后，原方程组将变为： （2.99）上述方程组一共有个未知量：、，但是方程却只有个。所以，需要补充如下条件来完成方程组的求解： （2.100）称为“弧长”，是附加的控制变量。这时，方程组变为。弧长法的迭代过程如图 所示。图2.14 弧长法迭代过程示意图Fig 2.14 The Arc-length method2.4.4 收敛准则基于切线刚度的Newton-Raphson迭代法的收敛准则有三种：力准则、位移准则以及能量准则。选取合适的收敛准则，有利于提高迭代法的计算效率，是增量求解策略中一个重要的组成部分。当迭代结束时，结构应当处于一个稳定的力平衡状态，力准则就是要求结构的恢复力和外力之差等于零。但是在数值计算中，不可能实现不平衡力完全为零，这就需要设定一个允许值。力准则的表达式为：其中，就是给残余力设置的一个允许值。在进行位移分析时，计算位移应该与真实值接近。位移准则的表达式为：其中，是矢量的欧几里德范数，是第步迭代的位移增量，为设置的允许值。能量准则是用来度量力和位移同它们的平衡值之间的接近程度的，其表达式为： （2.115)其中，不等式的左边表示残余力在位移增量上所做的功，不等式的右边是残余力所做的功的初始值，是为能量设置的允许值。邓兴龙和蒋晓华等编制的FEAPpv程序采用的是能量收敛准则。2.5 基于MCFT理论的割线刚度迭代法F. J. Vecchio于1990年在ACI发表的论文中采用了基于MCFT理论的割线刚度迭代法来进行钢筋混凝土的有限元分析。通过比较割线刚度迭代法与基于切线刚度的Newton-Raphson法，可以得出以下几点结论：1. 在结构的应力-应变曲线进入到下降段后，对于切线刚度法，刚度会出现负值，从而给计算处理带来难度。而割线刚度迭代法则没有这种负刚度问题，但是割线刚度迭代法会在卸载曲线穿越坐标轴时出现数值问题，这是由于此时割线刚度极小造成的。2. 以往采用切线刚度法的程序，在更新应变时，对于一个单元采用了4个高斯积分点，每次迭代得到了一个新的结构位移，就需要在4个积分点上分别求出它们的应变，再以此计算出4个新的割线刚度，用以判断收敛与否。而F. J. Vecchio采用的割线刚度迭代法只涉及了1个积分点，即中心积分点。根据附录1可知，单元应变和位移的关系为：其中：令，则可以得到割线刚度迭代法中心积分点的应变计算式为： 每次循环迭代出的新位移，仅需要换算成中心积分点上的应变，通过一次计算便可以得出割线刚度用于收敛判断。相比于切线刚度法，割线刚度迭代法简化了计算过程。3. 切线刚度法的过程，需要通过外荷载节点力和结构恢复力求解出一个残余力，再以此计算出增量位移来判断是否满足收敛条件，若不收敛，便更新结构位移并重新计算结构的切线刚度矩阵。对于割线刚度迭代法，并不需要求解残余力，而是直接利用外荷载节点力和结构割线刚度来求解出一个全量位移，之后通过计算出的应变来重新计算结构的割线刚度矩阵。二者相比起来，割线刚度迭代法不需要启用一个新的数组来储存残余力值，降低了对程序的存储空间要求并提高了程序的运算速度。从上面三点可以看出，切线刚度法和割线刚度迭代法各有特点和不足，在某些方面，割线刚度迭代法还体现出了一定的优势：4. 割线刚度迭代法相比于切线刚度法的收敛准则更为简单。5. 割线刚度迭代法相比于切线刚度法的有限元模型更为简单。6. 割线刚度迭代法的稳定性较好。所以对割线刚度迭代法进行研究佐证是十分有必要的。2.5.1 算法流程以图 所示的单元为例，Vecchio在其论文中所述的割线刚度迭代法的具体计算过程如下：图 3.5 单元尺寸Fig3.5 Element dimension第1步：输入结构和材料的属性数据第2步：计算外荷载节点力第3步：计算材料的刚度矩阵：1) 确定初始应变：，； 2) 计算主方向应变，以及主拉应变角；3) 根据混凝土的本构关系，确定混凝土主拉和主压应力，4) 根据钢筋的本构关系，确定钢筋的主拉和主压应力，；5) 计算混凝土和钢筋的割线刚度元素值： 第4步：计算混凝土和钢筋的材料割线刚度矩阵1) 计算混凝土主方向上的材料割线刚度矩阵： (3-2)2) 计算转换矩阵3) 计算整体坐标系（X-Y坐标）下的混凝土割线刚度矩阵 (3-3)4) 计算整体坐标系(X-Y坐标系)下的钢筋割线刚度矩阵： (3-4)第5步：计算组合材料刚度： (3-5)第6步：计算单元刚度矩阵： (3-6)第7步：计算单元预应变产生的位移第8步：计算单元预应变产生的节点力第9步：计算单元节点力：第10步：计算结构总刚矩阵： (3-7) 第11步：计算结构总刚矩阵的逆矩阵：第12步：计算节点位移 (3-8) 第13步：计算单元应变（取单元形心处的应变为平均应变）： (3-9) (3-10) (3-11)注：，为相应节点的水平位移和竖向位移。第14步：检查单元应力：第15步：重新计算单元割线刚度值：1) 根据求解出的位移，重新确定主方向上的应变值，；3) 根据混凝土本构模型，计算新的混凝土应力，；4) 根据钢筋本构模型，计算新的钢筋应力，；5) 重新计算钢筋和混凝土的割线刚度元素：，,。第16步：收敛判断详细的算法流程图如图 所示。割线刚度迭代法的迭代过程与基于切线刚度的Newton-Raphson迭代法有一定的可比性，其具体的迭代过程图线如图 所示。每次迭代通过外荷载与即时的割线刚度求解出全量位移，直到与真实位移的误差小于允许值，迭代计算结束。2.5.2 收敛准则不同于基于切线刚度的Newton-Raphson迭代法，割线刚度迭代法并没有采用前面所述及的三种收敛准则。由于这三种收敛准则都是贯彻了增量策略的思想，对于应用全量理论的割线刚度迭代法而言，显然不太合适。因此，割线刚度迭代法采用了一种全新的收敛准则：以单元的割线刚度作为判断收敛的条件。如果某前后两个迭代步得到的割线刚度误差小于预定的限制，那么迭代终止。对于多个单元的情况而言，割线刚度收敛准则可以表现为不同的形式。其中一种是以各个单元的割线刚度误差的最大值作为判断标准，即其中，是向量的最大范数，分别是第个单元的混凝土两个主方向上的迭代步前后两次割线刚度误差，是预设的误差限。另外一种形式是以各个单元的割线刚度误差之和作为判断准则，数学表达式为：其中，是向量的1-范数，是预设的误差限。 本文程序选用的是第二种形式的收敛准则。3.1 FEAPpv简介在过去的几十年里，有限元法从一种线性结构分析方法演变成了一种求解非线性微分方程的通用方法，FEAPpv便是应运而生的一个有限元分析程序（Finite Element Analysis Program）。在研究和应用领域中，需要不断对其进行改进来解决新的问题和满足新的分析要求。FEAPpv程序主要开发于UNIX平台和Window平台，其包括一系列模块用于实现：1.有限单元模型数据的输入；2. 用于处理大范围应用实例的求解方法的建立；3. 结果的图形和数值输出。一个问题的求解，是通过使用命令语句的形式来达到的，其中的求解算法则由用户定义并编写。因此，FEAPpv的这一特点使得每一个用户可以根据自己的具体需要来定义一种特定的求解方式，以此达到分析研究目的。FEAPpv程序自身包含了足够多的命令来应用于结构力学、流体力学、热传导以及许多其他需要用微分方程来模拟问题的领域，静态和动态问题都可以用该程序求解。根据不同的应用需求，用户也可以自行添加新的模型和命令语句，这些功能可以帮助用户对一些特殊结构进行单元的划分，更或者使用户能从其他程序中导入单元划分信息。FEAPpv程序包含一个有限的单元库。这些单元可以用于一维、二维和三维线性或非线性结构问题、固体力学问题以及线性热传导问题的建模，并调用一个材料模型库。所提供的材料模型适用于弹性、粘滞弹性、塑性以及热传导本构方程，用户也可以根据自己的需要添加另外的模型。在结构问题中，程序能够通过单元形成质量和几何刚度矩阵，并计算单元上的相关值，如应力和应变，并能够把相关值投影到结点上，输出图形结果。FEAPpv程序可分为三个基本模块：1. 数据输入模块和先处理器2. 求解模块3. 输出模块以下是简化的FEAPpv程序流程图：上图中的程序输入模块将读入输入文件中的材料（MATE）、几何（COOR， ELEM）、边界条件（BOUN）、荷载（FORC）以及位移（DISP）信息，用以建立所有的单元数组。之后，材料参数、结点荷载和结点位移、方程编号、温度、单元连接以及坐标等分别存储在程序中的不同数组里。表 列出了一些主要数组的名称、指针编号和它们各个维度的大小。表3.1 网格划分数组名称、指针编号和数组的大小32Table 3.1 Mesh Array Names, Numbers and Sizes32NAMENum.dim 1dim 2dim 3DescriptionD25nddnummat-Material parametersF27ndfnumnp2Force and DisplacementID31ndfnumnp2Equation nos.IE32nienummat-Element control, dofs, etc.IX33nen1numel-Element connectionsT38numnp-TemperatureU40ndfnumnp3Solution arrayVEL42ndfnumnpntSolution rate arrayX43ndmnumnp-Coordinates本文主要涉及和研究FEAPpv的控制程序，至于详细的前处理模块介绍，可以参见文献。下图为控制程序的简化流程图。可以看出，整个控制程序由网格划分（PMESH）和宏命令（PMACR）求解两大部分组成。FEAPpv的宏（MACR）命令流包含与6个子程序，分别是PMACR、PMACIO、PMACR1、PMACR2、PMACR3和PMACR4，其中PMACR4是一个用户自定义宏命令子程序，不使用时为空，且不会对整个程序运行有影响。这几个宏命令子程序的具体任务和作用是：PMACR位于流程上端，负责调用其他几个子程序；PMACIO负责读入命令流，并检查宏命令的循环和嵌套；PMACR1主要负责程序的计算部分，如计算刚度和残余力；PMACR2主要负责控制计算过程的循环、荷载步长以及收敛误差的大小等；PMACR3则负责完成求解和绘图等其他的功能。图 所示为整个宏命令的结构图：3.2 钢筋混凝土基于切线刚度法的非线性分析流程图FEAPpv的源程序采用的是典型的基于切线刚度的有限元计算方法。切线刚度法的特点是：1. 需要通过输入信息来求得单元的结点恢复力，并集成结构的结点恢复力，再与外荷载产生的结点力求残差，确定出残余力的大小；2. 根据残余力和形成的总刚计算结构的增量位移，并更新总位移；3. 以力、位移或者能量作为准则判断是否收敛。钢筋混凝土基于切线刚度法的二维非线性分析具体流程图如下所示：3.3 钢筋混凝土基于割线刚度法的非线性分析流程图本文所研究的钢筋混凝土基于割线刚度的二维非线性有限元分析流程图如下所示可以看出，割线刚度法和切线刚度法算法流程的差别在于：1. 割线刚度法不需要计算单元恢复力以及结构的残余力，从而不必分配大型的数组用以存取恢复力和残余力，节省了计算时间和储存空间，提高了效率；2. 割线刚度法直接利用外荷载在结点上产生的力和集成的结构总体刚度矩阵来求解结构的全量位移，而非切线刚度法中通过残余力求得的增量位移；3. 割线刚度法以单元割线刚度元素作为收敛条件，相比于切线刚度法的收敛判断准则更为简洁。3.4 单元子程序的数据传递和结构用户可以通过ELMTnn这个单元子程序来向FEAPpv的单元库添加自己所需要的单元模块：subroutine ELMTnn (d, ul, xl, ix, tl, s, r, ndf, ndm, nst, isw)其中nn的值从01到05。单元子程序中的各种信息，是通过哑元和公共变量进行传递的。表 列出了程序中的一些哑元以及它们所传递的数据信息，图 列出了一部分公共块的名称及其所含的参数变量。表3.2 FEAPpv单元子程序的哑元32Table 3.2 Arguments of FEAPpv Element Subprogram32ParameterDescriptiond(*)Element data parametersul(ndf,nen,j)Element nodal solution parameters nen is number of nodes on an element (max)xl(ndm,nen)Element nodal reference coordinatesix(nen)Element global node numberstl(nen)Element nodal temperature valuess(nst,nst)Element matrix (e.g., stiffness, mass)r(ndf,nen)Element vector (e.g., residual, mass)may also be used as r(nst)ndfNumber unknowns (max) per nodendmSpace dimension of meshnstSize of element arrays S and RN.B. Normally nst = ndf*neniswTask parameter to control computationinteger numnp, numel, nummat, nen, neq,iprcommon /cdata/ numnp, numel, nummat, nen, neq,iprreal*8 dminteger n, ma, mct, iel, nelcommon /eldata/ dm, n, ma, mct, iel, nelinteger nh1,nh2,nh3common /hdata/ nh1,nh2,nh3integer ior,iowcommon /iofile/ ior,iowinteger ndf,ndm,nen1,nst,nneqcommon /sdata/ ndf,ndm,nen1,nst,nneqreal*8 hrinteger mrcommon / / hr(1),mr(1000)图3.5 FEAPpv 单元公共块32Figure 3.5 FEAPpv Element Common Blocks32FEAPpv中单元子程序的简单结构图如图 所示。图 与切线刚度法对应，图 与割线刚度法对应。其中ISW是FEAPpv程序中起控制作用的任务参数，相应于不同取值的ISW的具体任务见表表3.3 FEAPpv单元子程序任务选项32Table 3.3 Task Options for FEAPpv Element Subprogram32 isw taskDescriptionAccess Command1Input d(*) parametersMesh:MATE,n2Check elementsSoln:CHECk3Compute tangent/residual Soln:TANGStore in S/rSoln:TANG4Output element variablesSoln:FORM,REACPlot:REAC6Compute residualSoln:FORM,REACPlot:REAC8Nodal projectionsSoln:STRE,NODEPlot:STRE,PSTR12History updateSoln:TIME14Initialize historyBATCh,INTEr通过对比两个结构图，也可以看出，割线刚度法省去了计算恢复力的步骤，更为简便。3.4.1 单元子程序中材料参数的存储在程序中，单元的各种材料参数存储的一维数组D(*)里。D(*)数组中各个元素的具体代表值见表表3.1 钢筋混凝土单元材料参数Table 3.1 Material Parameters of the RC ElementParameterDescriptionD(1)混凝土的初始弹性模量 D(2）混凝土的初始泊松比D(3)混凝土抗压峰值强度D(4)混凝土的受压峰值应变D(5)混凝土受拉开裂强度D(7)钢筋弹性模量D(8)钢筋的应变强化模量D(10)x向钢筋的屈服强度D(11)y向钢筋的屈服强度D(12)x向钢筋的配筋率D(13)y向钢筋的配筋率D(14）单元的厚度D(181)X向钢筋直径D(182)Y向钢筋直径D(183)混凝土骨料尺寸3.4.2 单元应力和材料矩阵的形成过程钢筋混凝土单元的应力包含两个部分：钢筋应力和混凝土应力。同样的，钢筋混凝土单元的材料矩阵也由钢筋的材料矩阵和混凝土的材料矩阵组合而成。图 所示的是单元应力和材料矩阵形成的过程图。对于切线刚度法而言，图示的流程是针对每个高斯积分点进行的，而在割线刚度法中，整个过程在单元的中心点得以实现。3.4.3 单元刚度矩阵的形成切线刚度法采用的是FEAPpv程序默认的高斯点数，即用22的高斯点进行积分，其具体的单元刚度矩阵形成过程见图对于割线刚度法，在每个单元上仅取其形心点处的应变作为整个单元的平均应变，一定程度上简化了单元刚度形成的计算过程。3.4.4 单元历史状态的存储 在进行低周反复加载的非线性分析时，对单元的历史状态进行存储甚为重要。比如结构在经过了加卸载之后，会产生不可恢复的塑性应变，而下一荷载步通过公式计算得到的塑性应变值很有可能小于上一荷载步计算所得值，若程序没有记忆功能，那么小值将会把大值覆盖，得出不合理的结果。所以，在程序中需要设置一些这种参数：它们负责把单元有用的历史变量存储起来，如上一荷载步的应变以及加载过程中最大的塑性应变等。混凝土单元在某个迭代收敛步处所设置的历史变量如下表所示：表3.2 混凝土单元材料参数h(*)Table 3.2 Material Parameters of the Concret ElementParameterDescriptionh(1)混凝土主应力角h(2)混凝土X方向塑性应变h(3)混凝土Y方向塑性应变h(4)混凝土塑性应变h(5)混凝土X方向最大压应变h(6)混凝土Y方向最大压应变h(7)混凝土最大压应变h(8)混凝土X方向最大拉应变h(9)混凝土Y方向最大拉应变h(10)混凝土最大拉应变h(11)混凝土X方向上一步应变h(12)混凝土Y方向上一步应变h(13)混凝土上一步应变钢筋单元在某个迭代收敛步处所设置的历史变量如下表所示：表3.3 钢筋单元材料参数h(*)Table 3.3 Material Parameters of the Steel ElementParameterDescriptionh(18)，h（19）X，Y向钢筋上一步路径h(20)，h(21)X,Y向钢筋上一步应变h(22)，h（23）X,Y向钢筋上一步应力h(24)，h（25）X,Y向钢筋上一步塑性应变h(26)，h（27）X,Y向钢筋渐近线交点应变h(28)，h（29）X,Y向钢筋渐近线交点应力h(30)，h（31）X,Y向钢筋应变反向点应变h(32)，h（33）X,Y向钢筋应变反向点应力h(34)，h