一种大坝温度相场位移的估测方法及计算机可读介质

本发明属于水利水电工程领域，尤其涉及一种大坝温度相场位移的估测方法及计算机可读介质。

背景技术：

1、我国大坝主要以混凝土坝为主，大坝混凝土是一种典型的准脆性材料，由于混凝土大坝易遭受到温度作用的影响，比如大体积混凝土浇筑过程中产生大量水化热导致早期混凝土出现开裂；其次在运行过程中坝体表面混凝土易遭受到寒潮作用的影响导致混凝土表面出现裂缝等，大量裂缝的出现则会对工程造成安全隐患。而目前尚未有一种成熟的技术手段来预测混凝土坝在温度作用下的裂缝发展，传统的裂缝检测技术均是在大坝出现裂缝后进行人工视察和无人机拍照，检测较为繁琐，并且耗时耗力，也无法预知裂缝扩展路径。因此需发明一种方法对混凝土坝在温度变化下是否出现裂缝以及裂缝如何扩展进行预报，进而根据裂缝的预测结果采取合理的工程措施防止安全事故发生。

技术实现思路

1、根据上述研究背景，本发明提供了一种大坝温度相场位移的估测方法及计算机可读介质，能模拟复杂应力条件下的混凝土热开裂情况。

2、本发明方法的技术方案为一种大坝温度相场位移的估测方法，其特征在于：

3、获取水利大坝截面的长度、宽度；

4、在水利大坝截面的长度、宽度范围内，通过实验法获取参数；

5、将水利大坝截面的长和宽进行单元划分得到多个四边形单元、多个节点，构建已知节点的位移边界条件、已知节点的力学荷载边界条件、已知节点的温度边界条件、已知节点的热流荷载边界条件；

6、采用一阶形函数对每个四边形单元内部的位移场、温度场、相场依次进行描述，结合能量最小原理建立热-力-相场耦合模型的控制方程和边界条件所满足的弱积分形式，计算每个四边形单元的刚度矩阵、每个四边形单元的右端矢量；

7、采用牛顿迭代法和欧拉法求解每个时刻每个节点的温度值、相场值、位移值，将相场结果即裂缝带入到大坝安全评价系统对大坝进行安全评估。

8、本发明方法具体包括以下步骤：

9、步骤1：获取水利大坝截面的长度、宽度；

10、步骤2：在水利大坝截面的长度、宽度范围内，通过钻芯取样的方式得到水利大坝截面样品，将水利大坝截面的样品采用抗剪实验获取水利大坝截面样品的抗剪强度τs，采用抗拉实验获取水利大坝截面样品的ft、水利大坝截面样品的弹性模量e和水利大坝截面样品的泊松比ν，采用三点弯曲梁实验获取水利大坝截面样品的断裂能gf，采用应变法测量实验获取水利大坝截面样品的热膨胀系数κ，采用热板法实验获取水利大坝截面样品的导热系数k，采用差式扫描量热实验实验获取水利大坝截面样品的比热容c，采用排水法测量实验获取水利大坝截面样品的密度ρ；

11、步骤3：将水利大坝截面的长和宽进行单元划分得到多个四边形单元、多个节点，构建已知节点的位移边界条件、已知节点的力学荷载边界条件、已知节点的温度边界条件、已知节点的热流荷载边界条件；

12、步骤4：定义每个四边形单元，采用一阶形函数对每个四边形单元内部的位移场、温度场、相场依次进行描述，采用能量最小原理建立热-力-相场耦合模型的控制方程和边界条件所满足的弱积分形式，将已经获得的每个单元材料参数和单元每个节点的坐标带入到弱积分形式得到每个四边形单元的刚度矩阵、每个四边形单元的右端矢量；

13、步骤5：求出每个四边形单元的刚度矩阵和荷载向量，按照节点对应关系形成所有节点的总刚度矩阵和荷载向量，引入步骤3的约束条件采用牛顿迭代法和欧拉法求解每个时刻每个节点的温度值、相场值、位移值，将相场结果即裂缝带入到大坝安全评价系统对大坝进行安全评估；

14、作为优选，步骤3所述将水利大坝截面的长和宽进行单元划分得到多个四边形单元、多个节点，具体如下：

15、长度方向划分为a份，宽度方向划分为b份，则可得到nelem＝a*b个四边形单元和nnode＝(a+1)(b+1)个节点；

16、步骤3所述已知节点的位移边界条件，定义如下：

17、若lable(xj,yj)∈nu,则ux(xj,yj)＝ux0，uy(xj,yj)＝uy0

18、其中，lable(xj,yj)为已知节点j的编号，j∈[1,nnode]，xj和yj为已知节点j的横坐标和纵坐标，其满足xj∈[0,m]，yj∈[0,n]；

19、步骤3所述已知节点的力学荷载边界条件，定义如下：

20、若lable(xj,yj)∈nf,loadx(xj,yj)＝fx0，loady(xj,yj)＝fy0

21、步骤3所述已知节点的温度边界条件，定义如下：

22、若lable(xj,yj)∈nθ,

23、第k个时刻已知节点的温度为θ(k,xj,yj)＝t(t),k∈[1,k]

24、步骤3所述已知节点的热流荷载边界条件，定义如下：

25、若lable(xj,yj)∈nq，q(xj,yj)＝q0

26、m为大坝的长度，n为大坝的宽度；

27、已知位移的节点编号组成的集合为nu,已知力学荷载的节点编号组成的集合为nf，已知温度的节点编号组成的集合为nθ,已知热流荷载的节点编号为nq，所有节点编号组成的集合为o；

28、作为优选，步骤4定义每个四边形单元，具体如下：

29、第ii行和第jj列个单元由四个节点构成，其中ii∈(1,a)，jj∈(1,b)，四个节点编号记为nlb，nrb，nrt，nlt，即左下，右下，右上，左上，且满足nlb,nrb,nrt,nlt∈[1,b]，则单元的质心坐标为：

30、

31、其中xcenter，ycenter分别为单元质心的横坐标和纵坐标。

32、则四个顶点的坐标则可通过如下表达式确定：

33、

34、

35、

36、

37、其中(xnlb,ynlb)(ii,jj)，(xnrb,ynrb)(ii,jj)，(xnrt,ynrt)(ii,jj)，(xnlt,ynlt)(ii,jj)分别为左下，右下，右上，左上节点的坐标，m为水利大坝截面的长，a为长边划分的份数，n为水利大坝截面的宽，b为宽所在边划分的份数；

38、步骤4所述一阶形函数，定义如下：

39、

40、

41、其中n1，n2，n3，n4为四个形函数，ξ和η为单元局部坐标系，以单元的质心为原点，两边的方向向量分别为横轴与纵轴，与笛卡尔坐标系存在如下转换关系：

42、

43、其中为雅可比矩阵，维数为2行2列，n1，n2，n3，n4为四个形函数，为形函数关于局部坐标系ξ方向的导数，为形函数关于局部坐标系η方向的导数，为形函数关于笛卡尔坐标系x方向的导数，为形函数关于笛卡尔坐标系y方向的导数；

44、步骤4所述对每个四边形单元内部的位移场进行描述，具体如下：

45、

46、其中ux，uy表示单元内部的x方向位移场以及y方向位移场，n1，n2，n3，n4为四个形函数，uxnlb，uynlb，uxnrb，uynrb，uxnrt，uynrt，uxnlt，uynlt分别为四个节点的x方向位移和y方向位移；

47、步骤4所述对每个四边形单元内部的温度场进行描述，具体如下：

48、

49、其中θ为单元内部的温度场，n1，n2，n3，n4为四个形函数，θnlb，θnrb，θnrt，θnlt分别为四个节点的温度；

50、步骤4所述对每个四边形单元内部的相场进行描述，具体如下：

51、

52、其中s为单元内部的相场，n1，n2，n3，n4为四个形函数，snlb，snrb，snrt，snlt分别为四个节点的相场值；

53、步骤4所述采用能量最小原理建立热-力-相场耦合模型的控制方程和边界条件所满足的弱积分形式，具体如下：

54、所述建立热-力-相场耦合模型的控制方程如下：

55、

56、其中，σ为应力张量；s为相场标量，s＝0表示没有破坏，s＝1表示完全破坏，l0为裂纹宽度特征参数，为单元尺寸大小的5-10倍，h为最大裂纹驱动力，ω(s)为能量退化函数，α(s)＝2s-s2为裂纹几何函数，ρ和c为材料的密度和比热容，θ为温度标量，qint为热源，k为热传导张量。

57、所述建立热-力-相场耦合模型的边界条件如下：

58、

59、其中,t是边界上的面力，u是位移场；为裂纹密度函数；为热流，为给定的边界热流；

60、所述弱积分形式为：

61、

62、其中，nu，nθ，ns分别为位移、温度、相场形函数矩阵，转置矩阵分别为，bu，bθ，bs分别为位移、温度、相场形函数导数矩阵，转置矩阵分别为σ为应力张量，为单元内部的体积力，如重力，为边上的面力，如水压力，j为热流张量，qint为单元内部的热源，为边上的热流，ω″(s)为能量退化函数二阶导数，为裂纹驱动力，具体表达式如后文所示，α″(s)为裂缝几何函数二阶导数，为-2，l0为裂纹宽度特征参数，为单元尺寸大小的5-10倍，

63、步骤4所述将已经获得的每个单元材料参数和单元每个节点的坐标带入到弱积分形式可得到一个单元的刚度矩阵和右端矢量，具体如下：

64、将本构关系代入到控制方程的弱积分形式并进行化简计算可得到如下刚度矩阵和右端矢量：

65、

66、其中，kuu，kus，kuθ，ksu，kss，kθθ分别为位移刚度矩阵，位移相场刚度矩阵，位移温度刚度矩阵，相场位移刚度矩阵，相场刚度矩阵，温度刚度矩阵，fu，fs，fθ为位移、相场、温度右端矢量。

67、为简化计算，令kus＝ksu＝0。

68、各项刚度矩阵表达式如下：

69、

70、其中，nu为2行8列的矩阵，为位移形函数矩阵，表达式如下：

71、

72、其中n1，n2，n3，n4为四个形函数。

73、其中，bu为3行8列的行列式，为位移形函数导数矩阵，表达式如下：

74、

75、其中，为形函数关于笛卡尔坐标系x方向的导数，其中为形函数关于笛卡尔坐标系y方向的导数。

76、其中，为材料的弹性矩阵，为与材料的弹性模量e和泊松比ν相关的矩阵，维数为3行3列，表达式如下：

77、

78、其中，e，ν分别为材料的弹性模量和泊松比。

79、其中，ω(s)为能量退化函数，其表达式为:

80、

81、其中，s为相场标量，e，gf，ft分别为材料的弹性模量，断裂能和抗拉强度，c0＝π，l0为相场特征参数，为单元尺寸的5-10倍。

82、

83、其中，κ为材料的热膨胀系数，e，ν分别为材料的弹性模量和泊松比。

84、其中，ns为1行4列的矩阵，为相场形函数矩阵，其表达式如下：

85、ns＝[n1 n2 n3 n4]

86、其中，bs为2行4列的矩阵，为相场形函数导数矩阵，其表达式过于复杂没有显示的具体给出，但是是为η和ξ的函数

87、

88、其中，为形函数关于笛卡尔坐标系x方向的导数，其中为形函数关于笛卡尔坐标系y方向的导数。

89、其中，nθ为1行4列的矩阵，为温度形函数矩阵，其表达式如下：

90、nθ＝[n1 n2 n3 n4]

91、其中，n1，n2，n3，n4为四个形函数。

92、其中，bθ为2行4列的矩阵，为温度形函数导数矩阵，其表达式过于复杂没有显示的具体给出，但是是为η和ξ的函数

93、

94、其中，为形函数关于笛卡尔坐标系x方向的导数，其中为形函数关于笛卡尔坐标系y方向的导数；

95、其中ρ，c，k分别为材料的密度和比热容以及导热系数；

96、右端矢量表达式如下：

97、

98、其中为位移形函数转置矩阵，为温度形函数转置矩阵，为单元内部的体积力，如重力，为边上的面力，如水压力，qint为单元内部的热源，为边上的热流。

99、将第ii行jj列单元的刚度矩阵和右端矢量采用高斯插值公式进行计算可得如下表达式：

100、

101、

102、各矩阵和函数以及参数定义如前文所述，采用四个高斯积分点对刚度矩阵表达式进行求和可得到每个矩阵的数值大小。

103、其中4个高斯积分点对应的ξ和η的组合如下：

104、

105、

106、

107、

108、其中(ξ,η)i,i∈(1,4)为高斯积分点在局部坐标系的坐标。

109、其中hi,i∈(1,4)为四个高斯积分点的权值系数，且h1＝h2＝h3＝h4＝1。

110、将四组ξ和η的组合带入到刚度矩阵表达式，这样即可得到的每项的数值，是一个8行8列的矩阵。是一个8行4列的矩阵，是一个4行4列的矩阵，是一个4行4列的矩阵；

111、步骤4所述每个四边形单元的刚度矩阵、每个四边形单元的右端矢量，具体如下：

112、得到一个单元刚度矩阵大小为16行16列的矩阵，待求变量为4个节点两个方向的位移值ux和uy，相场值s和温度值θ，右端矢量为监测到的力学荷载和热流荷载；

113、

114、其中，分别为节点nlb对应的ux，uy，s，θ刚度矩阵数值大小，uxnlb，uynlb，snlb，θnlb分别为节点nlb对应的ux，uy，s，θ数值大小，fxnlb，fynlb，fθnlb分别为节点nlb对应的x方向力学荷载，y方向力学荷载，以及热流荷载；

115、作为优选，步骤5所述求出每个四边形单元的刚度矩阵和荷载向量，具体如下：

116、

117、其中，分别为节点nlb对应的ux，uy，s，θ刚度矩阵数值大小，uxnlb，uynlb，snlb，θnlb分别为节点nlb对应的ux，uy，s，θ数值大小，fxnlb，fynlb，fθnlb分别为节点nlb对应的x方向力学荷载，y方向力学荷载，以及热流荷载；

118、步骤5所述按照节点对于关系形成所有节点的总刚度矩阵和荷载向量，具体如下：

119、

120、其中分别为节点kk对应的ux，uy，s，θ刚度矩阵数值大小，uxkk，uykk，skk，θkk，kk∈(1,nnode)分别为节点kk对应的ux，uy，s，θ数值大小，fxkk，fykk，fθkk，kk∈(1,nnode)分别为节点kk对应的x方向力学荷载，y方向力学荷载，以及热流荷载；总刚度矩阵得行和列均为4*nnode，所求变量矩阵行为4*nnode，列为1，右端矢量矩阵得行数为4*nnode，列为1；

121、步骤5所述引入步骤3的约束条件采用牛顿迭代法和欧拉法求解每个时刻每个节点的温度值、相场值、位移值，具体如下：

122、如果kk∈nu，则uxkk＝ux0，uykk＝uy0，

123、如果kk∈nf，则fxkk＝fx0，fykk＝fy0；

124、如果kk∈nθ，则θkk(k)＝t(t)，k∈(1,k)，

125、如果kk∈nq，则fθkk＝q0；

126、其中nu，nf，nθ，nq为已知位移，已知力学荷载，已知温度，已知温度荷载的节点编号集合，ux0，uy0，fx0，fy0，t(t)，q0为已知节点x方向位移值，y方向位移值，x方向力学荷载，y方向力学荷载，温度和热流荷载，uxkk，uykk，θkk，kk∈(1,nnode)分别为节点kk对应的ux，uy，θ数值大小，fxkk，fykk，fθkk，kk∈(1,nnode)分别为节点kk对应的x方向力学荷载，y方向力学荷载，以及热流荷载，为节点kk对应的ux，uy，θ的刚度矩阵数值大小；

127、将传感器所在节点得到的温度值和位移值带入到方程的未知变量矩阵，同时将对应的刚度矩阵乘以106，将监测到的力学荷载和热流荷载大小带入到右端矢量矩阵，采用牛顿迭代法和欧拉法求解每个时刻下所有节点的位移值、相场值和温度值；

128、构建出新的裂纹驱动力表达式，改进的能量密度函数可以表示为：

129、

130、其中，ωι(s)、ωп(s)分别为受拉部分和受剪部分的能量退化函数，分别为拉伸、剪切和压缩应变能，为简化起见，采用同一退化函数进行能量退化，为初始弹性应变能，可表示为：

131、

132、其中λ和μ为材料的拉美常数，可通测量得到的弹性模量e和泊松比ν得到：

133、

134、

135、拉伸应变能和剪切应变能可表示为：

136、

137、

138、其中，分别为破坏面上的正应力和剪应力，其大小与破坏面与有效主应力的夹角β有关系，e为材料的弹性模量；

139、根据统一破坏准则可得如下关于破坏角β的函数：

140、

141、其中为破坏面上的正应力和剪应力，ft，τs分别为材料的抗拉强度和抗剪强度，分别为应力圆的圆心和半径，和为有效应力张量的第一主应力和和第二主应力，β为破坏面和有效主应力间的夹角，对于空间任意一点，其断裂面发生在f(2β)≥1处，所以破坏面的方位角可通过求解下式关系得到：

142、

143、对f(2β)求导得到如下式子：

144、

145、然后得到如下五种情况，令χ＝τs/ft为剪拉比；

146、情况1：

147、f′(2β)＞0，然而f′(2β)＜0，则因此：

148、

149、情况2：并且

150、f′(2β)＜0,因此βr＝0。

151、情况3：并且

152、f′(2β)＞0，因此：

153、

154、情况4：并且

155、f′(2β)＞0，f′(2β)＜0，因此：

156、βr＝0or

157、情况5：

158、此情况下，材料点完全处于受压状态，假设在压应力状态下不会导致断裂，因此在此情况下和被采用。

159、采用改进的裂纹搜索模式对此情况下的裂纹方向进行追踪，有以下三种情形：

160、情形1：并且βr＝0；

161、情形2：并且βr＝π/4；

162、情形3：βr＝π/4

163、在求出裂纹面上的正应力和剪应力后则可得到裂纹驱动力的大小：

164、

165、其中为临界裂纹驱动力，为当前裂纹驱动力，gfп为二型断裂能，与材料的断裂能存在关系：

166、由于最大裂纹驱动力由两部分组成，因此还需要引入两个历史变量记录求解过程中的拉伸和剪切驱动力，具体表达式如下：

167、

168、其中为上一时刻的拉伸驱动力，为下一时刻的拉伸驱动力，为上一时刻的剪切驱动力，为下一时刻的剪切驱动力。

169、其中总的驱动力历史变量为：

170、

171、本发明优点在于，本发明能对温度引起的混凝土等准脆性材料的复杂开裂模式提供了一种新的计算方法，此方法能解决传统计算方法无法考虑材料因剪切导致的裂纹演化，为水工混凝土结构在温度变化下是否产生裂缝以及裂缝如何发展提供了指导意义。