基金项目:国家自然科学重点基金资助项目(51138001)、国家自然科学基金委托创新研究群体项目(51421064)和中央高校基本科研业务费专项基金(DUT15TD17)联合资助.
(1.大连理工大学 海岸与近海工程国家重点实验室,辽宁 大连 116024; 2.大连理工大学 工程抗震研究所建设工程学部水利工程学院,辽宁 大连 116024)
(1.State key Laboratory of Coastal and Offshore Engineering,Dalian Vniversityof Technology,Dalian 116024,Lianning,China)(2.Institute of Earthquake Engineering,Faculty of Infrastructure Engineering,Dalian University of Technology,Dalian 116024,Liaoning,China)
static and dynamic analysis; polygon element; singular stress field
备注
基金项目:国家自然科学重点基金资助项目(51138001)、国家自然科学基金委托创新研究群体项目(51421064)和中央高校基本科研业务费专项基金(DUT15TD17)联合资助.
阐述近期大连理工大学抗震研究所在大坝抗震分析方面所取得的研究成果,为了提高大坝结构静动力分析的计算精度与效率,提出多边形单元的计算技术与奇异应力场计算方法,该方法使应力计算成果较有限元法有很大的改进,且计算工作量比较节省。
The recent progress of the research which is conducted by the Faculty of Infrastructure Engineering,Dalian University of Technology is briefly introduced. To enhance the accuracy and efficiency of the static and dynamic stress analysis for dam structures,the calculation technique for exploiting salient functions of polygon elements and the method for evaluating singular stress field has been developed. Compared with the finite element analysis,these achievements results of stress analysis have been greatly improved. Meanwhile,the computational effort is reduced to a great extent.
引言
大坝结构的安全对城市抗震防灾有重要影响。2008年5月12日发生的汶川8.0级特大地震对位于成都市西北60 km岷江上游的紫坪铺大坝造成了一定的地震损伤,威胁到成都市的安全。吉林市的安全与丰满大坝的安全息息相关,密云水库、官厅水库的安全也对北京市的安全产生一定影响。研究城市防灾需要关注周边大坝的安全设防。
新中国成立以来,特别是改革开放以来,我国水利水电建设迅速发展。进入21世纪,世界大坝建设中心已经转向中国。一大批世界级的高坝已经或正在中国大地上兴建,坝高屡破世界纪录,其成就举世瞩目。我国地处世界两大活跃地震带——环太平洋地震带和欧亚地震带的交汇地段,地震活动频繁。为适应地震危险区大坝建设的需要,我国的大坝抗震技术也蓬勃发展,并进入世界先进行列。
自20世纪50年代以来,大连理工大学工程抗震研究所一直根据我国大坝抗震的需要,致力于研究如何提高大坝抗震分析的水平。大坝抗震的特点在于大坝所处环境复杂(图1),参照世界坝工专家Chopra(1992,2008)的意见,大坝抗震的关键技术包括坝与库水的动力相互作用分析,坝与无限地基的动力相互作用分析,筑坝材料的动态特性及大坝地震损伤发展与风险分析。在2013年中国工程院举办的“中国工程科技论坛——水安全与水利水电可持续发展”会议上,林皋(2014)提出了大坝抗震模型的改进意见。关于坝与库水的动力相互作用分析需要考虑水的可压缩性以及动水压力波在水库边界的吸收作用,Chopra(2008)采用有限元法(FEM)进行求解,动水压力按坝的振动模态展开,计算复杂、难以实际应用。Lin等(2007a,2012a)采用比例边界有限元法(SBFEM)进行求解,按库水的模态进行展开,计算简便、精度高、便于实际应用。关于坝与无限地基的动力相互作用,由于计算模型比较复杂,文献中一般都将无限地基简化为均质地基进行考虑,难以反映实际情况。我们将无限地基分为层状分布、分块不均质和近场杂乱不规则等几种类型,分别提出了不同的计算模型(Lin et al,2015,2007b; Yin et al,2013),计算结果都具有较好的精度与效率。为了提高坝—库水—地基系统时域分析的计算精度与效率,Lin等(2012b)提出以等几何分析方法(isogeometric analysis,简称IGA)代替FEM来进行坝的应力分析,并将IGA与SBFEM相结合来进行坝与库水和坝与无限地基的动力相互作用分析。关于坝材料的动力特性方面,混凝土是速率敏感性材料。自从第一颗原子弹1945年在日本广岛爆炸以后,为了防卫核爆炸造成的灾难,许多国家都加强了混凝土材料抗击核冲击的研究。但这种研究偏重单调瞬时冲击,而且抗压性能居多,抗拉性能偏少。地震荷载的特点是不规则的循环往复,而且加载速率也偏低。为此,我们开展了2000多试件混凝土在地震作用下的动态特性的研究,取得了一批有价值的成果(Yan,Lin,2006; Lin et al,2007c)。
笔者认为除了计算模型以外,数值计算方法也对大坝抗震分析的计算精度与效率产生影响。大坝地震响应分析中,有限元方法的应用比较普遍。有限元法可以适应各种不规则的几何边界条件和分布不规律的单元材料特性,可以应用于各种类型物理问题的求解,计算能力强,不受问题求解自由度的限制。但有限元法也有其不足之处,其主要采用低阶插值函数,连续性差,应力计算精度低。但在工程设计中应力是对大坝安全性评价的重要依据,所以,有必要寻求更为有效的其它数值计算方法。根据我们研究,提出以下计算模型和方法供参考,欢迎批评指正。
1 多边形单元技术
多边形单元(PE-Polygon Element)不是一个新概念,但经过澳大利亚新南威尔士大学(UNSW)宋崇民教授团队的发展,建立在SBFEM基础上的多边形单元已经成为与FEM一样具有广泛适应性的网格离散技术,对数值计算具有强大的功能。PE可称为广义的FE。
大连理工大学发展了PE技术在结构静动力分析中的应用,主要工作为:(1)单元细化与连续性升阶可以只在PE的边界上进行,不改变网格的拓朴结构,使应力计算精度显著提高,但计算工作量则很节约;(2)将PE与曲线边界相结合,既可构造纯直线边的PE,也可构造曲线边的PE,或直线边与曲线边相结合的PE;(3)将PE应用于结构的奇异应力场分析,可以达到很高的计算精度。
多边形单元按SBFEM进行建模(图2)。单元域内选择一点O作为比例相似中心,其要求是从O点可以看见边界上各点,或O点可与边界上任意点用直线连接,没有障碍。单元各边上可以任意划分子单元和结点,数目不限。如图2所示,PE为六边形,每边可以只划分一个子单元,每个子单元含3个结点,则具有2阶精度。AB边也可以划分2个子单元,每个子单元各3个结点共5个结点,这就相当于网格加密。各子单元用SBFEM进行建模,由于SBFEM为半解析方法,从O点至边界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单元位移场的连续性。
Fig.3 Polygon element with curved boundary
(a)typical PE element;(b)interpolation basis
function of AF with curved boundary
PE由域内三角形单元形心点连线加周边三角形单元中点连线构成(图4),所以建模非常方便。单元便于细分和拼合,如图5所示。
(a)任意物理域;(b)多边形网格Fig.4 Modeling of polygon elements
(a)arbitrary physical domain;(b)polygon grid
Fig.5 Refinement and combination of polygon element
(a)initial polygon grid;(b)refinement of polygon mesh;(c)merging of polygon grid
多边形单元具有以下特点:(1)适于各种几何形体的建模,具有高度灵活性;(2)单元具备单位分解性与线性完备性,通过分片检验;(3)可以将网格加密进行单元细化,也可在不改变网格拓扑结构的条件下,在单元各边上增加结点实现单元细化和连续性升阶;(4)可以方便地进行粗细网格过度,只需在不同网格交界面上增加结点即可(图6);(5)可以和FE和NURBS兼容。
mesh without introducing any additional complexity
2 数值算例
为了检验PE的效果,通过以下算例加以说明。
2.1 厚壁圆筒承受均匀内压作用本例有解析解。利用对称性,只取厚壁圆筒的1/4进行离散。假设内压力p=1,径向与环向应力只和内径和外径的比值有关,假设r/R=1/3。将PE与FE(二阶单元)的计算结果进行对比,两种方法的网格划分如图7所示。两种方法计算精度与解析解相比的误差随网格结点数的变化的对比如图8所示,图中PE代表直线边界网格,CPE代表内外曲线边界网格,由图可见PE的计算优越性是明显的。
2.2 曲线形溢流坝受静水压力作用某工程溢流坝,高H=66 m,受静水压力作用。坝材料参数假设为:弹性模量E=36.5 GPa; 泊松比γ=0.167。采用PE与FE两种方法对该工程溢流坝进行计算对比。坝体两种方法网格划分如图9所示,廊道部位网格划分如图 10所示。图 10图题中p表示单元离散所用的插值阶数。廊道上部曲线边界的应力分量对比如图 11所示,廊道周边第一主应力分布如图 12所示。由图可见,有限元网格即使在廊道角点O附近细分,应力奇异性仍无法反映。
Fig.9 Mesh generation for overflow dam
(a)PE model;(b)finite element model(mesh4);(c)the local map of dam crest;(d)the local map of dam crest
(e)PE mesh1(p=2);(f)PE mesh2(p=3);(g)PE mesh3(p=4)
Fig.10 Mesh generation around the gallery
Fig.12 Singular stress distribution around the corner of the galleary
(a)the whole distribution;(b)distribution along the vertical
free surface;(c)distribution along the horizonal free surface
3 异应力场分析
有限元方法无法处理应力奇异性问题。SBFEM由于径向具有解析解,所以特别有利于处理应力奇异性问题。关键是须将应力奇异点设置
于相似中心处。但是SBFEM中,通过相似中心的边界不能进行离散,只能将此边界上作用的外力或变形转化为相应的边界条件,在计算方程中加以体现。我们通过算例来加以说明。
研究坝踵点在水压力作用下的奇异应力分布。图 13表示PE和FE两种网格划分。
(f)FEM mesh1;(g)FEM mesh2;(h)FEM mesh3
Fig.14 Various mesh generation for evaluating singular stress field around dam heel
图 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边上的水压力也对奇异应力场影响比较大的缘故。因此,为了取得良好的奇异应力场的效果,建立边界面的准确边界条件仍然很重要。
and SBFEM(a),between PE and FEM(a)
动力问题奇异应力场的分析也可采用相同的方法进行处理。图 16~20为地震波作用下坝踵点奇异应力场的分析结果。采用的材料参数列入表1,为了表示PE的优越性,将PE结果与比较准确的IGA的计算结果进行比较。IGA的各级网格均可自动加密,网格采用2阶单元。图 18表示奇异应力场分析结果的比较。由图可见,在距坝踵点水平距离x=1.0 m以外,各种网格的计算结果都收敛而且相符性比较好。但距离坝踵点x=0.5 m范围以内,IGA1-IGA3偏差较大,IGA4结果与本文方法吻合良好,达到收敛要求。图 19为坝踵点的位移与加速度响应时程比较。由于位移不存奇性,与应力结果不同,各网格计算结果相差不大。
Tab.1 Calculation parameters for the modelsfor the inputting seismic wave at dam foundation
Fig.18 Results and comparison of thesingular stress field
point of x=0.5 m(a)and x=1.0 m(b)
4 结语
按比例边界有限元法(SBFEM)原理建立的多边形单元(PE)是一种提高应力计算精度的高效数值方法。多边形单元可以进一步发展为多面体单元,即将二维问题发展到三维问题,具有广阔的发展前景。
- 林皋.2014.混凝土高坝抗震分析的新技术[C]//中国工程科技论坛:水安全与水利水电可持续发展.北京:高等教育出版社,49-79.
- Chopra A K.1992.Earthquake analysis,design,and safety evaluation of concrete arch dams[Z]//Proceedings of the Tenth World Conference on Earthquake.Rotterdarm:Balkema A A.publisher,6763.
- Chopro A k.2008.Earthquake analysis of arch dams:factors to be considered[C]//The 14WCEE.Bejing,China.
- Lin G,Du J G,Hu Z Q.2007a.Dynamic dam-reservoir interaction analysis including effect of reservoir boundary absorption[J].Science in China Series E:Technological Sciences,50(1):1-10.
- Lin G,Du J,Hu Z.2007b.Earthquake analysis of arch and gravity dams including the effects of foundation inhomogeneity[J].Frontiers of Architecture and Civil Engineering in China,1(1):41-50.
- Lin G,Han Z,Li J.2015.General formulation and solution procedure for harmonic response of rigid foundation on isotropic as well as anisotropic multilayered half-space[J].Soil Dynamics and Earthquake Engineering,70:48-59.
- Lin G,Wang Y,Hu Z.2012a.An efficient approach for frequency-domain and time-domain hydrodynamic analysis of dam-reservoir systems[J].Earthquake Engineering & Structural Dynamics,41(13):1725-1749.
- Lin G,Yan D,Yuan Y.2007c.Response of concrete to dynamic elevated-amplitude cyclic tension[J].ACI Material Journal,104(6):561-566.
- Lin G,Zhang Y,Wang Y,et al.2012b.A time-domain coupled scaled boundary isogeometric approach for earthquake response analysis of dam-reservoir-foundation system[C]//15th World Conference on Earthquake Engineering Lisbon,Dortugal.
- Yan D,Lin G.2006.Dynamic properties of concrete in direct tension[J].Cement and Concrete Research,36(7):1371-1378.
- Yin X,Li J,Wu C,et al.2013.NSYS implementation of damping solvent stepwise extraction method for nonlinear seismic analysis of large 3-D structures[J].Soil Dynamics and Earthquake Engineering,44:139-152.