1、l 绪绪 论论l 有限元法正分析简要有限元法正分析简要l 线弹性位移反分析线弹性位移反分析l 粘弹性位移反分析粘弹性位移反分析l 工程应用工程应用主要内容：主要内容：5.1 绪绪 论论 随着上世纪六十年代，电子计算机的问世和快速随着上世纪六十年代，电子计算机的问世和快速发展，数值方法成为岩体工程问题分析的主要手段，发展，数值方法成为岩体工程问题分析的主要手段，如何确定本构模型和输入参数就成为这种手段能否成如何确定本构模型和输入参数就成为这种手段能否成功应用的关键，在通过试验手段获得参数比较困难的功应用的关键，在通过试验手段获得参数比较困难的背景下，通过现场测定位移反求地应力和岩体力学参背景下，
2、通过现场测定位移反求地应力和岩体力学参数的数的“反分析方法反分析方法”被提出，经过多年的研究，目前被提出，经过多年的研究，目前已成为岩石力学一个独立分支。已成为岩石力学一个独立分支。 5.1.1 问题产生及发展历史问题产生及发展历史 岩石力学中的反分析最早由岩石力学中的反分析最早由Kavangh（1973）、）、Gioda和和Maier（1980）等人提出，）等人提出，Sakurai（1983）完成了岩体弹性模量和初始地应力的线弹性有限元位完成了岩体弹性模量和初始地应力的线弹性有限元位移反分析，此后又发展了弹塑性、粘弹性、粘塑性等移反分析，此后又发展了弹塑性、粘弹性、粘塑性等非线性位移反分析，
3、并引入了误差分析、优化技术等非线性位移反分析，并引入了误差分析、优化技术等一些手段，以求获得非线性反分析中的最佳值。一些手段，以求获得非线性反分析中的最佳值。5.1.2 正分析与反分析正分析与反分析已知系统的模型已知输入求输出)(tq)(t 1 正分析正分析2 反分析反分析已知系统的模型)(tq求输入已知输出量测值)()(tt )(tq已知输入已知输出量测值已知系统的模型结构 求模型未知参数)()(tt 岩石力学中的反分析主要有以下几种类型：岩石力学中的反分析主要有以下几种类型：1. 已知岩体的本构模型、初始地应力和位移量测值，已知岩体的本构模型、初始地应力和位移量测值，求岩体物理力学参数求岩
4、体物理力学参数;2. 已知岩体的本构模型、物理力学参数和位移量测已知岩体的本构模型、物理力学参数和位移量测值，求初始地应力值，求初始地应力;3. 已知岩体的本构模型、物理力学参数、初始地应已知岩体的本构模型、物理力学参数、初始地应力和位移量测值，求开挖空间最佳几何形状；力和位移量测值，求开挖空间最佳几何形状；4. 已知初始地应力和位移量测值，求岩体的本构模已知初始地应力和位移量测值，求岩体的本构模型及模型参数，即系统辨识。型及模型参数，即系统辨识。3 反分析问题的特点反分析问题的特点多解性、无解性、不稳定性。多解性、无解性、不稳定性。5.1.3 反分析中的几个要素反分析中的几个要素1 模模 型
5、型 模型是模型是 “原型原型”的一种的一种“类似类似”，任何模型都，任何模型都不能反映出原型的一切特征。不能反映出原型的一切特征。 模型的表达形式可以是概念的、物理的或数学的，模型的表达形式可以是概念的、物理的或数学的，用数学描述形式建立的模型为数学模型。用数学描述形式建立的模型为数学模型。2 参数和状态参数和状态 参数是系统的内部状态变量，反映了系统的本参数是系统的内部状态变量，反映了系统的本质，是不可测量的；状态是系统的外部表现，是可质，是不可测量的；状态是系统的外部表现，是可以测量的。以测量的。 在岩石力学数学模型中，因变量，如位移、应力、在岩石力学数学模型中，因变量，如位移、应力、应变
6、均为外部状态变量，弹性模量、泊松比、内粘结应变均为外部状态变量，弹性模量、泊松比、内粘结力等均为参数。力等均为参数。3 准则函数准则函数 由于模型的近似性和量测误差的存在，在已知量由于模型的近似性和量测误差的存在，在已知量和待求量之间对等的情况下，求出的结果往往不能很和待求量之间对等的情况下，求出的结果往往不能很好地反映系统的本质。好地反映系统的本质。 可行的方法就是增加已知量的数量，求待求量的可行的方法就是增加已知量的数量，求待求量的最优值，为此需要引入一个准则函数。最优值，为此需要引入一个准则函数。 准则函数有两类：以量测值为基础的第一类准则准则函数有两类：以量测值为基础的第一类准则函数；
7、以量测误差及其统计特性为基础的第二类准则函数；以量测误差及其统计特性为基础的第二类准则函数。常用准则函数。函数。常用准则函数。21 )()(iinittyJ 常规最小二乘法（常规最小二乘法（OLS法）法） 高斯高斯马尔可夫估计（马尔可夫估计（GM法）法） 最大似然估计（最大似然估计（ML法）法） 贝叶斯估计（贝叶斯估计（MAP法）法） niiiittyJ122 )()( niiypJ1)/( dYpJ)/()(2 5.1.4 反分析求解方法反分析求解方法1 逆法逆法 将模型输出表达成待求量的显函数，与量测值构成准则函将模型输出表达成待求量的显函数，与量测值构成准则函数直接求解。数直接求解。2
8、正法正法 当模型输出不能表达成待求量的显函数时，先给出待求量的当模型输出不能表达成待求量的显函数时，先给出待求量的初值，计算出模型的输出，与量测值一起代入准则函数求出准初值，计算出模型的输出，与量测值一起代入准则函数求出准则函数值，按一定的路径待求量的值，可计算出一系列准则函则函数值，按一定的路径待求量的值，可计算出一系列准则函数值，使得准则函数值达到最小的待求量值即为最优值。数值，使得准则函数值达到最小的待求量值即为最优值。 该方法是由一系列正算过程构成，故名正法。其适用范围较该方法是由一系列正算过程构成，故名正法。其适用范围较逆法更广。逆法更广。 正法中要用到最优化方法，最常用的有模式搜索
9、法、变量轮正法中要用到最优化方法，最常用的有模式搜索法、变量轮换法、单纯形法、鲍威尔法。换法、单纯形法、鲍威尔法。5.2.1 有限单元法的基本思路有限单元法的基本思路 将连续求解域离散为有限个、按一定方式将连续求解域离散为有限个、按一定方式相互连接在一起的单元组合体，在每个单元内相互连接在一起的单元组合体，在每个单元内用一假设的位移函数来表示待求的未知位移场用一假设的位移函数来表示待求的未知位移场函数，而假设的位移函数用单元节点上的未知函数，而假设的位移函数用单元节点上的未知位移来表示，以此可导出单元内以未知节点位位移来表示，以此可导出单元内以未知节点位移所表示的应力、应变，最后通过最小势能原
10、移所表示的应力、应变，最后通过最小势能原5.2 有限元法正分析简要有限元法正分析简要理导出单元每个节点上以未知节点位移表示的理导出单元每个节点上以未知节点位移表示的平衡方程，整个求解域所有节点平衡方程将构平衡方程，整个求解域所有节点平衡方程将构成一方程组，通过求解该方程组可求得各节点成一方程组，通过求解该方程组可求得各节点上的位移，从而求得单元内的位移、应力、应上的位移，从而求得单元内的位移、应力、应变。变。 总之，有限单元法就是将连续无限自由度总之，有限单元法就是将连续无限自由度问题转变为求离散的有限自由度问题，将偏微问题转变为求离散的有限自由度问题，将偏微分方程组的求解转化为代数方程组的求
11、解。分方程组的求解转化为代数方程组的求解。 对求解连续体的边值问题，有限单元法是一种近对求解连续体的边值问题，有限单元法是一种近似方法，近似的程度随单元划分加密而提高，但会带似方法，近似的程度随单元划分加密而提高，但会带来计算工作量的增加。来计算工作量的增加。5.2.2 有限单元法求解的一般过程有限单元法求解的一般过程1. 确定计算模型，包括计算坐标系、模型整体尺寸、确定计算模型，包括计算坐标系、模型整体尺寸、边界条件、计算参数和地应力场；边界条件、计算参数和地应力场；2. 单元划分，平面三角形（单元划分，平面三角形（3、6、7节点）、平面节点）、平面四边形（四边形（4、8、9节点）、四面体、
12、六面体（节点）、四面体、六面体（8、20、21节点）；节点）；3. 选择位移模式（形函数）；选择位移模式（形函数）；yx4. 单元刚度分析；单元刚度分析；5. 总刚度分析；总刚度分析；6. 节点等效载荷分析；节点等效载荷分析；7. 整体平衡方程建立；整体平衡方程建立；8. 已知边界条件引入；已知边界条件引入；9. 求解平衡方程，获得节点位移；求解平衡方程，获得节点位移；10.根据几何和物理方程求单元高斯点应变和应力；根据几何和物理方程求单元高斯点应变和应力；11.计算结果后处理（绘制位移、应力、塑性区曲线计算结果后处理（绘制位移、应力、塑性区曲线 图或等值线图图或等值线图 位移模式位移模式 1
13、88212eN5.2.3 有限单元法正分析基本方程有限单元法正分析基本方程以平面四节点的等参单元为例。以平面四节点的等参单元为例。等参数单元解决了矩形复杂性问题，坐标变换与位移等参数单元解决了矩形复杂性问题，坐标变换与位移变幻的形函数一样。变幻的形函数一样。其中：其中： Tvu, ININININN4321 Tevuvuvuvu44332211,00114/1iN 几何方程几何方程 188313 eB 其中：其中： Txyyx , 4321 BBBBB xNyNyNxNBiiiii00（i=1、2、3、4）根据等参单元的坐标变换式：根据等参单元的坐标变换式： 4141iiiiiiyNyxNx
14、yyNxxNNyyNxxNNiiiiii得：得： yNxNJyNxNyxyxNNiiiiii 其中：其中： 41414141iiiiiiiiiiiiyNxNyNxNyxyxJ 用矩阵表示：用矩阵表示： 4433221143214321yxyxyxyxNNNNNNNN 本构方程本构方程 188333133313 eeBDD 其中：其中： Txyyx , 单元势能分析单元势能分析a 单元应变能单元应变能 eeTeTeeTeTeTeTexyxyyyxxedxdyBDBhdxdyBDBhdxdyDhdxdyhdxdyhV 2 2 2 2 2令令 1111 ddJBDBhdxdyBDBhKTeTe则则
15、eeeTeKV 21 b 体积力势能体积力势能 设体积力设体积力 Tyxeppp, 则体积力势能为：则体积力势能为： 1111 dhdJpNdxdypNhdxdyphTeTeTeTeTc 表面力势能表面力势能 设面力设面力 Tyxeqqq, 则面力势能为：则面力势能为： sTeTsTeTsThdspNdsqNhdsqh d 集中力势能集中力势能 设集中力设集中力 TyxeRRR, 则集中力势能为：则集中力势能为： eeTR e 单元总势能单元总势能 seeTeeTeTeeeTeRhdsqNhdxdypNK) (21 将所有单元势能叠加得系统总势能将所有单元势能叠加得系统总势能 FKTTe 21
16、根据最小势能原理，真实解应使系统势能最小根据最小势能原理，真实解应使系统势能最小即即 0 由此得系统总平衡方程由此得系统总平衡方程 FK 引入位移边界条件，得最终求解方程：引入位移边界条件，得最终求解方程： FK 5.3 线弹性位移反分析线弹性位移反分析5.3.1 反分析基本公式反分析基本公式 以隧道开挖平面应变问题为例。设开挖边界上的以隧道开挖平面应变问题为例。设开挖边界上的初始地应力为：初始地应力为：Txyyx, 将该初始地应力转化开挖边界上的等效节点力：将该初始地应力转化开挖边界上的等效节点力：dshBFsT 或写成：或写成：321BBBFxyyx 在整个求解域上：在整个求解域上：321
17、BBBFxyyx 根据有限元求解基本方程：根据有限元求解基本方程： FK 令：令：*KEK 则：则：321 *BBBKExyyx 假定量测点与有限元网格节点重合，则可把节点位假定量测点与有限元网格节点重合，则可把节点位移分成已知和待求量部分：移分成已知和待求量部分： TNM , 相应的平衡方程写为：相应的平衡方程写为： 23132212211122211211*BBBBBBKKKKExyyxNM 将未知位移消去：将未知位移消去：xyxyyyxxMNBBBKE * 其中：其中：211221211 *BKKBBx 221221212 *BKKBBy 231221213 *BKKBBxy *211*
18、22*12*11*KKKKKN 0 AM 上式可简写为：上式可简写为： 其中：其中： , ,1*1*1*xyNyNxNBKBKBKA TxyyxEEE/,/,/0 因为测量都是两点相对位移，绝对位移与相对位因为测量都是两点相对位移，绝对位移与相对位移之间转化关系：移之间转化关系： 22111212vuvuTvvuu MMT cossincossinsincossincosT T 为转换矩阵：为转换矩阵： 0* AM 则：则： *ATA 其中其中 上式中待求量上式中待求量 为为3个，若量测值个，若量测值 刚刚好为好为3个，则可从上式中求出唯一的个，则可从上式中求出唯一的0 0 MA* 10M M
19、TTAAA*1*0 若量测值若量测值 多于多于3个，则通过最小二乘法个，则通过最小二乘法求得求得 ，构造以下目标函数，构造以下目标函数0 M ) () ()(*MTMAAF 000令令022000MTTAAAF )(* MTMTMTTAAA)( )(* 0002 得得 的最小二乘估计为：的最小二乘估计为：0 5.3.2 实例实例 根据围岩内部位移求得：根据围岩内部位移求得：T0003. 0,0609. 0,0489. 00 根据围岩收敛变形求得：根据围岩收敛变形求得：T0001. 0,0587. 0,0490. 00 隧道埋深隧道埋深400m，可求得：，可求得：Mpa8 . 94000245.
20、 0 Hy 则：则：160.9Mpa0.06099.80 yyE Mpa87. 70 xxE Mpa048. 00 xyxyE 5.4 粘弹性位移反分析粘弹性位移反分析5.4.1 反分析基本公式 ttve 实际测试的位移是一组随时间变化的值，是由于实际测试的位移是一组随时间变化的值，是由于释放荷载随时间的改变及围岩的蠕变造成，采用粘弹释放荷载随时间的改变及围岩的蠕变造成，采用粘弹性反分析可以很好反映时间因素。性反分析可以很好反映时间因素。 任意时刻任意时刻 t 的应变可以看做由瞬时弹性应变和蠕的应变可以看做由瞬时弹性应变和蠕变应变两部分构成：变应变两部分构成： 其中弹性应变：其中弹性应变： 1
21、00CEe 粘弹性应变：粘弹性应变： tCtEtv 10 则：则： tCtEEtEEtCtEEtEEtCtECEttve 11000000000 tDtEEtCtEEtEEt /1101000 得：得： 其中：其中： 为弹性矩阵为弹性矩阵100DCE 任意时刻任意时刻 t 的有限元平衡方程：的有限元平衡方程： 0 tPdtBT 将将 代入上述平衡方程：代入上述平衡方程： t 0/110 tPdtDBtEET tPtEtEEdtUBDBT 0 tPtEtEEtUK 0 tPtEtEEtUKE 0*0 tPEtPtEEtEEtUKt100* tEEtEEEt 00 其中：其中： 称之为综合模量。称
22、之为综合模量。 常用流变模型的粘弹性模量和综合模量常用流变模型的粘弹性模量和综合模量粘弹性模型粘弹性模量综合模量MaxwellKelvinKelvin-VoigtPoynting-ThomsonBurgerst /1 1/100tEE tEE111exp1/ tEE111exp1/ tEEEEEE102102exp1/ tEEEEEE1021012exp1/ tEE111exp1/ tEEEEE110011exp/ tEtEE11211exp1/ tEtEEEE1121011exp1/ 将位移将位移 U( t ) 分为已知和未知两部分分为已知和未知两部分 TnmtUtUtU, 仿照线弹性位移反
23、分析公式推导，得粘弹性反分仿照线弹性位移反分析公式推导，得粘弹性反分析基本方程析基本方程 tAtUm0 tATtUm0 TxyyxtEt0000,1 其中：其中：5.4.2 参数回归与优化 实际量测过程中，设置测点时间为实际量测过程中，设置测点时间为 t0，到，到 t 时刻时刻量测的位移是量测的位移是 t0t 之间发生的位移，为此把综合模量之间发生的位移，为此把综合模量也用也用t0t 之间的变化来代替之间的变化来代替0111tttEEE 设置测点时间t0时间t设从设从 t0t 之间共进行了之间共进行了n次量次量 ，每次量测有每次量测有m个测点，量测结果为：个测点，量测结果为：),321nttt
24、t（ )(,)(,)(,)(321nmmmmtUtUtUtU 则利用每次量测的则利用每次量测的m个测点已知位移，可反求该个测点已知位移，可反求该时刻的综合模量时刻的综合模量 ，n 次量测可以反求出不同时刻次量测可以反求出不同时刻的综合模量：的综合模量：tE nttttEEEE ,321 根据综合模量与模型参数之间的关系，比如根据综合模量与模型参数之间的关系，比如P-T模型：模型：)(10021iatatteeEEEE 1021 EEEa 式中：式中：iiatteEEEEE02111 将上式变换为以下形式：将上式变换为以下形式：两边取对数：两边取对数：itatEEEEEi 021ln)11ln(
25、令：令：021ln,)11ln(EEEbEEyiti 得：得：batyii 用矩阵表示：用矩阵表示：BXY TnyyyyY321 其中：其中： 111nttX TbaB 最小二乘法解为：最小二乘法解为： YXXXBTT1 根据求得的根据求得的a、b、yi 再求再求E0、E1、E2、1 设设 时刻的量测值为：时刻的量测值为：5.5 工程应用工程应用5.5.1 量测数据预处理量测数据预处理 量测数据带有误差，甚至有奇异数据，为此需对量测数据带有误差，甚至有奇异数据，为此需对其进行一定的数学处理，改善其规律性。最常用的就其进行一定的数学处理，改善其规律性。最常用的就是多项式拟合。是多项式拟合。ntt
26、tt,321nuuuu,321 令：令：nntatataatu 2210)( 建立如下目标函数：建立如下目标函数： niiinutuaaaf1210)(),(0)( iiaaf求解方程组求解方程组得得 ai从而得到：从而得到：nntatataatu 2210)(利用回归公式求得任意时刻的利用回归公式求得任意时刻的u(ti)进行反分析。进行反分析。5.5.2 反分析结果检验反分析结果检验 方法一：将反分析得到的参数作为输入，在相同方法一：将反分析得到的参数作为输入，在相同模型下进行计算，通过计算结果与量测结果的比模型下进行计算，通过计算结果与量测结果的比较来衡量反分析结果的正确性。较来衡量反分析
27、结果的正确性。 方法二：将反分析得到的参数作为输入，对尚未方法二：将反分析得到的参数作为输入，对尚未发生的位移进行预测，根据预测的精度来衡量反发生的位移进行预测，根据预测的精度来衡量反分析结果的正确性。分析结果的正确性。5.5.3 测点布置测点布置收敛测线收敛测线内部位移收敛测线 （a） 3条收敛测线 （b） 6条收敛测线 （c）收敛测线加内部位移5.5.4 围岩变形预测围岩变形预测 利用小规模实验工程量测结果进行反分析，再用反分析结利用小规模实验工程量测结果进行反分析，再用反分析结果对实际工程围岩变形进行预测。果对实际工程围岩变形进行预测。 利用正在施工工程的量测结果进行反分析，在利用反分析利用正在施工工程的量测结果进行反分析，在利用反分析结果对施工效果和后期围岩变形、支护状态进行预测。结果对施工效果和后期围岩变形、支护状态进行预测。
