轴心受压构件的弹塑性屈曲荷载的计算方法

技术领域

本发明涉及一种结构设计中的数值计算方法，具体涉及一种轴心受压构件 的弹塑性屈曲荷载的计算方法。

背景技术

关于轴心受压构件的弹塑性屈曲问题，1889年Engesser.F.提出了切线模量 理论，建议用变化的变形模量E，代替欧拉公式中的弹性模量E，从而获得轴心 受压构件的弹塑性屈曲荷载。然而，当构件微弯时，凹面的压应力增加，而凸 面的应力减少，它们遵循着不同的应力-应变关系。因而，1891年Considere.A. 在论文中阐述了双模量的概念，在此基础上1895年Engesser.F.提出了双模量理 论，建议用与E_t和E都有关的折算模量E_r计算屈曲荷载。但是，试验资料表明， 实际的屈曲荷载介于两者之间而更接近于切线模量屈曲荷载。直到1946年 Shanley.F.R.提出构件在微弯状态下加载时凸面可能不卸载的概念，并用力学模 型证明了切线模量屈曲荷载是弹塑性屈曲荷载的下限，而双模量屈曲荷载是其 下限。轴心受压构件在微弯状态下可以继续加载的概念与前面已经阐明的大挠 度理论是一致的。新的切线模量理论可以广泛地用于解决稳定分岔失稳类型构 件或板的非弹性屈曲问题。

对于轴心受压构件，切线模量屈曲荷载的基本假定如下：

(1)构件是挺直的。

(2)构件两端铰接，荷载沿构件轴线作用。

(3)构件的弯曲变形很微小。

(4)弯曲前的平截面在弯曲后仍为平面。

(5)在弯曲时全截面没有出现反号应变。

作用于截面的压力和内力矩分别为

P＝∫_AσdA

M＝∫_AΔσzdA＝E_tIΦ＝-E_tIy″

构件的平衡方程为

EI_ty″+Py＝0

20世纪中，美国Lehigh大学Fritz工程结构实验室Huber.A.W.和Beedle， L.S.等人对构件中残余应力的分布和其对轴心受压构件屈曲荷载的影响进行了 系统研究，发现和初始几何缺陷对轴心受压构件的影响一样，残余应力也是影 响屈曲荷载的重要因素。残余应力将使构件提前进入弹塑性状态，降低构件的 屈曲载荷。不同残余应力分布对于轴心受压构件屈曲荷载的影响有很大差异， 其中以位于截面外侧且具有很高残余压应力峰值时对屈曲荷载的影响最为显 著。

一般双轴对称截面轴心受压构件，可能绕截面的两个对称轴发生弯曲屈曲， 但是有些抗扭刚度和抗翘曲刚度很弱的轴心受压构件，如双轴对称的十字型截 面轴心受压构件，除了有可能发生绕对称轴x或y的弯曲屈曲外，还有可能发 生绕截面纵轴z转动的扭转屈曲，此时纵轴本身不发生弯曲变形，只是截面绕 纵轴旋转一个角度。

另外，对于单轴对称截面轴心受压构件，如等边角钢轴心受压构件，除了 可能绕截面的非对称轴x发生弯曲屈曲外，还可能在绕截面的对称轴y弯曲的 同时又绕通过截面剪心S的纵轴扭转而发生弯扭屈曲。截面的位移有侧移u和 扭转角如果截面不具有对称轴，如不等边角钢截面(其截面的形心和剪心不 重合)，这种截面的轴心受压构件只可能发生弯扭屈曲，如图1(a)和1(b) 所示，此时截面不仅绕两个主轴x和y弯曲，而且还绕通过剪心的纵轴扭转。

理想的轴心受压构件的弯曲屈曲、扭转屈曲和弯扭屈曲都属于稳定分岔屈 曲失稳问题。除了截面的残余应力分布比较简单的可以用解析法得到轴心受压 构件的屈曲荷载外，一般都要用数值法求解。

如图2所示，现有的数值法计算流程一般分为两个运算过程。第一个运算 过程是：先假定屈曲荷载之值为P，计算目的是要确定截面的弹性单元面积和屈 服区单元面积，并判断荷载计算值与给定值是否一致，如果满足条件即进入下 一环节；第二运算过程是：先算出弹性区的惯性矩，根据弯曲屈曲(扭转屈曲 或弯扭屈曲)平衡方程得到屈曲荷载F，再检验F与前面假定的P是否吻合，如 果误差较大，可以用平均值P＝(P+F)/2作为第二轮计算的初始值，重复前面的 计算步骤，最终达到收敛要求。

然而，一方面，由于现有的数值计算法包含两个嵌套循环，在计算中需要 反复迭代，因而过程复杂，迭代循环多，数据处理量大且收敛性较差；另一方 面，由于现有的数值计算法的最终目的是为了得到弹塑性屈曲荷载，而计算中 对荷载的调整依靠式P＝(P+F)/2调整速度较慢，因而导致现有程序的运行时间 过长，不利于快速计算出结果，使得数值模型在工程中的应用受到很大的限制。 再一方面，现有的数值计算法没有考虑截面非对称因素对轴压构件的影响，因 而通过该方法计算得到的弹塑性屈曲荷载的值不准确。

因此，有必要提供一种轴心受压构件的弹塑性屈曲荷载的计算方法来克服 上述缺陷。

发明内容

本发明的目的是提供一种计算量小、计算效率高、速度快、应用范围广且 准确度高的轴心受压构件的弹塑性屈曲荷载的计算方法。

为了实现上述目的，本发明提供了一种轴心受压构件的弹塑性屈曲荷载的 计算方法包括以下步骤：

一种轴心受压构件的弹塑性屈曲荷载弦截法，所述方法包括以下步骤：

(1)判断轴心受压构件的截面信息，并根据截面信息将轴心受压构件的截 面划分为多个单元，形成截面网格；

(2)判断轴心受压构件的屈服类型，屈服类型包括弯曲屈曲、扭转屈曲和 弯扭屈曲；

(3)根据判断的屈服类型，计算对应的弦截区间点的值、屈曲荷载初始值 上的屈曲函数值和弦截区间点上的屈曲函数值；

(4)计算割线点的值；

(5)根据收敛条件其中，δ为小常数， 判断迭代是否结束，如果是，则将割线点的值作为弹性屈曲荷载输出，否则， 重新设定屈曲荷载初始值和弦截区间点的值并返回步骤(3)。

优选地，截面网格中各个单元的面积不同。

优选地，所述步骤(3)具体为：

如果判断轴心受压构件的屈服类型为弯曲屈曲，则根据公式(1)-(3)分 别计算弦截区间点的值P₁、屈曲荷载初始点上的屈曲函数值F(P₀)和弦截区间点 上的屈曲函数值F(P₁)：

P₁＝(∑σ_i(P₀)A_i+π²EI_e(P₀)/l²)/2  (1)

F(P₀)＝∑σ_i(P₀)A_i-π²EI_e(P₀)/l²  (2)

F(P₁)＝∑σ_i(P₁)A_i-π²EI_e(P₁)/l²  (3)

其中，P₁为弦截区间点的值，P₀为给定的屈曲荷载初始点的值，σ_i(P₀)为初 始荷载P₀下的第i个截面单元的应力值，σ_i(P₁)为荷载P₁下的第i个截面单元的 应力值，E为轴心受压构件的弹性模量，I_e(P₀)和I_e(P₁)分别为荷载P₀、P₁下轴心 受压构件的截面处于弹性状态区域的抗弯惯性矩，l为轴心受压构件的长度，A_i为第i个单元的面积，F(P₀)为屈曲荷载初始点上的屈曲函数值，F(P₁)为弦截区 间点上的屈曲函数值；

如果判断轴心受压构件的屈服类型为扭转屈曲，则根据公式(4)-(6)分 别计算弦截区间点的值P₁、屈曲荷载初始点上的屈曲函数值F(P₀)和弦截区间点 上的屈曲函数值F(P₁)：

P₁＝[π²EI_ω/l²+GI_et(P₀)+G_tI_pt(P₀)+K(P₀)]/2  (4)

F(P₀)＝π²EI_ω/l²+GI_et(P₀)+G_tI_pt(P₀)-K(P₀)  (5)

F(P₁)＝π²EI_ω/l²+GI_et(P₁)+G_tI_pt(P₁)-K(P₁)  (6)

其中，P₁为弦截区间点的值，P₀为给定的屈曲荷载初始点的值，EI_ω为翘曲 刚度，l为轴心受压构件的长度，G为剪切模量，I_et(P₀)和I_et(P₁)分别为荷载P₀、P₁下轴心受压构件的截面处于弹性状态区域的抗扭惯性矩，I_pt(P₀)和I_pt(P₁)分别为 荷载P₀、P₁下轴心受压构件的截面处于塑性状态区域的抗扭惯性矩，K(P₀)和 K(P₁)分别为荷载P₀、P₁下的Wagner系数,G_t为剪变模量，F(P₀)为屈曲荷载初始 点上的屈曲函数值，F(P₁)为弦截区间点上的屈曲函数值；

如果判断轴心受压构件的屈服类型为弯扭屈曲，则根据公式(7)-(9)计算弦 截区间点的值P₁、屈曲荷载初始点上的屈曲函数值F(P₀)和弦截区间点上的屈曲 函数值F(P₁):

$P_1=\left\{\right[{\pi }^2{\mathrm{EI}}_{\mathrm{ey}}\left(P_0\right)/l^2-P_0\left]\right[{\pi }^2{\mathrm{EI}}_{\omega }/l^2+{\mathrm{GI}}_{\mathrm{et}}\left(P_0\right)+G_t{I}_{\mathrm{pt}}\left(P_0\right)-K\left(P_0\right)]+{P_0}^2{y}_0^2\}/4---\left(7\right)$

$F\left(P_0\right)=[{\pi }^2{\mathrm{EI}}_{\mathrm{ey}}\left(P_0\right)/l^2-P_0][{\pi }^2{\mathrm{EI}}_{\omega }/l^2+{\mathrm{GI}}_{\mathrm{et}}\left(P_0\right)+G_t{I}_{\mathrm{pt}}\left(P_0\right)-K\left(P_0\right)]-{P_0}^2{y}_0^2---\left(8\right)$

F(P₁)＝[π²EI_ey(P₁)/l²-P₁][π²EI_ω/l²+GI_et(P₁)+G_tI_pt(P₁)-K(P₁)]-P₁²y₀²  (9)

其中，P₁为弦截区间点的值，P₀为给定的屈曲荷载初始点的值，E为弹性模 量，G_t为剪变模量，y₀为剪心距，I_ey(P₀)和I_ey(P₁)分别为荷载P₀、P₁下弹性区域关 于y轴的惯性矩，l为轴心受压构件的长度，G为剪切模量，I_et(P₀)和I_et(P₁)分别 为荷载P₀、P₁下轴心受压构件的截面处于弹性状态区域的抗弯惯性矩，I_pt(P₀)和 I_pt(P₁)分别为荷载P₀、P₁下轴心受压构件的截面处于塑性状态区域的抗扭惯性矩， K(P₀)和K(P₁)分别为荷载P₀、P₁下的Wagner系数，EI_ω为翘曲刚度，F(P₀)为屈曲 荷载初始点上的屈曲函数值，F(P₁)为弦截区间点上的屈曲函数值。

优选地，所述步骤(4)具体为：

根据公式计算割线点的值，其中，P₂为割线 点的值，P₀为给定的屈曲荷载初始点的值，P₁为弦截区间点的值，F(P₀)为屈曲 荷载初始点上的屈曲函数值，F(P₁)为弦截区间点上的屈曲函数值。

优选地，截面信息包括截面形状、宽度和厚度。

与现有技术相比，一方面，本发明的计算方法简化了截面的塑性平衡分析， 较现有方法少了一个迭代循环，无需嵌套循环计算，大大减少了计算量。另一 方面，本发明的计算方法采用图3所示的弦截法，也称割线法，用解曲线上的 两点确定的直线来近似曲线求得方程的近似根且收敛阶数为1.618，虽然较二阶 收敛的牛顿法为低，但其不用计算函数的导数，且每个迭代步只需计算一次函 数值，计算量小，因此计算效率较高，适于进行非线性方程求根。再一方面， 本发明的计算方法提供了对于短柱的弹塑性屈曲快速计算方法，为以后两端铰 接轴心受压构件的快速计算提供了基础，而且本发明的方法可应用于H型、T型、 L型、十字型等不同截面型式构件的轴压屈曲荷载的计算。

通过以下的描述并结合附图，本发明将变得更加清晰，这些附图用于解释 本发明的实施例。

附图说明

图1(a)为轴心受压构件扭转屈曲变形时的示意图。

图1(b)为轴心受压构件弯扭屈曲变形时的示意图。

图2为现有轴心受压构件屈曲荷载的计算方法的流程图。

图3为弦截法的示意图。

图4为本发明轴心受压构件的弹塑性屈曲荷载的计算方法的流程图。

具体实施方式

现在参考附图描述本发明的实施例，附图中类似的元件标号代表类似的元 件。

如图4所示，本实施例的轴心受压构件的弹塑性屈曲荷载的计算方法包括 以下步骤：

(1)判断轴心受压构件的截面信息，截面信息包括截面形状、宽度和厚度， 并根据截面信息将轴心受压构件的截面划分为多个单元，形成截面网格；

(2)判断轴心受压构件的屈服类型，屈服类型包括弯曲屈曲、扭转屈曲和 弯扭屈曲；

(3)如果判断轴心受压构件的屈服类型为弯曲屈曲，则根据公式(1)-(3) 分别计算弦截区间点的值P₁、屈曲荷载初始点上的屈曲函数值F(P₀)和弦截区间 点上的屈曲函数值F(P₁)：

P₁＝(∑σ_i(P₀)A_i+π²EI_e(P₀)/l²)/2  (1)

F(P₀)＝∑σ_i(P₀)A_i-π²EI_e(P₀)/l²  (2)

F(P₁)＝∑σ_i(P₁)A_i-π²EI_e(P₁)/l²  (3)

其中，P₁为弦截区间点的值，P₀为给定的屈曲荷载初始点的值，σ_i(P₀)为初 始荷载P₀下的第i个截面单元的应力值，σ_i(P₁)为荷载P₁下的第i个截面单元的 应力值，E为轴心受压构件的弹性模量，I_e(P₀)和I_e(P₁)分别为荷载P₀、P₁下轴心 受压构件的截面处于弹性状态区域的抗弯惯性矩，l为轴心受压构件的长度，A_i为第i个单元的面积，F(P₀)为屈曲荷载初始点上的屈曲函数值，F(P₁)为弦截区 间点上的屈曲函数值；

如果判断轴心受压构件的屈服类型为扭转屈曲，则根据公式(4)-(6)分 别计算弦截区间点的值P₁、屈曲荷载初始点上的屈曲函数值F(P₀)和弦截区间点 上的屈曲函数值F(P₁)：

P₁＝[π²EI_ω/l²+GI_et(P₀)+G_tI_pt(P₀)+K(P₀)]/2  (4)

F(P₀)＝π²EI_ω/l²+GI_et(P₀)+G_tI_pt(P₀)-K(P₀)  (5)

F(P₁)＝π²EI_ω/l²+GI_et(P₁)+G_tI_pt(P₁)-K(P₁)  (6)

其中，P₁为弦截区间点的值，P₀为给定的屈曲荷载初始点的值，EI_ω为翘曲 刚度，l为轴心受压构件的长度，G为剪切模量，I_et(P₀)和I_et(P₁)分别为荷载P₀、P₁下轴心受压构件的截面处于弹性状态区域的抗扭惯性矩，I_pt(P₀)和I_pt(P₁)分别为 荷载P₀、P₁下轴心受压构件的截面处于塑性状态区域的抗扭惯性矩，K(P₀)和 K(P₁)分别为荷载P₀、P₁下的Wagner系数,G_t为剪变模量，F(P₀)为屈曲荷载初始 点上的屈曲函数值，F(P₁)为弦截区间点上的屈曲函数值；

如果判断轴心受压构件的屈服类型为弯扭屈曲，则根据公式(7)-(9)计算 弦截区间点的值P₁、屈曲荷载初始点上的屈曲函数值F(P₀)和弦截区间点上的屈 曲函数值F(P₁):

$P_1=\left\{\right[{\pi }^2{\mathrm{EI}}_{\mathrm{ey}}\left(P_0\right)/l^2-P_0\left]\right[{\pi }^2{\mathrm{EI}}_{\omega }/l^2+{\mathrm{GI}}_{\mathrm{et}}\left(P_0\right)+G_t{I}_{\mathrm{pt}}\left(P_0\right)-K\left(P_0\right)]+{P_0}^2{y}_0^2\}/4---\left(7\right)$

$F\left(P_0\right)=[{\pi }^2{\mathrm{EI}}_{\mathrm{ey}}\left(P_0\right)/l^2-P_0][{\pi }^2{\mathrm{EI}}_{\omega }/l^2+{\mathrm{GI}}_{\mathrm{et}}\left(P_0\right)+G_t{I}_{\mathrm{pt}}\left(P_0\right)-K\left(P_0\right)]-{P_0}^2{y}_0^2---\left(8\right)$

F(P₁)＝[π²EI_ey(P₁)/l²-P₁][π²EI_ω/l²+GI_et(P₁)+G_tI_pt(P₁)-K(P₁)]-P₁²y₀²  (9)

其中，P₁为弦截区间点的值，P₀为给定的屈曲荷载初始点的值，E为弹性模 量，G_t为剪变模量，y₀为剪心距，I_ey(P₀)和I_ey(P₁)分别为荷载P₀、P₁下弹性区域关 于y轴的惯性矩，l为轴心受压构件的长度，G为剪切模量，I_et(P₀)和I_et(P₁)分别 为荷载P₀、P₁下轴心受压构件的截面处于弹性状态区域的抗弯惯性矩，I_pt(P₀)和 I_pt(P₁)分别为荷载P₀、P₁下轴心受压构件的截面处于塑性状态区域的抗扭惯性矩， K(P₀)和K(P₁)分别为荷载P₀、P₁下的Wagner系数，EI_ω为翘曲刚度，F(P₀)为屈曲 荷载初始点上的屈曲函数值，F(P₁)为弦截区间点上的屈曲函数值；

(4)根据公式计算割线点的值，其中，P₂为 割线点的值，P₀为给定的屈曲荷载初始点的值，P₁为弦截区间点的值，F(P₀)为 屈曲荷载初始点上的屈曲函数值，F(P₁)为弦截区间点上的屈曲函数值；

(5)根据收敛条件|P₂-P₁|/P₁＜δ，其中，P₂为割线点的值，P₁为弦截区间点 的值，δ为小常数，判断迭代是否结束，如果是，则将割线点的值作为弹塑性屈 曲荷载输出，否则，重新设定屈曲荷载初始值和弦截区间点的值并返回步骤(3)。

较佳地，δ计算上根据需要可取为1×10^-3或1×10^-6等小常数。重新设定屈 曲荷载初始值和弦截区间点的值时，令屈曲荷载初始值等于弦截区间点的值， 而弦截区间点的值等于下一弦截区间点的值。

详细地，当轴心受压构件的屈服类型为弯曲屈曲时，轴心受压构件的弹塑 性屈曲荷载的推导过程如下：

假设截面的残余应力分布为σ_r。

因为截面各点的应力不同，在计算时可将截面划分许多单元，不同区域(如 翼缘和腹板等)的单元面积可不相同。以A_i表示任一单元的面积，在压力P的作 用下，此单元的应变和应力分别为ε_i和σ_i，该单元的残余应力为σ_ri。

因P值之不同，此单元的应力σ_i可能小于屈服强度σ_y，处在弹性范围；也 可能等于σ_y，即已经屈服。为了确定σ_i之值，可以先从计算单元的应变ε_i着手。 ε_i是外力作用的轴向应变ε₀和残余应变ε_ri＝σ_ri/E之代数和。为便于运算，可取 压应力为正，拉应力为负。

计算时假定材料为理想的弹塑性体，屈服强度σ_y。这样截面中任一单元的 应变为

ε_i＝ε₀+ε_ri  (11)

当-σ_y/E≤ε_i≤σ_y/E时，

σ_i＝Eε₀+σ_ri

当ε_i＞σ_y/E时，

σ_i＝σ_y

当ε_i＜-σ_y/E时，

σ_i＝-σ_y

短柱的轴向压力为

P＝∑σ_iA_i

短柱的截面积为A＝∑A_i，其中包括了部分已经屈服的单元面积。如果以A_ei表示弹性状态的单元面积，则总的弹性区的截面积为A_e＝∑A_ei。

本发明所涉及是两端铰接的轴心受压构件，按照切线模量理论求解。

当构件有微小弯曲时，截面任一点纤维的应变ε_i是轴向应变ε₀、残余应变ε_ri和弯曲应变z_iΦ的三项代数和；z_i是单元面积A_i之中心至弯曲轴的距离。Φ是曲 率，Φ＝-y″。如单元面积的应力在弹性状态，则其面积用A_ei表示。这样截面上 任一点的应变为

ε_i＝ε₀+ε_ri-z_iy″  (12)

在临界条件下-z_iy″与ε₀+ε_ri相比是一个小量，故仍可以式(11)来确定单元 的应变和应力，从而断定孰是截面的弹性单元和屈服单元。

因此，构件内力矩可表示为

$\left(\begin{array}{c}{M}_i={\int }_A\mathrm{\sigma zdA}={\int }_{A_e}E({ϵ}_0+{ϵ}_r-{\mathrm{zy}}^{\prime \prime })\mathrm{zdA}+{\int }_{(A-A_e)}{\sigma }_y\mathrm{zdA}\\ ={\int }_{A_e}E({ϵ}_0+{ϵ}_r)\mathrm{zdA}-{\int }_{A_e}{\mathrm{Ez}}^2{y}^{\prime \prime }\mathrm{dA}+{\sigma }_y{\int }_{(A-A_e)}\mathrm{zdA}\\ =-y^{\prime \prime }·E{\int }_{A_e}{z}^2\mathrm{dA}+E({ϵ}_0+{ϵ}_r)·{\int }_{A_e}\mathrm{zdA}+{\sigma }_y{\int }_{(A-A_e)}\mathrm{zdA}\\ =-{\mathrm{EI}}_e·y^{\prime \prime }+M_{\mathrm{ep}}\left(P\right)\end{array}\right)---\left(13\right)$

其中$I_e={\int }_{A_e}{z}^2\mathrm{dA},M_{\mathrm{ep}}=E({ϵ}_0+{ϵ}_r)·{\int }_{A_e}\mathrm{zdA}+{\sigma }_y{\int }_{(A-A_e)}\mathrm{zdA}.$

从上式中可以看出与曲率相关的项为弹性区域的惯性矩，M_ep与弹性区轴向 荷载及弹塑性区域的大小相关，是轴压P的函数。当构件截面及残余应力分布 对称时，可以得出在截面的塑性区域有M_ep＝0。但对于截面或残余应力分布非对 称的情况，此式并不成立。

而外力矩可以表示为

M_e＝P·y  (14)

则内外力矩的平衡方程为

EI_e·y″+Py-M_ep＝0  (15)

此时方程与受初偏心影响的轴心受压构件的屈曲平衡方程形式相同。

令k²＝P/EI_e，可以得到通解

$y=-(\frac{1-\cos \mathrm{kl}}{\sin \mathrm{kl}}\sin \mathrm{kx}+\cos \mathrm{kx}-1)\frac{M_{\mathrm{ep}}}{P}---\left(16\right)$

构件的最大挠度为

$y\left(\frac{l}{2}\right)=y_{\max }=-(\sec \frac{\mathrm{kl}}{2}-1)\frac{M_{\mathrm{ep}}}{P}---\left(17\right)$

由此可以得出构件的荷载-挠度曲线。此时构件不出现分岔失稳，而出现极 值点失稳现象。由于此类问题属于压弯构件的极限荷载计算，不属于本发明所 涉及内容。

当截面及残余应力对称时，平衡方程为

EI_ey″+Py＝0  (18)

解此方程可得到弹塑性屈曲荷载

$P_{\mathrm{cr}}=\frac{{\pi }^2{\mathrm{EI}}_e}{l^2}---\left(19\right)$

从上式可以得到另一个结论：残余应力将使构件提前进入弹塑性状态，降 低构件的屈曲载荷，构件截面的有效惯性矩只是截面弹性区的惯性矩I_e，而构件 的抗弯刚度由EI降至EI_e。

根据以上各式可以认为截面各单元的应变、应力及弹性区的惯性矩均可表 示为轴向荷载P的函数，分别为ε_i(P)，σ_i(P)和I_e(P)。因此本问题可表示为以下 非线性方程的求解

F(P)＝∑σ_i(P)A_i-π²EI_e(P)/l²＝0  (20)

该非线性方程的求解主要有二分法及牛顿迭代法。但由于应力-应变间弹塑 性本构关系的影响，该方程不能表示为等价形式P＝g(P)，因此不能使用不动点 迭代法求解；同样也不能得到函数F(P)的导数F′(P)。因此使用牛顿迭代法中的 弦截法进行数值计算。

详细地，当轴心受压构件的屈服类型为扭转屈曲时，轴心受压构件的弹塑 性屈曲荷载的推导过程如下：

一般双轴对称截面轴心受压构件，可能绕截面的两个对称轴发生弯曲屈曲， 但是有些抗扭刚度和抗翘曲刚度很弱的轴心受压构件，如双轴对称的十字型截 面轴心受压构件，除了有可能发生绕对称轴x或y的弯曲屈曲外，还有可能发 生绕截面纵轴z转动的扭转屈曲，此时纵轴本身不发生弯曲变形，只是截面绕 纵轴旋转一个角度。

考虑具有残余应力分布σ_r的轴心受压构件，如果按照弹性公式得到的扭转 屈曲应力σ_ω超过了构件的有效比例极限，此时构件将在弹塑性状态屈曲。计算 时，假定材料为理想弹塑性，在屈服区，变形模量取E_t＝0，而剪变模量可取G_t＝G/4。

将截面划分为n个单元，每个单元面积为A_i。这样截面中任一单元的应变 为

ε_i＝ε₀+ε_ri＝ε₀+σ_ri/E  (21)

当-σ_y/E≤ε_i≤σ_y/E时，

σ_i＝Eε₀+σ_ri，G_t＝G。

当ε_i＞σ_y/E时，

σ_i＝σ_y，G_t＝G/4。

当ε_i＜-σ_y/E时，

σ_i＝-σ_y，G_t＝G/4。

记处于弹性阶段的单元为A_ei，处于弹塑性阶段的单元为A_pi。

构件的截面压力为

P＝∑σ_iA_i  (22)

构件任意截面的扭转屈曲控制方程为

其中EI_ω为翘曲刚度；I_et为截面弹性区的抗扭惯性矩，I_pt为屈服区的抗扭惯性 矩；Wagner效应系数为K，即截面内单位应力绕截面剪心形成的扭矩。

K＝∫_Aσ[(x-x_s)²+(y-y_s)²]dA＝∑σ_i[(x_i-x_s)²+(y_i-y_s)²]A_i  (24)

对于两端简支的轴压构件，可得到第一阶扭转屈曲模态为代 入(23)可得

$K=\frac{{\pi }^2}{l^2}{\mathrm{EI}}_{\omega }+{\mathrm{GI}}_{\mathrm{et}}+G_t{I}_{\mathrm{pt}}---\left(25\right)$

根据以上各式，可以认为截面各单元的应力、Wagner效应系数均可表示为 轴向荷载P的函数，分别为σ_i(P)和K(P)。因此本问题可表示为以下非线性方程 的求解。

$F\left(P\right)=\frac{{\pi }^2}{l^2}{\mathrm{EI}}_{\omega }+{\mathrm{GI}}_{\mathrm{et}}\left(P\right)+G_t{I}_{\mathrm{pt}}\left(P\right)-K\left(P\right)=0---\left(26\right)$

该非线性方程的求解主要有二分法及牛顿迭代法。但由于应力-应变间弹塑 性本构关系的影响，计算函数F(P)的导数F′(P)非常困难，因此使用牛顿迭代法 中的弦截法进行数值计算。

详细地，当轴心受压构件的屈服类型为弯扭屈曲时，轴心受压构件的弹塑 性屈曲荷载的推导过程如下：

对于单轴对称截面轴心受压构件，如等边角钢轴心受压构件，除了可能绕 截面的非对称轴x发生弯曲屈曲外，还可能在绕截面的对称轴y弯曲的同时又 绕通过截面剪心S的纵轴扭转而发生弯扭屈曲。截面的位移有侧移u和扭转角如果截面不具有对称轴，如不等边角钢截面，截面的形心和剪心不重合，这种 截面的轴心受压构件只可能发生弯扭屈曲。此时截面不仅绕两个主轴x和y弯 曲，而且还绕通过剪心的纵轴扭转。

对于单轴对称轴压构件，其弹性屈曲方程为

其中y₀为剪心距。

对于弹塑性问题，当截面形状及残余应力分布对于弯扭轴对称时，其屈曲 控制方程可表示为

其中I_ey为弹性区域关于y轴的惯性矩，P＝∫_AσdA＝∑σ_iA_i，

K＝∫_Aσ[(x-x_s)²+(y-y_s)²]dA＝∑σ_i[(x_i-x_s)²+(y_i-y_s)²]A_i。

对于两端简支轴心受压构件的弯扭屈曲，根据边界条件可以确定其屈曲模 态为u＝C₁sin(πz/l)和代入屈曲方程，可得

$\left(\begin{array}{c}({\pi }^2{\mathrm{EI}}_y/l^2-P)C_1+{\mathrm{Py}}_0{C}_2=0\\ -{\mathrm{Py}}_0{C}_1+(K-{\pi }^2{\mathrm{EI}}_{\omega }/l^2-{\mathrm{GI}}_{\mathrm{et}}-G_t{I}_{\mathrm{pt}})C_2=0\end{array}\right)---\left(29\right)$

C₁和C₂有非零解的条件为其系数行列式为零，即

$({\pi }^2{\mathrm{EI}}_{\mathrm{ey}}/l^2-P)({\pi }^2{\mathrm{EI}}_{\omega }/l^2+{\mathrm{GI}}_{\mathrm{et}}+G_t{I}_{\mathrm{pt}}-K)-P^2{y}_0^2=0---\left(30\right)$

同样可以认为截面各单元的应力、Wagner效应系数、惯性矩和扭转矩均可 表示为轴向荷载P的函数，分别为σ_i(P)、K(P)和I_ey(P)、I_et(P)。因此本问题可 表示为以下非线性方程的求解。

$F\left(P\right)=[{\pi }^2{\mathrm{EI}}_{\mathrm{ey}}\left(P\right)/l^2-P][{\pi }^2{\mathrm{EI}}_{\omega }/l^2+{\mathrm{GI}}_{\mathrm{et}}\left(P\right)+G_t{I}_{\mathrm{pt}}\left(P\right)-K\left(P\right)]-P^2{y}_0^2=0---\left(31\right)$

由于截面存在台阶式的残余应力分布，随着轴向压力的增大，上述非线性 方程会出现突变，造成方程导数的不连续。因此使用牛顿迭代法中的弦截法进 行数值计算。

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限 制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人 员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未 脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利 要求范围当中。