大坝结构静动力分析的精细化模型
林皋 庞林
摘要：阐述近期大连理工大学抗震研究所在大坝抗震分析方面所取得的研究成果，为了提高大坝结构静动力分析的计算精度与效率，提出多边形单元的计算技术与奇异应力场计算方法，该方法使应力计算成果较有限元法有很大的改进，且计算工作量比较节省。
关键词：静动力分析；多边形单元；奇异应力场
中图分类号：TU45 文献标识码：A 文章编号：1000-0666（2016）01-0001-09
0 前言
大坝结构的安全对城市抗震防灾有重要影响。2008年5月12日发生的汶川8.0级特大地震对位于成都市西北60km岷江上游的紫坪铺大坝造成了一定的地震损伤，威胁到成都市的安全。吉林市的安全与丰满大坝的安全息息相关，密云水库、官厅水库的安全也对北京市的安全产生一定影响。研究城市防灾需要关注周边大坝的安全设防。
新中国成立以来，特别是改革开放以来，我国水利水电建设迅速发展。进入21世纪，世界大坝建设中心已经转向中国。一大批世界级的高坝已经或正在中国大地上兴建，坝高屡破世界纪录，其成就举世瞩目。我国地处世界两大活跃地震带——环太平洋地震带和欧亚地震带的交汇地段，地震活动频繁。为适应地震危险区大坝建设的需要，我国的大坝抗震技术也蓬勃发展，并进入世界先进行列。
自20世纪50年代以来，大连理工大学工程抗震研究所一直根据我国大坝抗震的需要，致力于研究如何提高大坝抗震分析的水平。大坝抗震的特点在于大坝所处环境复杂（图1），参照世界坝工专家Chopra（1992，2008）的意见，大坝抗震的关键技术包括坝与库水的动力相互作用分析，坝与无限地基的动力相互作用分析，筑坝材料的动态特性及大坝地震损伤发展与风险分析。在2013年中国工程院举办的“中国工程科技论坛——水安全与水利水电可持续发展”会议上，林皋（2014）提出了大坝抗震模型的改进意见。关于坝与库水的动力相互作用分析需要考虑水的可压缩性以及动水压力波在水库边界的吸收作用，Chopra（2008）采用有限元法（FEM）进行求解，动水压力按坝的振动模态展开，计算复杂、难以实际应用。Lin等（2007a，2012a）采用比例边界有限元法（SBFEM）进行求解，按库水的模态进行展开，计算简便、精度高、便于实际应用。关于坝与无限地基的动力相互作用，由于计算模型比较复杂，文献中一般都将无限地基简化为均质地基进行考虑，难以反映实际情况。我们将无限地基分为层状分布、分块不均质和近场杂乱不规则等几种类型，分别提出了不同的计算模型（Linet al，2015，2007b；Yin et al，2013），计算结果都具有较好的精度与效率。为了提高坝-库水-地基系统时域分析的计算精度与效率，Lin等（2012b）提出以等几何分析方法（isogeometric a-nalysis，简称IGA）代替FEM来进行坝的应力分析，并将IGA与SBFEM相结合来进行坝与库水和坝与无限地基的动力相互作用分析。关于坝材料的动力特性方面，混凝土是速率敏感性材料。自从第一颗原子弹1945年在日本广岛爆炸以后，为了防卫核爆炸造成的灾难，许多国家都加强了混凝土材料抗击核冲击的研究。但这种研究偏重单调瞬时冲击，而且抗压性能居多，抗拉性能偏少。地震荷载的特点是不规则的循环往复，而且加载速率也偏低。为此，我们开展了2000多试件混凝土在地震作用下的动态特性的研究，取得了一批有价值的成果（Yan，Lin，2006；Lin et al，2007c）。
笔者认为除了计算模型以外，数值计算方法也对大坝抗震分析的计算精度与效率产生影响。大坝地震响应分析中，有限元方法的应用比较普遍。有限元法可以适应各种不规则的几何边界条件和分布不规律的单元材料特性，可以应用于各种类型物理问题的求解，计算能力强，不受问题求解自由度的限制。但有限元法也有其不足之处，其主要采用低阶插值函数，连续性差，应力计算精度低。但在工程设计中应力是对大坝安全性评价的重要依据，所以，有必要寻求更为有效的其它数值计算方法。根据我们研究，提出以下计算模型和方法供参考，欢迎批评指正。
1 多边形单元技术
多边形单元（PE-Polygon Element）不是一个新概念，但经过澳大利亚新南威尔士大学（UN-SW）宋崇民教授团队的发展，建立在SBFEM基础上的多边形单元已经成为与FEM一样具有广泛适应性的网格离散技术，对数值计算具有强大的功能。PE可称为广义的FE。
大连理工大学发展了PE技术在结构静动力分析中的应用，主要工作为：（1）单元细化与连续性升阶可以只在PE的边界上进行，不改变网格的拓朴结构，使应力计算精度显著提高，但计算工作量则很节约；（2）将PE与曲线边界相结合，既可构造纯直线边的PE，也可构造曲线边的PE，或直线边与曲线边相结合的PE；（3）将PE应用于结构的奇异应力场分析，可以达到很高的计算精度。
多边形单元按SBFEM进行建模（图2）。单元域内选择一点0作为比例相似中心，其要求是从0点可以看见边界上各点，或0点可与边界上任意点用直线连接，没有障碍。单元各边上可以任意划分子单元和结点，数目不限。如图2所示，PE为六边形，每边可以只划分一个子单元，每个子单元含3个结点，则具有2阶精度。AB边也可以划分2个子单元，每个子单元各3个结点共5个结点，这就相当于网格加密。各子单元用SBFEM进行建模，由于SBFEM为半解析方法，从0点至边界AM上任意点相连的射线上，变量具有解析解，这就使PE计算域具有很高的连续性，使应力计算达到很高的计算精度。
含有曲线边的多边形单元如图3所示。这时，含曲线边AF的子单元1e按等几何分析（IGA）建模，采用非均匀有理B样条（NURBS）作为插值基函数，即使在很粗网格的条件下也可以准确地描述几何外形。IGA充分利用了NURBS的优点，可以在保持几何外形不变的条件下方便地进行网格细分，提升插值的阶次，提高计算的精度，所以很适于曲线形结构的建模。NURBS基函数由控制点进行描述，不具有插值性。图3中BE、CD、DE各边相应的子单元由SBFEM进行建模，AB、EF边相对应的两个子单元2e和6e按SBFEM建模，但是由于和1#相连需作一定调整。亦即A点的位移需用控制点1、2、3的位移加以表示，F点的位移需用控制点2、3、4的位移加以表示，这样就保持了PE单元位移场的连续性。
PE由域内三角形单元形心点连线加周边三角形单元中点连线构成（图4），所以建模非常方便。单元便于细分和拼合，如图5所示。
多边形单元具有以下特点：（1）适于各种几何形体的建模，具有高度灵活性；（2）单元具备单位分解性与线性完备性，通过分片检验；
（3）可以将网格加密进行单元细化，也可在不改变网格拓扑结构的条件下，在单元各边上增加结点实现单元细化和连续性升阶；（4）可以方便地进行粗细网格过度，只需在不同网格交界面上增加结点即可（图6）；（5）可以和FE和NuRBS兼容。
2 数值算例
为了检验PE的效果，通过以下算例加以说明。
2.1 厚壁圆筒承受均匀内压作用
本例有解析解。利用对称性，只取厚壁圆筒的1/4进行离散。假设内压力p=1，径向与环向应力只和内径和外径的比值有关，假设r/R=1/3。将PE与FE（二阶单元）的计算结果进行对比，两种方法的网格划分如图7所示。两种方法计算精度与解析解相比的误差随网格结点数的变化的对比如图8所示，图中PE代表直线边界网格，CPE代表内外曲线边界网格，由图可见PE的计算优越性是明显的。
2.2 曲线形溢流坝受静水压力作用
某工程溢流坝，高H=66m，受静水压力作用。坝材料参数假设为：弹性模量E=36.5GPa；泊松比γ=0.167。采用PE与FE两种方法对该工程溢流坝进行计算对比。坝体两种方法网格划分如图9所示，廊道部位网格划分如图10所示。图10图题中p表示单元离散所用的插值阶数。廊道上部曲线边界的应力分量对比如图11所示，廊道周边第一主应力分布如图12所示。由图可见，有限元网格即使在廊道角点O附近细分，应力奇异性仍无法反映。
3 异应力场分析
有限元方法无法处理应力奇异性问题。SBFEM由于径向具有解析解，所以特别有利于处理应力奇异性问题。关键是须将应力奇异点设置于相似中心处。但是SBFEM中，通过相似中心的边界不能进行离散，只能将此边界上作用的外力或变形转化为相应的边界条件，在计算方程中加以体现。我们通过算例来加以说明。
研究坝踵点在水压力作用下的奇异应力分布。图13表示PE和FE两种网格划分。
图14表示PE与FE在坝踵点附近的网格形状和节点分布。图15表示各种方法的计算结果。图中SBFEM相似中心位于域内，但采用不同的结点数进行离散。PE为相似中心位于坝踵A点，但采用不同的边界条件。由图可见，FEM即使进行网格细分都无法表示奇异应力场的分布，采用SBFEM时，如不将相似中心置于应力奇异点处也难以描述坝踵点附近的奇异应力分布。采用PE方法同时将相似中心置于坝踵点有可能获得接近奇异应力场分布的结果。这时需通过其它方法求得AB和AC两边界面上的位移，以建立相应的AB和AC边的Dirichlet边界条件。但处理方式不同所取得的效果也有差别。笔者采用了两种方式处理。PE1中AB和AC两边界均按其它方法计算出的位移建立Dirichlet边界条件进行计算，而PE2中边界AC采用已知水压力建立Neumann边界条件，AB边界仍采用其它方法计算的位移建立Dirichlet边界条件。可以看出PE2的效果比较好，这是因为AC边的边界条件比较准确，同时AC边上的水压力也对奇异应力场影响比较大的缘故。因此，为了取得良好的奇异应力场的效果，建立边界面的准确边界条件仍然很重要。
动力问题奇异应力场的分析也可采用相同的方法进行处理。图16～20为地震波作用下坝踵点奇异应力场的分析结果。采用的材料参数列入表1，为了表示PE的优越性，将PE结果与比较准确的IGA的计算结果进行比较。IGA的各级网格均可自动加密，网格采用2阶单元。图18表示奇异应力场分析结果的比较。由图可见，在距坝踵点水平距离x=1.0m以外，各种网格的计算结果都收敛而且相符性比较好。但距离坝踵点x=0.5m范围以内，IGA1-IGA3偏差较大，IGA4结果与本文方法吻合良好，达到收敛要求。图19为坝踵点的位移与加速度响应时程比较。由于位移不存奇性，与应力结果不同，各网格计算结果相差不大。
4 结语
按比例边界有限元法（SBFEM）原理建立的多边形单元（PE）是一种提高应力计算精度的高效数值方法。多边形单元可以进一步发展为多面体单元，即将二维问题发展到三维问题，具有广阔的发展前景。
摘要：阐述近期大连理工大学抗震研究所在大坝抗震分析方面所取得的研究成果，为了提高大坝结构静动力分析的计算精度与效率，提出多边形单元的计算技术与奇异应力场计算方法，该方法使应力计算成果较有限元法有很大的改进，且计算工作量比较节省。
关键词：静动力分析；多边形单元；奇异应力场
中图分类号：TU45 文献标识码：A 文章编号：1000-0666（2016）01-0001-09
0 前言
大坝结构的安全对城市抗震防灾有重要影响。2008年5月12日发生的汶川8.0级特大地震对位于成都市西北60km岷江上游的紫坪铺大坝造成了一定的地震损伤，威胁到成都市的安全。吉林市的安全与丰满大坝的安全息息相关，密云水库、官厅水库的安全也对北京市的安全产生一定影响。研究城市防灾需要关注周边大坝的安全设防。
新中国成立以来，特别是改革开放以来，我国水利水电建设迅速发展。进入21世纪，世界大坝建设中心已经转向中国。一大批世界级的高坝已经或正在中国大地上兴建，坝高屡破世界纪录，其成就举世瞩目。我国地处世界两大活跃地震带——环太平洋地震带和欧亚地震带的交汇地段，地震活动频繁。为适应地震危险区大坝建设的需要，我国的大坝抗震技术也蓬勃发展，并进入世界先进行列。
自20世纪50年代以来，大连理工大学工程抗震研究所一直根据我国大坝抗震的需要，致力于研究如何提高大坝抗震分析的水平。大坝抗震的特点在于大坝所处环境复杂（图1），参照世界坝工专家Chopra（1992，2008）的意见，大坝抗震的关键技术包括坝与库水的动力相互作用分析，坝与无限地基的动力相互作用分析，筑坝材料的动态特性及大坝地震损伤发展与风险分析。在2013年中国工程院举办的“中国工程科技论坛——水安全与水利水电可持续发展”会议上，林皋（2014）提出了大坝抗震模型的改进意见。关于坝与库水的动力相互作用分析需要考虑水的可压缩性以及动水压力波在水库边界的吸收作用，Chopra（2008）采用有限元法（FEM）进行求解，动水压力按坝的振动模态展开，计算复杂、难以实际应用。Lin等（2007a，2012a）采用比例边界有限元法（SBFEM）进行求解，按库水的模态进行展开，计算简便、精度高、便于实际应用。关于坝与无限地基的动力相互作用，由于计算模型比较复杂，文献中一般都将无限地基简化为均质地基进行考虑，难以反映实际情况。我们将无限地基分为层状分布、分块不均质和近场杂乱不规则等几种类型，分别提出了不同的计算模型（Linet al，2015，2007b；Yin et al，2013），计算结果都具有较好的精度与效率。为了提高坝-库水-地基系统时域分析的计算精度与效率，Lin等（2012b）提出以等几何分析方法（isogeometric a-nalysis，简称IGA）代替FEM来进行坝的应力分析，并将IGA与SBFEM相结合来进行坝与库水和坝与无限地基的动力相互作用分析。关于坝材料的动力特性方面，混凝土是速率敏感性材料。自从第一颗原子弹1945年在日本广岛爆炸以后，为了防卫核爆炸造成的灾难，许多国家都加强了混凝土材料抗击核冲击的研究。但这种研究偏重单调瞬时冲击，而且抗压性能居多，抗拉性能偏少。地震荷载的特点是不规则的循环往复，而且加载速率也偏低。为此，我们开展了2000多试件混凝土在地震作用下的动态特性的研究，取得了一批有价值的成果（Yan，Lin，2006；Lin et al，2007c）。
笔者认为除了计算模型以外，数值计算方法也对大坝抗震分析的计算精度与效率产生影响。大坝地震响应分析中，有限元方法的应用比较普遍。有限元法可以适应各种不规则的几何边界条件和分布不规律的单元材料特性，可以应用于各种类型物理问题的求解，计算能力强，不受问题求解自由度的限制。但有限元法也有其不足之处，其主要采用低阶插值函数，连续性差，应力计算精度低。但在工程设计中应力是对大坝安全性评价的重要依据，所以，有必要寻求更为有效的其它数值计算方法。根据我们研究，提出以下计算模型和方法供参考，欢迎批评指正。
1 多边形单元技术
多边形单元（PE-Polygon Element）不是一个新概念，但经过澳大利亚新南威尔士大学（UN-SW）宋崇民教授团队的发展，建立在SBFEM基础上的多边形单元已经成为与FEM一样具有广泛适应性的网格离散技术，对数值计算具有强大的功能。PE可称为广义的FE。
大连理工大学发展了PE技术在结构静动力分析中的应用，主要工作为：（1）单元细化与连续性升阶可以只在PE的边界上进行，不改变网格的拓朴结构，使应力计算精度显著提高，但计算工作量则很节约；（2）将PE与曲线边界相结合，既可构造纯直线边的PE，也可构造曲线边的PE，或直线边与曲线边相结合的PE；（3）将PE应用于结构的奇异应力场分析，可以达到很高的计算精度。
多边形单元按SBFEM进行建模（图2）。单元域内选择一点0作为比例相似中心，其要求是从0点可以看见边界上各点，或0点可与边界上任意点用直线连接，没有障碍。单元各边上可以任意划分子单元和结点，数目不限。如图2所示，PE为六边形，每边可以只划分一个子单元，每个子单元含3个结点，则具有2阶精度。AB边也可以划分2个子单元，每个子单元各3个结点共5个结点，这就相当于网格加密。各子单元用SBFEM进行建模，由于SBFEM为半解析方法，从0点至边界AM上任意点相连的射线上，变量具有解析解，这就使PE计算域具有很高的连续性，使应力计算达到很高的计算精度。
含有曲线边的多边形单元如图3所示。这时，含曲线边AF的子单元1e按等几何分析（IGA）建模，采用非均匀有理B样条（NURBS）作为插值基函数，即使在很粗网格的条件下也可以准确地描述几何外形。IGA充分利用了NURBS的优点，可以在保持几何外形不变的条件下方便地进行网格细分，提升插值的阶次，提高计算的精度，所以很适于曲线形结构的建模。NURBS基函数由控制点进行描述，不具有插值性。图3中BE、CD、DE各边相应的子单元由SBFEM进行建模，AB、EF边相对应的两个子单元2e和6e按SBFEM建模，但是由于和1#相连需作一定调整。亦即A点的位移需用控制点1、2、3的位移加以表示，F点的位移需用控制点2、3、4的位移加以表示，这样就保持了PE单元位移场的连续性。
PE由域内三角形单元形心点连线加周边三角形单元中点连线构成（图4），所以建模非常方便。单元便于细分和拼合，如图5所示。
多边形单元具有以下特点：（1）适于各种几何形体的建模，具有高度灵活性；（2）单元具备单位分解性与线性完备性，通过分片检验；
（3）可以将网格加密进行单元细化，也可在不改变网格拓扑结构的条件下，在单元各边上增加结点实现单元细化和连续性升阶；（4）可以方便地进行粗细网格过度，只需在不同网格交界面上增加结点即可（图6）；（5）可以和FE和NuRBS兼容。
2 数值算例
为了检验PE的效果，通过以下算例加以说明。
2.1 厚壁圆筒承受均匀内压作用
本例有解析解。利用对称性，只取厚壁圆筒的1/4进行离散。假设内压力p=1，径向与环向应力只和内径和外径的比值有关，假设r/R=1/3。将PE与FE（二阶单元）的计算结果进行对比，两种方法的网格划分如图7所示。两种方法计算精度与解析解相比的误差随网格结点数的变化的对比如图8所示，图中PE代表直线边界网格，CPE代表内外曲线边界网格，由图可见PE的计算优越性是明显的。
2.2 曲线形溢流坝受静水压力作用
某工程溢流坝，高H=66m，受静水压力作用。坝材料参数假设为：弹性模量E=36.5GPa；泊松比γ=0.167。采用PE与FE两种方法对该工程溢流坝进行计算对比。坝体两种方法网格划分如图9所示，廊道部位网格划分如图10所示。图10图题中p表示单元离散所用的插值阶数。廊道上部曲线边界的应力分量对比如图11所示，廊道周边第一主应力分布如图12所示。由图可见，有限元网格即使在廊道角点O附近细分，应力奇异性仍无法反映。
3 异应力场分析
有限元方法无法处理应力奇异性问题。SBFEM由于径向具有解析解，所以特别有利于处理应力奇异性问题。关键是须将应力奇异点设置于相似中心处。但是SBFEM中，通过相似中心的边界不能进行离散，只能将此边界上作用的外力或变形转化为相应的边界条件，在计算方程中加以体现。我们通过算例来加以说明。
研究坝踵点在水压力作用下的奇异应力分布。图13表示PE和FE两种网格划分。
图14表示PE与FE在坝踵点附近的网格形状和节点分布。图15表示各种方法的计算结果。图中SBFEM相似中心位于域内，但采用不同的结点数进行离散。PE为相似中心位于坝踵A点，但采用不同的边界条件。由图可见，FEM即使进行网格细分都无法表示奇异应力场的分布，采用SBFEM时，如不将相似中心置于应力奇异点处也难以描述坝踵点附近的奇异应力分布。采用PE方法同时将相似中心置于坝踵点有可能获得接近奇异应力场分布的结果。这时需通过其它方法求得AB和AC两边界面上的位移，以建立相应的AB和AC边的Dirichlet边界条件。但处理方式不同所取得的效果也有差别。笔者采用了两种方式处理。PE1中AB和AC两边界均按其它方法计算出的位移建立Dirichlet边界条件进行计算，而PE2中边界AC采用已知水压力建立Neumann边界条件，AB边界仍采用其它方法计算的位移建立Dirichlet边界条件。可以看出PE2的效果比较好，这是因为AC边的边界条件比较准确，同时AC边上的水压力也对奇异应力场影响比较大的缘故。因此，为了取得良好的奇异应力场的效果，建立边界面的准确边界条件仍然很重要。
动力问题奇异应力场的分析也可采用相同的方法进行处理。图16～20为地震波作用下坝踵点奇异应力场的分析结果。采用的材料参数列入表1，为了表示PE的优越性，将PE结果与比较准确的IGA的计算结果进行比较。IGA的各级网格均可自动加密，网格采用2阶单元。图18表示奇异应力场分析结果的比较。由图可见，在距坝踵点水平距离x=1.0m以外，各种网格的计算结果都收敛而且相符性比较好。但距离坝踵点x=0.5m范围以内，IGA1-IGA3偏差较大，IGA4结果与本文方法吻合良好，达到收敛要求。图19为坝踵点的位移与加速度响应时程比较。由于位移不存奇性，与应力结果不同，各网格计算结果相差不大。
4 结语
按比例边界有限元法（SBFEM）原理建立的多边形单元（PE）是一种提高应力计算精度的高效数值方法。多边形单元可以进一步发展为多面体单元，即将二维问题发展到三维问题，具有广阔的发展前景。
