拉压杆模型(strut-and-tie models，STM)是一种广泛使用的钢筋混凝土结构塑性设计方法，是类桁架模型的一般化。由塑性理论可知，STM的理论依据是极限分析下限定理[1]。不同于传统的截面设计方法，该方法将受压区的混凝土抽象为压杆，将受拉区的钢筋抽象为拉杆，将整个结构的力学行为抽象为由拉杆、压杆和结点组成的拉压杆模型。STM的空间架构特性使得其特别适用于混凝土结构中由几何或静力不连续所产生的D(discontinuity)区域的设计，但Schlaich指出STM经过适当调整亦可适用于符合伯努力假定的B(Bernoulli)区域的设计[2]。
混凝土结构中梁柱节点就是一个典型的D区域。由于梁柱节点复杂的应力分布，传统的弹性应力迹线法或荷载路径法[3]很难生成合适的拉压杆模型，且由于设计者工程经验的差异，可能生成多种形式的STM。随着拓扑优化方法[4-5]的不断发展，基于拓扑优化的拉压杆模型逐渐成为一个热门的研究方向。最初基结构法[6]被用来产生拉压杆模型，但是该方法可能产生局部或者奇异最优解。此后基于性能的优化方法(PBO)[7]和固体各向同性惩罚方法(SIMP)[8]等其他连续体拓扑优化方法也相继被用来产生拉压杆模型。然而，在上述两种拓扑优化方法中，拓扑优化结果的边界不能用表达式显式地描述，即边界采用基于单元的黑白像素的方式表示，因此具有隐式的特性。这会造成以下两方面的问题：一方面，锯齿状的模糊边界难于形成明确的边界，从而不利于STM的形成；另一方面，隐式的拓扑优化方法中，优化变量数与有限单元数相等。较多的优化变量导致优化效率较低，这一点在三维构件中尤为明显。
2014年，大连理工大学郭旭提出了基于可移动变形组件(moving morphable components，MMC)的拓扑优化方法[9]。该方法采用拓扑描述函数来显式表达优化结果的边界，并通过一组可移动、变形、交叉和重叠的组件来表示结构的拓扑。本质上，MMC拓扑优化方法将结构的拓扑优化转化为组件的形状优化和尺寸优化。近年来，该方法与等几何分析、增材制造技术和比例边界元相结合，获得了较大的成功。作为一种显式的拓扑优化方法，该方法具有边界相对光滑、设计变量数较少、不存在灰度单元等优点。ZHANG et al[10]提出了一种基于MMC的多种材料的拓扑优化方法，并通过一组二维和三维的算例验证了该方法的有效性和高效性。SUN et al[11]用一组MMC来模拟复合加筋板的加劲肋，将基于MMC的拓扑优化方法用于加劲肋的布局优化设计。CUI et al[12]采用符合基尔霍夫板理论的板单元，研究了基于MMC 拓扑优化方法的薄板结构设计方法。
综上所述，基于MMC的拓扑优化方法虽然在多个领域获得了成功应用，但是在混凝土结构中的应用研究较少。鉴于组件和拉压杆的相似性，本文将显式的拓扑优化方法与钢筋混凝土结构中的拉压杆模型相结合，以期准确、高效地生成可靠的STM.此外，在实际工程中，结构上的多种荷载往往不是同时作用的，因此多荷载工况的分析尤为重要。在文献[13-14]的基础上，本文通过3种不同荷载作用下混凝土结构梁柱节点算例，研究了基于MMC的拓扑优化的拉压杆模型，考虑了单荷载和多荷载两种工况，证实了该方法的可行性、有效性和高效性。
1 基于MMC拓扑优化方法
1.1 结构的拓扑描述函数
先给出一个组件的拓扑描述函数为：
(1)
式中：x表示设计域D0中点的位置坐标；χi(x)表示第i个组件的拓扑描述函数，Ωi和∂Ωi分别是第i个组件的内部区域和边界。
为了简化，本文仅讨论平面问题。骨架为直线型且厚度呈二次函数变化的组件如图1，组件更多的几何形式参见文献[14].
图1 单个组件的几何描述
Fig.1 Geometric description of a single component
从而可得到χi(x)的具体表达式：
(2)
其中
(3)
式中：x和x'分别为设计域内在整体和局部坐标系下点的横坐标；y为设计域内在整体坐标系下点的纵坐标；x0i和y0i分别为第i个组件在整体坐标系下中心点的横、纵坐标；Li和θi分别为第i个组件长轴的一半和从x轴正方向逆时针旋转的倾斜角；参数p取为6.
从而χS(x)可由下式表达：
(4)
其中
χS(x)=max(χ1(x),…,χnc(x)) .
(5)
式中：nc表示所有组件的个数，ΩS是所有组件所占的区域。
因此，列向量
可以惟一表达结构的拓扑，其中di=(x0i,y0i,Li,h1,h2,h3,θi)T.
1.2 问题列式
按照惯例，将产生最优拉压杆模型转换成求解一个体积约束下最小柔度优化问题。
寻找
最小化
使得
Kui-Fi=0 ,
d⊆Ud,
(6)
式中：C和V分别为该问题的目标函数和约束函数；Fi和ui分别为第i个荷载列向量和位移列向量，K为结构的整体刚度矩阵；M是参与荷载组合的个数，对于单荷载工况为1，对于多荷载工况为大于1的整数；αi为柔度的组合系数，对于单荷载工况取1，多荷载工况根据实际情况选取；V0和
分别是设计域总体积和给定的体积百分比上限；
表示狄利克雷边界Γui上给定的位移；Ud为设计变量d的容许集。
2 数值实现
2.1 有限单元法
有限单元分析中，采用四结点四边形单元求解结构在荷载作用下的位移响应，弹性模量E采用混凝土的弹性模量。混凝土的虚拟杨氏模量Ee为：
(7)
其中正则化函数H(x)为：
H(x)=
(8)
式中：
为结构拓扑描述函数在单元e的四个结点处的值；q为大于1的幂指数，本文中取2；α=0.001为一个非0的较小的正数；ω为正则化参数，其取值与单元大小有关。
式(6)中的表达式可进一步整理为：
(9)
(10)
(11)
式中：
为采用混凝土材料的单元刚度矩阵；ue为单元位移列向量；NE为有限单元的总数。
2.2 设计灵敏度
本文采用移动渐近线法[15]作为最优化算法。该算法需要提供灵敏度信息。这里给出目标函数和不等式约束函数的灵敏度，分别为：
(12)
(13)
式中：a表示任意一个设计变量。
3 数值算例
本节给出MMC拓扑优化的混凝土结构梁柱节点在单荷载和多荷载两种工况下的算例，说明用MMC拓扑优化方法产生拉压杆模型的可行性和有效性。单荷载工况考虑了两种不同的受力情况，分别为梁柱节点仅受弯矩作用和受弯矩和剪力共同作用；多荷载工况为梁柱节点受到弯矩、剪力和轴力的共同作用。梁柱节点的几何尺寸见图2.在所有工况中，梁柱节点初始的组件布置都是相同的，均采用24个组件，见图3.本文的工作是在MATLAB优化代码[13]基础上进行的。根据拉压杆模型的通用做法，粗实线代表拉杆，粗虚线代表压杆。
图2 梁柱节点几何尺寸
Fig.2 Geometric dimensions of beam-column connections
图3 梁柱节点初始组件布置
Fig.3 Initial layout of components in beam-column connections
为了与其他研究者的工作进行比较，混凝土的杨氏模量取为28 567 MPa，泊松比为0.15.为了研究的简化，梁和柱的厚度均取为300 mm，体积百分比的上限取为0.4.有限元分析中均采用四结点平面正方形单元，单元的边长为25 mm.
3.1 梁柱节点仅受弯矩作用
梁柱节点仅受弯矩作用的情况属于单荷载工况，其所受荷载为M=40 kN·m，FV=0 kN，FN=0 kN.弯矩的施加方式为在距梁上下表面各50 mm处施加一对水平向右和向左的力，大小为100 kN.受弯矩作用的梁柱节点拓扑过程和相应的拉压杆模型见图4。经过239次迭代后，达到最优拓扑，最小柔度目标函数值为51.68 J.
在优化过程中，采用了24个组件，设计变量共168个。由图4可知，优化结果中边界并非呈现理想的直线，一方面是由于梁柱节点区存在凹角，即所谓的几何不连续；另一方面由于计算机绘图的原因，在组件相交处会存在不同程度的弯曲。柔度目标和体积约束函数的收敛过程见图5.荷载通过受拉和受压两条路径传递到固端支座，在节点区有两个拉杆连接这两条荷载路径。
为了进行直观的比较，采用了SIMP的拓扑优化方法对该算例进行了研究。图6为SIMP方法的最优拓扑结果。在优化结果中存在灰度单元，而且边界是锯齿状的。表1为一个典型迭代步MMC和SIMP方法的CPU时间比较。在传统的SIMP方法中，设计变量数与有限单元数量相同，即1 360个。由表1可知，两种方法中有限元分析和灵敏度分析所用时间相差不大，MMA所用时间则相差较大。可见MMC方法使得优化求解效率会得到极大的提升。
图4 仅受弯矩作用梁柱节点优化过程和拉压杆模型
Fig.4 Optimization process and STM in beam-column connections under bending moment
图5 目标和约束函数收敛过程曲线
Fig.5 Convergence history curves of objective and constraint functions
图6 SIMP方法最优拓扑结果
Fig.6 Optimal topology by SIMP method
表1 一个典型迭代步MMC和SIMP方法CPU时间比较
Table 1 A comparison of CPU time for a typical iterative step between MMC and SIMP approaches
方法CPU时间/s有限元分析灵敏度分析MMA合计SIMP0.169 30.010 79.582 29.762 2MMC0.216 10.028 20.001 90.246 2
3.2 梁柱节点受弯矩和剪力共同作用
考虑到梁端在竖向荷载下可以产生弯矩和剪力，梁柱节点受弯矩和剪力共同作用仍可按单荷载工况考虑，其所受荷载为M=40 kN·m，FV=80 kN，FN=0 kN.剪力施加位置同3.1节中力偶的施加位置，施加方式为竖直向上各40 kN.受弯矩和剪力共同作用的梁柱节点拓扑过程和相应的拉压杆模型见图7.经过693次迭代后，达到最优拓扑，最小柔度目标函数值为284.80 J.
图7 受弯矩和剪力共用作用梁柱节点优化过程和拉压杆模型
Fig.7 Optimization process and STM in beam-column connections under bending moment and shear force
随着荷载数目的增多，结构中应力分布更加复杂，迭代所需次数也增多。与3.1节不同的是，剪力的存在，使得结构受力机理也发生变化。在拉压杆模型中表现为增加了一个倾斜的拉杆来抵抗竖向的剪力。
3.3 梁柱节点受弯矩、剪力和轴力共同作用
考虑到梁柱节点在水平荷载作用下会产生轴向力，梁柱节点受弯矩、剪力和轴力共同作用的情况可视为多荷载工况，其所受荷载为M=40 kN·m，FV=80 kN，FN=80 kN.为了简化，柔度组合系数取1.弯矩和剪力的施加方法同3.2节，轴力的施加位置同3.2节中剪力的施加位置，施加方式为水平向右各40 kN.受弯矩、剪力和轴力共同作用的梁柱节点拓扑过程和相应的拉压杆模型分别见图8和图9.经过396次迭代后，达到最优拓扑，最小柔度目标函数值为366.62 J.
图8 受弯矩、剪力和轴力共用作用梁柱节点优化过程
Fig.8 Optimization process in beam-column connections under bending moment, shear force, and axial force
图9 受弯矩、剪力和轴力共用作用梁柱节点拉压杆模型
Fig.9 Strut-and-tie models of beam-column connections under bending moment, shear force, and axial force
为了使用的方便，在形成拉压杆模型过程中对模型的结点位置进行了一定程度的简化。图9(b)中数字0表示该杆件所受轴力为0.由于轴力的存在，荷载传递路径有所不同。相比于3.2节的STM，用位于柱上端一个斜向的拉杆来承受由轴力引起的柱中剪切应力。
4 结论
将基于MMC的拓扑优化方法应用到混凝土结构拉压杆模型中，以梁柱节点的拉压杆模型为例研究了3种不同荷载情况下的节点拉压杆模型，考察了这一方法的可行性、有效性和高效性，得到以下结论：
1) 作为一种基于构件的拓扑优化方法，基于MMC的拓扑优化方法从机理上适合于形成可靠的梁柱节点拉压杆模型。
2) 由于显式的特性，基于MMC拓扑优化的结果克服了传统的隐式拓扑优化的缺点，如灰度单元、网格依赖性、锯齿的边界等。此外，由于设计变量数的大幅减少，极大地提高了优化求解效率。
3) 基于MMC的拓扑优化方法不仅适用于单荷载工况，也适用于多荷载工况，所产生的梁柱节点拉压杆模型符合力学原理，合理可靠。
[1] DRUCKER D C.On structural concrete and the theorems of limit analysis[J].IABSE Publication,1961,21:49-59.
[2] SCHLAICH J,SCH
FER K,JENNEWEIN M.Toward a consistent design of structural concrete[J].PCI Journal,1987,32(3):75-150.
[3] HE Z Q,LIU Z,WANG J Q,et al.Development of strut-and-tie models using load path in structural concrete[J].Journal of Structural Engineering,2020,146(5):06020004.
[4] GUO X,CHENG G D.Recent development in structural design and optimization[J].Acta Mechanica Sinica,2010,26(6):807-823.
[5] 杜建镔.结构优化及其在振动和声学设计中的应用[M].北京：清华大学出版社,2015.
[6] ALI M A,WHITE R N.Automatic generation of truss model for optimal design of reinforced concrete structures[J].ACI Structural Journal,2001,98(4):431-442.
[7] LIANG Q Q,UY B,STEVEN G P.Performance-based optimization for strut-tie modeling of structural concrete[J].Journal of Structural Engineering,2002,128(6):815-823.
[8] BRUGGI M.Generating strut-and-tie patterns for reinforced concrete structures using topology optimization[J].Computers & Structures,2009,87(23-24):1483-1495.
[9] GUO X,ZHANG W S,ZHONG W L.Doing topology optimization explicitly and geometrically-a new moving morphable components based framework[J].Journal of Applied Mechanics-Transactions of the ASME,2014,81(8):081009.
[10] ZHANG W S,SONG J F,ZHOU J H,et al.Topology optimization with multiple materials via moving morphable component (MMC) method[J].International Journal for Numerical Methods in Engineering,2018,113(11):1653-1675.
[11] SUN Z,CUI R H,CUI T C,et al.An optimization approach for stiffener layout of composite stiffened panels based on moving morphable components (MMCs)[J].Acta Mechanica Solida Sinica,2020,33(5):650-662.
[12] CUI T C,SUN Z,LIU C,et al.Topology optimization of plate structures using plate element-based moving morphable component (MMC) approach[J].Acta Mechanica Sinica,2020,36(2):412-421.
[13] QIAO W Z,CHEN G R.Generation of strut-and-tie models in concrete structures by topology optimization based on moving morphable components[J].Engineering Optimization,2020,doi:10.1080/ 0305215X.2020.1781843.
[14] ZHANG W S,YUAN J,ZHANG J,et al.A new topology optimization approach based on Moving Morphable Components (MMC) and the ersatz material model[J].Structural and Multidisciplinary Optimization,2016,53(6):1243-1260.
[15] SVANBERG K.The method of moving asymptotes-a new method for structural optimization[J].International Journal for Numerical Methods in Engineering,1987,24(2):359-373.
