基于BIM建模的坝体结构设计及快速有限元分析
2019-11-28张道法巩绪威
张道法，巩绪威
(上海勘测设计研究院有限公司，上海 200434)
数值计算辅助设计已成为水工结构设计过程中的重要环节，对于坝体结构设计而言，通常需要进行结构优化设计，对不同方案进行对比，选择一个最佳设计方案。坝体结构优化分析时需要进行大量结构计算，其中最主要内容之一就是有限元分析。水工结构通常形态各异、体形复杂，难以采用有限个参数表征其形状信息，即难以采用参数化建模的方式进行快速有限元分析[1]。不同方案进行比选时，计算工况众多，模型形状多样、无规律，有限元数值模型建立过程费时费力，数值分析效率低下。在生成算稿时又需要将大量计算云图整理成计算报告，工作量巨大。如何快速高效地进行有限元数值计算成为一个现实问题，郭瑞[2]、于彦伟[3]、蒋绘静[4]、王文进[5]等人对重力坝参数化设计、体形优化及参数化有限元的实现方法进行了大量研究；卢羽平[6]、孙斌[7]对水闸结构的参数化表征方法进行了系统分析，并采用参数化方式快速实现水闸的结构建模及有限元分析；陈艳[8]、伍鹤皋[9]实现了水电站钢岔管等结构的参数化设计，加快了岔管设计效率；钟登[10]深入研究了水电站地下厂房三维参数化方法；李松辉[11]对水工结构采用APDL和遗传算法优化设计进行了初步分析。
以上研究主要是采用参数化设计方式来加快设计或计算效率，对于结构规则构件相对易于实现，而对于复杂水工结构的快速有限元分析的研究相对较少。随着BIM技术在水利工程领域的快速应用，BIM三维模型的创建过程日趋方便[12]。从挖掘BIM模型应用潜力的角度出发，本文在研究三维BIM模型与有限元数值计算耦合应用时，总结了一套基于BIM模型和ABAQUS计算平台快速进行有限元分析的一般流程，实现了BIM模型与有限元计算之间的信息共享，有效提升了数值模型创建、数值分析及计算报告整理的工作效率。
1 快速有限元分析流程
a) 采用BIM软件创建三维实体模型，导入ABAQUS软件创建三维数值模型。利用BIM软件创建三维模型是十分方便的，其建模效率比传统有限元软件建模高出数倍。本文采用Microstation快速建立三维几何模型，并将其存储为sat中间格式，导入ABAQUS软件中进行数值模型的转换[13]。划分有限元网格是建立有限元数值模型、进行有限元计算工作量最大的部分，而水工结构常常由于体型各异、形态复杂很难使用参数化有限元方式进行处理，目前还是以手工划分网格单元为主，费时费力。本文基于这个问题，采用10节点二次单元代替传统8节点单元进行网格剖分处理，因10节点四面体二次单元网格划分可采用自由剖分方法，无需大量切割实体，对实体形状要求低，单元生成速度快。BIM模型仅需稍加处理，添加地基结构即可快速剖分网格，生成有限元模型。这种方式能够大量节省网格剖分时间，尤其在多方案比选时，具有良好效果。当然，也应该看出采用四面体二次单元会导致网格畸形概率增大、单元数量增多、计算时间加长的现象[14]。但通过合理设置网格密度，网格单元尽量分区设置，在结构复杂区域单元尺寸设置相对小一些，结构简单区域单元尺寸设置相对大一些，远离结构地基单元尺寸可由小尺寸逐渐过渡到较大尺寸。同时应适当修改BIM模型，如结构细小尖角等易引起应力集中部位可作适当简化处理，网格畸形问题可以较好地解决，而现阶段随着计算机性能的不断提高，计算时间长的问题也越来越不是限制性问题。
b) 采用inp文件、荷载子程序for文件及py命令流快速进行边界条件施加、截图及生成图表等后处理操作。网格剖分后，需要施加荷载及边界条件，水工结构有限元计算时，此设置过程基本一致。而针对不规则分布的扬压力、土压力等特殊荷载设置过程则较为复杂，公式难以直接在界面中施加在结点或单元上，且各工况难以重复利用。采用ABAQUS命令流与荷载子程序能够较好地解决公式编写困难、荷载施加繁琐等问题，对于不同工程仅需要简单复制代码并修改即可。材料参数、位移边界条件、初始条件等方便修改的参数在inp文件中设置，复杂荷载建议在荷载子程序中设置，荷载子程序在设置复杂荷载时具有较大的优势，公式在程序中的编写与实际公式基本一致，采用参数化编制的公式针对不同工况仅需要修改参数值即可，能够显著提高施加复杂荷载的效率。ABAQUS一般荷载施加程序接口如下[15]。
SUBROUTINE LOAD(F,KSTEP,KINC,TIME,NOEL,NPT,LAYER, &
KSPT,COORDS, JLTYP,SNAME)
INCLUDE 'ABA_PARAM.INC'
DIMENSION TIME(2), COORDS(3)
CHARACTER*80 SNAME
REAL X,Y,Z
X=COORDS(1)
Y=COORDS(2)
Z=COORDS(3)
! 编写荷载公式
RETURN
END
针对有限元后处理复杂繁琐等问题，本文采用脚本批量处理模式，录制一个工况下的计算云图及报表的截取保存过程，生成批量化py脚本文件，其他工况直接载入该脚本文件，可迅速生成命名有规律的云图及图表，操作过程详见参考文献[16]。
c) 利用宏语言将计算云图批量导入报告并批量进行格式调整。有限元计算完毕后会产生大量的计算云图，手工将这些图片插入word报告中需要耗费大量的时间精力，而且操作过程简单重复、枯燥乏味、容易出错。利用宏可以方便地将计算云图及报表批量导入word计算报告中，且能对导入的图片自动进行居中对齐、大小调整等格式化处理，操作简单、方便、高效。主要操作代码如下。
Sub ImportPicturesInBatch()
FilePath = “D:
PicNum = 2′每个工况图片张数
CaseNum = 7′工况数量
StartIndex=
ActiveDocument.InlineShapes.Count + 1
For kk = 1 To PicNum ′导入图片
For cc = 1 To CaseNum
Selection.InlineShapes.AddPictureFilePath& cc &"_"& kk &".png"
Selection.TypeParagraph
Selection.TypeParagraph
Selection.ParagraphFormat.Alignment=
wdAlignParagraphCenter
Next cc
Next kk
EndIndex=
ActiveDocument.InlineShapes.Count
For i = StartIndex To EndIndex ′设置图片大小及格式
With ActiveDocument.InlineShapes(i)
.LockAspectRatio = msoTrue ′锁定图片纵横比
.Width = 250
End With
ActiveDocument.InlineShapes(i).Select
Selection.ParagraphFormat.Alignment=
wdAlignParagraphCenter
Next i
End Sub
2 工程案例
2.1 工程概况
国外某水电站位于严寒地区，采用河床式布置，厂房与坝体共同挡水，分为15个坝段，总长度为291.25 m，从左岸向右岸依次为：左岸连接坝段(1—3号)、左岸扶壁挡墙砂砾石混合坝段(4—9号)、安装间上游挡水重力坝坝段(10号)、厂房坝段(11号)、泄洪深孔坝段(12号)、溢流表孔坝段(13号)、右岸重力坝段(14—15号)。左岸扶壁挡墙砂砾石混合坝段最大坝高30.50 m，坝顶宽度4.50 m。扶壁式挡墙底板厚度2.00～3.00 m，墙身面板厚度2.50 m，扶壁厚1.00～2.00 m，扶壁中心距5.00～7.00 m。挡墙墙身上、下游侧采用卵砾石回填，上游侧为顺岸坡向回填，回填高程571.00～546.27 m，坡比为1.0∶4.5；下游侧回填顶高程568.00 m，垂直坝轴线方向坡比1∶2，护坡采用干砌块石，厚度400 mm。
本文坝体结构设计为左岸9号扶壁挡墙坝段，扶壁挡墙底板顺水流方向长24.40 m，垂直水流方向宽18.00 m，挡墙顶高程为571.00 m，建基面高程为540.50～543.00 m，挡墙顶宽4.50 m。坝体墙身及底板采用C30W6F200混凝土，扶壁采用C25W4F200混凝土。右侧岸扶壁与中间扶壁之间设置横向支撑挡墙，由此形成的中间框格内回填砂砾石，梁扶壁顶设置楼梯通向坝顶，右侧扶壁及支撑挡墙下游坡比1.0∶0.4，以便与10号坝段平顺连接，外形协调美观。中间扶壁坡比1∶2，下游另设挡墙以便侧向支挡坝体填料。左侧扶壁坡比1.000∶0.735，9号坝段结构相对复杂，其典型纵断面见图1。
2.2 坝体比选方案
对9号坝段右岸扶壁与中间扶壁之间的横向支撑进行优化比选设计。考虑右岸扶壁与10号重力坝段衔接的功能需要，坝后回填料仅填至中间扶壁处，中间扶壁一侧挡土，一侧临空，由此填料对中间扶壁产生了较大的侧向土压力，鉴于中间扶壁结构的安全性，需要在中间扶壁与右岸扶壁之间设置横向支撑， 共拟定2种支撑方案进行比选。方案一为小挡墙支撑方案，小挡墙顶高程为556.00 m，顶部宽0.8 m，以1.0∶0.4坡与坝底板衔接，横断面见图2。方案二为全断面支撑墙方案，支撑墙厚度为1 m，形状与右岸扶壁一致，以便与10号安装间重力坝段衔接，横断面见图3。由图2—3可知，2个方案形状差异较大，难以采用参数化方式提高效率，优化比选时建模、计算、后处理等过程十分繁琐。
图1 9号坝段典型纵断面(mm)
图3 方案二全断面支撑墙方案(mm)
2.3 计算条件
采用上述建立的有限元快速计算流程，快速对2个支撑方案建立了三维模型，进行结构网格剖分、荷载及边界条件施加及有限元分析，坝体结构有限元计算模型见图4。
计算时以正常蓄水工况为例进行分析，正常蓄水位工况下上游水位为567.50 m，下游无水，坝体及地基材料均采用线弹性本构，采用的材料参数见表1。
因坝后回填料仅在中间扶壁与左岸扶壁间设置，即中间扶壁右岸侧为临空侧，左岸侧为临土侧，中间扶壁承受较大的不均衡土压力。
a)方案一
b)方案二
表1 计算参数
扶壁切向倾斜面采用库伦土压力理论，作用在某一点处的土压力为：
pA=γhKA
(1)
式中pA——某一点处的土压力强度，kN/m；γ——土的重度，对于地下水位以下部分取浮容重γ'；h——填土高；KA——库伦主动土压力系数，计算过程见下式。
KA=
(2)
式中φ——土的内摩擦角；θ——墙背与竖直向夹角；δ——填土与墙背摩擦角，取δ=1/3φ；β——填土表面坡角。
其余直立面采用静止土压力计算，静止土压力系数为：
K0=1-sinφ
(3)
由以上公式可以看出，土压力公式较为复杂，土压力作用面较多且倾斜面与直立面采用不同土压力计算方法，界面施加过程十分复杂。本文采用荷载子程序方式进行施加，根据上述土压力计算公式编写荷载子程序，可以方便地为各土压力作用面施加土压力，方便高效，且不同工况仅需修改控制参数，代码可重复使用，效率较高。
2.4 计算结果
以正常蓄水位工况为例，2个方案的计算结果见图5—7。方案一的最大拉应力为3.426 MPa，方案二的最大拉应力为3.298 MPa，均发生在中间扶壁左岸侧下部位置，两方案的最大拉应力已超出混凝土抗拉强度，需要配置抗拉钢筋。方案一的最大压应力-10.230 MPa，方案二的最大压应力为-8.894 MPa，均发生在中间扶壁右岸侧下部位置，未超过混凝土抗压强度。方案一的最大水平侧向位移为7.795 mm，方案二的最大水平侧向位移为7.668 mm，均发生在中间扶壁顶部，指向临空侧。2种支撑形式的方案受力与变形规律基本一致，但全断面支撑墙方案的最大拉应力比方案一小0.128 MPa，最大压应力小1.336 MPa，最大水平侧向位移小0.127 mm。因而，方案二相较方案一在结构安全上更优，全断面支撑墙方案能够给中间扶壁提供更好的支撑作用，两方案在工程量基本相当的情况下，显然方案二的受力与变形更小一些。同时，内部隔腔内回填砂砾石对上游挡水面墙有很好的保温作用，因此设计时最终选择了方案二。
a)方案一
b)方案二
图5最大拉应力计算结果对比
a)方案一
b)方案二
图6最大压应力计算结果对比
a)方案一
b)方案二
图7最大水平侧向位移计算结果对比
3 结论
a) 本文在研究BIM模型与有限元计算耦合应用的基础上，总结了快速进行有限元分析的一般流程，使用BIM软件建立三维模型，利用ABAQUS软件提供的计算及批量化处理功能进行有限元分析，采用自编程序进行计算报告整理，该有限元分析流程可快速进行数值模型建立、有限元分析及计算报告生成，有效提升了复杂结构数值分析的效率。
b) 将建立的快速有限元分析方法应用于某坝体结构支撑挡墙优化设计过程中，针对该坝体结构的多种比选方案，使用本文的有限元快速分析流程准确建立了多种方案的数值计算模型，快速进行坝体复杂荷载的施加，并进行多方案的比选分析。实际应用结果表明，本文的快速有限元分析法可显著提升不规则坝体结构优化设计的效率。
c) 对坝体结构的优化分析表明，设计中采用的全断面支撑墙方案相比小挡墙支撑方案能够将中间扶壁的最大拉应力降低0.128 MPa，最大压应力降低1.336 MPa，最大侧向位移降低0.127 mm，全断面支撑墙方案中间的受力与变形相对更小一些，能够给予中间扶壁更好的支撑作用。因此，设计选择了全断面支撑墙方案。分析结果论证了结构设计的合理性，为坝体设计方案的选择提供了依据。
