发明内容

本发明的目的是提供一种地下水库增湿作用下土-结构相互作用计算方法，该方法能够实现地下水库未渗漏阶段、渗漏阶段的数值计算，为黄土湿陷性地区的非饱和渗流工程问题提供模拟预估，进一步为黄土地区的地下调蓄水库结构设计优化提供理论指导和技术支持。

本发明所采用的技术方案是，一种地下水库增湿作用下土-结构相互作用计算方法，具体按照以下步骤实施：

步骤1：构建调蓄水库混凝土几何结构、土体几何结构模型；

步骤2：设置土体、水库混凝土材料各参数；

步骤3：定义“广义拉伸”、压力函数变量、土的强度指标和模量关于有效饱和度函数变量；

步骤4：设置“固体力学”并设置各边界条件；增加“理查兹方程”并划分网格；

步骤5：设置水体“流体属性”“基本属性”各参数并划分网格；

步骤6：添加“研究”并计算“研究1”；

步骤7：对水库边界再次进行“布尔操作和分割”；

步骤8：设定“理查兹方程”参数；

步骤9：基于“研究1”，使用“辅助扫描”添加水库注水功能；添加并设定“研究2”时间单位、时间步；

步骤10：添加并设定“固体力学”并设置各参数条件并划分网格；

步骤11：基于“研究2”，添加并设定“研究3”；

步骤12：结果计算后数据成图。

本发明的特征还在于，

步骤1具体按照以下步骤实施：

步骤1.1：打开COMSOL软件，在“新建”点击“模型向导”，接下来在“选择空间维度”点击“二维”；

步骤1.2：接上一步，进入到“选择物理场”界面，选择“结构力学”模块，在该模块树中选择并添加物理场接口——“固体力学(solid)”；

步骤1.3：接上一步，接着点击“研究”，在“选择研究”点击“一般研究”，并在该板块树中选择“稳态”；后点击“完成”；

步骤1.4：接上一步，为构建水库、土体几何模型进入下一界面；在模型开发器中选择“几何”，右键“几何”选择“矩形”，该矩形A作为土体结构；在“几何设置”的“大小和形状”中分别设置“宽度”为a，“高度”为b，在“设置”窗口点击“构建选定对象”；再绘制矩形B(水库结构的框架)，右键“几何”选择“矩形”在“几何设置”的“大小和形状”中分别设置“宽度”为c，“高度”为d(上述的a，b为土体模拟尺寸，c，d为水库框架模拟尺寸)；在“图形”窗口栏选择“编辑”将矩形B移动到矩形A中部，在菜单栏选择“几何”中“布尔操作和分割”下的“差集”；然后在“差集”中将“要添加的对象”选为矩形A并激活，将“要减去的对象”选为矩形B并激活(其余选择为默认状态)，点击“构建选定对象”；再回到“几何”右键“矩形”绘制“宽度”为c，“高度”为d的矩形B₂和B大小一致，点击绘制；在“图形”栏选择“编辑”将B₂一定到矩形A的空白位置处(其余选项默认)；采取同上述步骤相同的原理，再绘制宽e，高f的矩形C(e，f是水库内框架，尺寸小于c，d值一个混凝土墙厚的值)，然后将矩形C移动到矩形B的中部，再使用“差集”，用矩形B减去矩形C(取消“保留内部对象”)，并点击“构建选定对象”得到水库内无立柱框架结构D；添加立柱，仍旧右键“几何”选择“矩形”，依次设置若根个宽为g，高度为h的立柱E、F、G…；均匀分布这些立柱在矩形D结构中的位置，再利用“布尔操作和分割”下的“并集”，将水库框架D和立柱E、F、G…作为并集对象(取消“保留内部边界”)，点击“构建选定对象”；以上得到周围土体——域1、水库结构——域2。

步骤2具体按照以下步骤实施：

步骤2.1：在模型开发器中，右键“材料”选择“空材料”修改为“土”，将“几何实体层”“域”设定为“1”并激活；再从材料属性中找到“粘聚力”、“内摩擦角”增加到该“材料属性明细”，依次输入土的“杨氏模量”、“泊松比”、“密度”、“粘聚力”、“内摩擦角”；右键“材料”选择“空材料”修改为“混凝土”，将“几何实体层”“域”设定为“2”并激活；再从材料属性中找到“单轴抗拉强度”“单轴抗压强度”“双轴抗压强度”增加到该“材料属性明细”，依次输入混凝土的“杨氏模量”、“泊松比”、“密度”、“单轴抗拉强度”、“单轴抗压强度”、“双轴抗压强度”(根据实际数据输入参数)；

步骤2.2：在模型开发器中，“几何”树的最后的“形成装配体”并选择全部区域。

步骤3具体按照以下步骤实施：

步骤3.1：在模型开发器中，右键单击“定义”在“组件耦合”下选择“广义拉伸”，在“图形”界面“几何实体层”选择“边界”，将矩形D结构的上下边界作为“源选择”，并将“目标映射”“x，y”均修改为大写；将“源”中“源坐标系”修改为“材料(X,Y,Z)”；再次右键单击“定义”在“组件耦合”下选择“广义拉伸”，在“图形”界面“几何实体层”选择“边界”，将矩形D结构的左右边界作为“源选择”，并将“目标映射”“x，y”均修改为大写；将“源”中“源坐标系”修改为“材料(X,Y,Z)”；

步骤3.2：在“模型开发器”下“组件”“定义”右键选择“变量”，设置侧边和底边的水压力函数，即侧边函数为随着水库高度变化的压力函数——load₁＝H*9.8*(h₁<y<h₂)，并设定水库底部随水位高度变化的竖直向压力，load₂＝water_h*9.8；添加“C＝-15.8*log(dl.Se)+98.3”、 采用插值法，在“全局定义”下选用“插值”，输入t值及其函数f(t),即“f(0.6)＝11000000，f(0.7)＝10000000，f(0.75)＝9500000，f(0.8)＝9000000，f(0.85)＝8500000，f(0.9)＝8000000，f(0.95)＝7500000，f(1.0)＝7000000，f(1.05)＝6900000”；在“组件”的“定义”下的“变量”中写入“E＝Eint(dl.Se)”；在“全局定义”右键“变量”设定“H＝(y-h₁)*(h₁<y<h₂)”和water_h＝0(上式中h₁是水库底边注水起始位置，h₂是水库注水上限，water_h为水库储水高度)。

步骤4具体按照以下步骤实施：

步骤4.1：“模型开发器”下“固体力学”，“线弹性材料”右键选择“土壤塑性”、“混凝土”；其中，“土壤塑性”选择“域1”部分，“混凝土”选择“域2”部分，“屈服准则”采用“Drucker-Prager”并选择“匹配莫尔-库伦准则”，内聚力和内摩擦角均选自“来自材料”；

步骤4.2：在菜单栏“物理场”“域”中选择“重力”其作用域为全部域；“边界”中选择“固定约束”作用边界为域1的下边界；并选择“对称”作用边界为域2的左右边界；在“边界”中选择“指定位移”，依次作为x、y方向的指定位移；添加第一组“指定位移”，边界选择域2的左右外边界，在“设置”中“指定位移”勾选“在x方向指定”并写入前面“定义”中的广义拉伸的算子名称“genext1(x)”；添加第二组“指定位移”边界选择域2的上下外边界，在“设置”中“指定位移”勾选“在y方向指定”并写入前面“定义”中的广义拉伸的算子名称“genext1(y)”；

步骤4.3：在菜单栏选择“物理场”中的“边界荷载”，添加第一组“边界荷载”选择水库内部的沿高度方向的内壁，并在“力”中添加前面“定义”中的“load₁”，添加第二组“边界荷载”选择水库内部的底部内壁，并在“力”中添加前面“定义”中的“load₂”；

步骤4.4：在菜单栏中选择“物理场”“添加物理场”，在“流体流动”接口选择“多孔介质和地下水流”中的“理查兹方程”。

步骤5具体按照以下步骤实施：

步骤5.1：在“理查兹方程”下，“设置”中，选定域1、域2；在“流体属性”下，将密度设为“1000kg/m³”；在“基本属性”设置“多孔材料”为“土体”材料，“饱和液体体积分数”、残余液体体积分数，根据该地区地质勘察报告获取；“渗透率模型”选定“水力传导率”(各向同性)；根据地质勘察报告获；并在“储水模型”中选择“线性储水”，设定“流体的可压缩性”为“4.4e-10[1/Pa]”，“基体有效压缩率”为“10e-6[1/Pa]”；“滞留模型”选定“vanGenuchten”，“本构关系常数”的“α、n、l”根据勘察报告获取；

步骤5.2：定位到“模型开发器”的“初始值”；定位到“初始值”的“设置”，选择全部区域；选定“压力头”设定“H_p＝-0.1m”(该数值可以根据实际调整，-0.1m为举例)；定位到“模型开发器”下的“理查兹方程”右键选择“压力头”。

步骤6具体按照以下步骤实施：

基于上述“压力头”设定；回到“几何”设置；计算“固体力学”下的“研究1”；定位到菜单栏的“结果”，选定“二维绘图组”，右键单击选择“表面”，定位到“设置”，选定“数据”的“数据集”为“固体力学”，点击“模型开发器”的“表面”，定位到“设置”中的“表达式”，输入“solid.pm”，点击“绘制”；得到立柱在受周围土体重力、自身重力的弹性受力下的应力云图，获得立柱底部受力分布图，大概量测得到立柱底端受力从坡段到平缓过渡的水平投影长度，将该长度作为前“理查兹方程”立柱渗流区域的入渗端；

步骤7具体按照以下步骤实施：

再定位回“模型开发器”的“几何”部分，在菜单栏选择“布尔操作和分割”下的“分割边”；回到图形区域，基于前面选定的渗流投影区长度，分割出对应长度的渗流端，分别获得各立柱下的渗流边(例如立柱下渗流区、水库两角点处渗流区)；

步骤8具体按照以下步骤实施：

定位回“理查兹方程”，右键单击“模型开发器”下的“理查兹方程”选择“压力头”，定位到“设置”部分的“边界选择”选择上一步设置好的渗流边，在“压力头”处输入“H_p0＝0.2m”(0.2m为举例，可根据实测压力头输入该值)。

步骤9具体按照以下步骤实施：

步骤9.1：定位回“固体力学”，在对应计算“研究1”；定位到“设置”，在“物理场和变量选择”只勾选“固体力学”；勾选“研究扩展”下的“辅助扫描”；选定“water_h”设定“参数值列表”即“range(h₁,1,h₂)”(h₁是水库底边注水起始位置，h₂是水库注水上限)；单击“计算”；

步骤9.2：定位到菜单栏：单击“研究”选择“添加研究”，选定“一般研究”下的“瞬态”；定位到“设置”，选择“时间单位”为“d”，设置“时间步”为“range(0,3,90)”(本例按照三天为一个计算周期，计算三个月的渗流场)；在“物理场和变量选择”只勾选“理查兹方程”；在“因变量值”下的“设置”选择“用户控制”，“方法”选择“解”，“研究”选择“研究1”，“参数值water_h”选择“自动”，其余默认；单击“计算”。

步骤10具体按照以下步骤实施：

步骤10.1：定位回菜单栏的“物理场”，选择“添加物理场”，在“结构力学”下选择“固体力学”；添加该物理场，并重命名为“固体力学(二)”；定位到“模型开发器”的“固体力学(二)”，右键“线弹性材料”，选定“土壤塑性”、“混凝土”、“外部应力”；定位到“线弹性材料”，将“域1、2”选中作为域选择，“杨氏模量”选用“用户定义”输入“E”，其余默认；定位到“土壤塑性”将区域选择为“域1”，“屈服准则”选定为“Drucker-Prager准则”，并勾选“匹配莫尔-库伦准则”，“内聚力”选定“用户定义”填入“C”，“内摩擦角”选用“用户定义”填入其余默认；定位到“混凝土”选定“域2”，“混凝土准则”选用“Bresler-Pister准则”，“单轴抗拉强度”、“单轴抗压强度”、“双轴抗压强度”均选择“来自材料”(根据实际材料设定)；定位到“外部应力”，选择“域1”，“应力输入”选定“孔隙压力”，“绝对压力”选定“用户定义”并输入“p*(p>0)”，其余默认；

步骤10.2：定位回“模型开发器”，右键“固体力学(二)”，选定“指定位移”，并同样设定“对称”、“重力”，同前述“固体力学”一样设置步骤。

步骤11具体按照以下步骤实施：

定位回菜单栏，点击“研究”并“添加研究”，选定“一般研究”下的“稳态”；定位到“设置”的“物理场和变量选择”，勾选“修改研究步骤的模型配置”，只禁用“固体力学”；定位到“因变量值”在“求解变量的初始值”下“设置”选定“用户控制”，“方法”选定“解”，研究选定“研究2”(即理查兹物理场对应的研究)，“解”选择“解2”，“时间(d)”选择“30d”(举例为30d)；点击“计算”；

以上计算基于相同网格划分，如图2所示。

步骤12具体按照以下步骤实施：

按照上述方案，后处理绘制土体的压力、杨氏模量、内摩擦角、杨氏模量、粘聚力的表面分布图；结构的压力表面图，结构的位移变化曲线。

本发明的有益效果是：本发明在COMSOL软件平台，建立简化的二维平面内土-结构联结物理模型，实现水库不同储水高度下结构的力学响应；建立土体中采用粘聚力、杨氏模量、内摩擦角分别关于有效饱和度的函数作为固体力学的相应参数即考虑增湿作用下土-结构相互作用计算，对结构在“破坏-安全情况”即水库基底破坏发生渗漏的情况下进行渗流模拟，实现了结构破坏、土体渗流后水分场和力学场的耦合作用下结构构件的力学响应。该发明提出的该种计算方法实现地下水库未渗漏阶段、渗漏阶段的数值计算，为黄土湿陷性地区的非饱和渗流工程问题提供模拟预估，进一步为黄土地区的地下调蓄水库结构设计优化提供理论指导和技术支持。