基于变密度法的实体重力坝拓扑优化

， ，

(西华大学 能源与环境学院，成都 610039)

1 研究背景

实体重力坝在我国在建或拟建的重力坝中所占比例最大，但因其体积庞大，使其修建成本较高，如何通过结构断面优化达到降低实体重力坝的体积是目前研究的重点。结构优化是在满足结构安全性和稳定性的前提下，使结构的材料分布更均匀合理，而结构优化中以结构拓扑优化[1]为当前研究的热点。结构拓扑优化主要分为连续性拓扑优化和离散性拓扑优化[2]，目前研究较多的为连续性拓扑优化，其主要方法有均匀化方法(Homogenization Method)[3]、变密度法(Variable Density Method)[4]、渐进结构优化法(Evolutionary structural optimization)[5]等。变密度法作为拓扑优化中的一种优化方法，假定材料的相对密度和弹性模量之间建立某种映射关系，以连续变量的函数形式将相对密度与弹性模量间的假定关系表达出来，然后通过对材料惩罚的方式使材料进行重分布，最后通过优化准则，将相对密度较小、对结构刚度影响较小的材料删除，使结构达到优化的目的。常见的密度刚度插值方式有2种：一种是SIMP(solid isotropic microstructures with penalization)插值模型；另一种是RAMP(rational approximation of material properties)插值模型。

变密度法通过引入惩罚因子对连续体结构的中间密度进行惩罚，使密度向0或1两端靠拢，形成分明的实体和空洞。但在拓扑优化设计中往往会遇到各种问题[6-7]，比如多空材料、棋盘格、网格依赖性、局部极限值等。为了解决拓扑优化存在的此类问题，本文通过改变惩罚因子，以寻找最适合的惩罚因子使优化模型符合设计要求。

2 变密度插值模型

变密度法是将0.0～1.0的离散变量通过引入连续变量X、惩罚因子p(q)等，转变为连续体优化问题。其中连续变量X取值为[0,1]，惩罚因子p(q)取值为大于1的自然数。

2.1 SIMP插值模型

对于SIMP密度刚度插值模型，通过惩罚因子惩罚后弹性模量的函数关系式为

(1)

式中：Ep为插值后的弹性模量；E0，Emin分别为实体和空洞部分的的弹性模型。

对不同的惩罚因子p，惩罚后材料的设计变量与弹性模量之间的关系如图1所示。

图1 SIMP模型设计变量与弹性模量关系

由材料力学知，在给定荷载和边界条件下使结构最优，则结构的柔度要最小化(或刚度最大化，应变能力最小化)，约束条件为模型的体积，建立拓扑优化的数学模型为

(2)

式中：[K]为结构的刚度矩阵；{U}为结构的位移矩阵；C(X)为结构的柔度矩阵;S.T.表示“使…满足”之意。

SIMP插值模型的刚度矩阵、柔度矩阵和敏度矩阵函数分别为

(3)

(4)

(5)

式中：Ki为第i个单元的“单元”刚度矩阵；ΔE=E0-Emin；C′(X)为结构柔度的敏度矩阵；Xi为单元的设计变量，为了刚度矩阵不产生奇异矩，所以取最小值为0.001。

2.2 RAMP插值模型

RAMP插值模型与SIMP插值模型的基本形式相同，只是优化后设计变量和弹性模量之间的关系式不同。其弹性模量函数为

(6)

对不同的惩罚因子q，惩罚后的材料弹性模量与设计变量之间的关系如图2所示。

图2 RAMP模型设计变量与弹性模量关系

RAMP插值模型的的刚度矩阵、柔度矩阵和敏度矩阵函数分别为

(7)

(8)

(9)

3 实例应用及成果分析

3.1 实例概况

某重力坝[8-9]上游水位95 m，下游水位5 m，坝体混凝土材料参数[10]见表1；上游面泥沙淤积高30 m，泥沙浮重度为8 kN/m3，内摩擦角为20°；坝基主排水孔扬压力折减系数0.3，距上游面7 m，水重度10 kN/m3；滑动面的抗剪摩擦因数0.7，坝体与基岩连接抗剪断摩擦因数1.45，基岩凝聚力1.084 MPa；根据《混凝土重力坝设计规范》(SL 349—2005)中规定按抗剪强度计算的允许抗滑稳定安全系数K=1.10，按抗剪断强度计算的允许抗滑稳定安全系数K′= 3.0。

表1 坝体混凝土材料参数

3.2 参数化建模及优化方法

(1) 首先根据基本资料建立参数化模型，拟定坝底宽度75 m，坝顶宽度7 m，距上游宽度7 m区域为非优化区，其余为优化区。坝基宽度分别在上游、下游取150 m，地基深度取150 m；模型施加荷载：水压力、泥沙压力、扬压力。地基底边和左右边为固定边界。模型单元网格尺寸为2 m×2 m。

(2) 之后建立优化任务，设定响应为应变能和体积，将优化目标设为应变能最小，体积作为约束条件，设定初始优化体积为50%。将非优化区域冻结，使得优化过程中区域内的材料不被删除。设定不同的惩罚因子，分别采用SIMP密度刚度插值模型和RAMP密度刚度插值模型的2种变密度法对模型进行优化。

3.3 优化结果与分析

根据SIMP插值模型的设计变量与弹性模量的关系(见图1)，设定对应的p值，得到不同的优化模型如图3所示，不同p值对应不同的迭代步见表2。

图3 SIMP插值优化结果

表2 SIMP插值法中不同p值的优化特征值

同理，根据RAMP插值模型的设计变量与弹性模量的关系(见图2)，设定对应的q值，得到不同的优化模型如图4所示，不同q值对应不同的迭代步如表3所示。

图4 RAMP插值优化结果

表3 RAMP插值法中不同p值的优化特征值

由图3可知，对SIMP插值优化后得到的优化模型都不规则，且在连续体中间出现大量空洞，完全达不到实体重力坝的剖面形状。由图4可知，对RAMP插值模型：当惩罚因子q较小时，优化的模型虽然与实体重力坝断面形状相似，但是坝体中也出现大量空洞，如果用该优化模型进行静力分析，需要进行人工填补模型空洞，将增加模型的后期处理；当q=10时，优化模型基本符合实体重力坝体断面形状，只需要对下游模型进行数据拟合即可；当惩罚因子q继续增大，优化模型无空洞且优化后模型和实体重力坝断面形状大体相同，但优化模型下游出现凸点，并且随惩罚因子的增大下游坝体形状变得不规则，因此无法准确拟定下游几何尺寸。

由表2可知，对于SIMP模型，除了惩罚因子p=1外，其他惩罚因子的优化迭代步都大于50步。目标应变能随着惩罚因子的增大而增大，当惩罚因子p= 5，6时，目标应变能数量级增大，是因为在有限元计算时出现刚度矩阵数值异常。

由表3可知，对于RAMP模型,惩罚因子从1≤q<20，20≤q<45，45≤q<90迭代步分阶段逐渐缩小。当1≤q≤20时目标应变能逐渐增大，当20≤q≤45时目标应变能在逐渐缩小，当q=90时目标应变能再次增大。由此可知，优化迭代步随着惩罚因子增加呈现周期性变化，而目标应变能随着惩罚因子先增大，当达到20时又逐渐缩小，但仍大于q<20时的目标应变能。为了使目标应变能较小，能充分利用结构材料，惩罚因子应越小，所以选择惩罚因子应小于20较合理。而根据图4可知，惩罚因子q≤5时图形出现空洞现象，因此惩罚因子的最佳选择范围为5

从以上分析可得：基于变密度法拓扑优化，对不同的插值模型，不同的惩罚因子所得出的优化模型差别很大。对于实体重力坝断面优化，RAMP插值模型比SIMP插值模型更符合优化结果，并且收敛速度更快，更节约计算资源，RAMP模型惩罚因子的最佳范围为5

3.4 优化模型静力分析

通过对第3.3节的分析，最符合的优化断面为RAMP密度刚度插值模型中惩罚因子q=10时的优化模型，由于优化后的模型下游断面线出现不规则现象，所以对最优模型下游断面进行线性拟合和非线性拟合，拟合方程分别为

(10)

对拟合后的坝体断面进行静力分析得如下压应力云图和静力分析结果分别见图5和表4。

图5 优化拟合静力分析压应力云图

由图5可知，坝体无论拟合为直线还是曲线压应力均为正，未出现拉应力，表示坝体处于安全状态。通过进一步观察压应力云图可看出，坝体拟合为曲线断面时受力均匀，而传统直线断面在下游坝面折线转折点处出现较明显的应力集中现象。

表4 优化模型静力分析结果

由表4可知，无论坝体拟合为直线还是曲线，坝体的抗剪强度抗滑稳定性和抗剪断强度抗滑稳定性均满足规范要求。坝体拟合为直线时断面面积比原始设计方法[8]减小18.209%，比模糊优化方法[8]减小0.37%，比进化策略方法[9]增大0.065%，证明拟合直线断面已达到了优化效果。而拟合为曲线断面的面积比拟合为直线断面的面积还减少了1.287%。通过以上分析可知，坝体无论拟合为直线还是曲线均符合规范要求，并达到了优化效果，且拟合曲线断面时较直线断面更优。

4 结 论

本文以大型通用商业有限元软件ABAQUS为平台，对变密度法的SIMP密度刚度插值法和RAMP密度刚度插值法2种插值模型进行了理论分析，并通过改变这2种模型的惩罚因子对实体重力坝进行优化，得出以下结论：

(1) SIMP插值模型优化的坝体出现大量空洞且形状不规则，并且迭代步均较长，优化过程中消耗资源较高；RAMP插值模型，在惩罚因子较小时，坝体出现空洞，惩罚因子较大时，下游出现凸起点，拟合时难以确定下游坝体尺寸。通过对惩罚因子与迭代步和目标应变能的分析，当惩罚因子为5

(2) 优化后的坝体各项应力指标和稳定性指标均符合规范；优化拟合的坝体断面面积比原始设计断面面积和其他优化方法[8-9]优化的坝体断面面积均有一定程度的优化。证明：基于RAMP插值模型的变密度法对实体重力坝的优化符合实际要求，具有工程意义。

(3) 下游断面为曲线时，缓解了下游折坡点处出现的应力集中现象。证明：将实体重力坝下游断面设计为曲线型更具有实际工程意义。在实际工程中，可根据实际地形环境和施工的难易程度进行安全性与经济性分析，选择出最安全且经济成本较低的坝型。

参考文献：

[1] MICHELL A G M.The Limits of Economy of Materials in Frame Structures[J].Philosophical Magazine, 1904，8(47)：589-597.

[2] BENDOSE M P.Optimization of Structural Topology，Shape and Material[M].Berlin：Springer-Verlag , 1995.

[3] BENDSOE M P，KIJUCHI N.Generating Optimal Topologies in Structural Design Using a Homogenization Method[J].Computer Methods in Applied Mechanics and Engineering，1988, 71(2): 197-224.

[4] BENDSOE M P.Optimal Shape Design As a Material Distribution Problem[J].Structural Optimization，1989, 1(4): 193-202.

[5] XIE Y M, STEVEN G P.A Simple Evolutionary Procedure for Structural Optimization[J].Computers & Structures，1993, 49(5): 885-896.

[6] SIGMUND O，PETERSSON J.Numerical Instabilities in Topology Optimization: A Survey on Procedures Dealing with Checkerboards，Mesh-dependencies and Local Minima[J].Structural Optimization，1988，16(1):68-75.

[7] 孙 蓓.连续体结构拓扑优化理论研究及其在水利工程中的应用[D].南京：河海大学，2005.(SUN Bei.Topology Optimum Design Theory of Continuum Structure and Its Application on Hydraulic Engineering[D].Nanjing：Hohai University，2005.(in Chinese))

[8] 蔡 新,王德信.模糊优化在坝工结构设计中的应用[J].河海大学学报(自然科学版),1995，23(6)：10-15.(CAI Xin，WANG De-xin.Application of Fuzzy Optimization to Dam Structural Design[J].Journal of Hohai University (Natural Science)，1995,23(6)：10-15.(in Chinese))

[9] 孙林松，王德信，许世刚.进化策略在重力坝优化设计中的应用[J].河海大学学报(自然科学版)，2000，28 (4)：104-106.(SUN Lin-song, WANG De-xin, XU Shi-gang.Application of Evolutionary Strategies to Gravity Dam Optimization[J].Journal of Hohai University (Natural Science)，2000，28(4)：104-106.(in Chinese))

[10]王金昌，陈页开.ABAQUS在土木工程中的应用[M].杭州：浙江大学出版社, 2006.(WANG Jin-chang，CHEN Ye-kai.Application of ABAQUS to Civil Engineering[M].Hangzhou: Zhejiang University Press,2006.(in Chinese))