地下水库增湿作用下土-结构相互作用计算方法与流程

地下水库增湿作用下土
‑
结构相互作用计算方法
技术领域
1.本发明属于多场耦合下土
‑
结构相互作用计算方法技术领域，具体涉及一种地下水库增湿作用下土
‑
结构相互作用计算方法。


背景技术：

2.为建设绿色“海绵城市”和解决城市雨季内涝问题，修建城市地下调蓄水库成为了解决办法之一。土体渗流
‑
应力耦合的数值模拟技术在地下工程、基坑等学科研究中发挥着不可或缺的作用。通过建立渗流
‑
应力耦合数值计算平台，地下结构稳定性分析过程中存在的复杂三维渗流、应力场问题得到科学分析与合理预测。此项技术在实际工程应用中保证了工程结构后期使用的安全和可靠。
3.我国的黄土土质分布广泛，黄土的湿陷性是不同于其他土质条件，遇水会导致结构强度发生改变。这种特性在工程实际中会引发一系列问题，比如建筑结构的不均匀沉降导致结构开裂等。在湿陷性黄土地区修建地下调蓄水库更易遇到此类工程问题。而传统工程模拟中并未能考虑在土、结构互相作用在非饱和渗流水分场、力学场应力耦合的力学响应情况，这与实际现象相去甚远。在实际工程实践中并未能起到真实模拟、实际防范、靶向指导的作用。
4.本发明针现有数值模拟存在的对以上缺陷，在comsol平台不仅实现了土
‑
结构在增湿前的相互作用下的力学响应即地下水库基底未漏水阶段土体、结构相互作用的力学响应而且实现了增湿作用即非饱和土渗透作用下，水分场、力学场耦合对土
‑
结构相互作用的影响即地下水库发生基底漏水阶段土体、结构相互作用对最关注的结构的变形和其最不利状态的力学响应。更合理地模拟了非饱和土发生渗流时水分场和力学场的耦合状态，更符合真实力学变化地模拟了增湿作用下土
‑
结构相互作用结构和土的力学响应。因此通过此软件仿真来实现地下水库未渗漏阶段、渗漏阶段的土
‑
结构相互作用的数值计算，实现地下调蓄水库底板渗漏情况下最合理的地下水库结构和土体的力学响应，为地下水库的结构设计和使用维护提供理论指导和技术支持。


技术实现要素：

5.本发明的目的是提供一种地下水库增湿作用下土
‑
结构相互作用计算方法，该方法能够实现地下水库未渗漏阶段、渗漏阶段的数值计算，为黄土湿陷性地区的非饱和渗流工程问题提供模拟预估，进一步为黄土地区的地下调蓄水库结构设计优化提供理论指导和技术支持。
6.本发明所采用的技术方案是，一种地下水库增湿作用下土
‑
结构相互作用计算方法，具体按照以下步骤实施：
7.步骤1：构建调蓄水库混凝土几何结构、土体几何结构模型；
8.步骤2：设置土体、水库混凝土材料各参数；
9.步骤3：定义“广义拉伸”、压力函数变量、土的强度指标和模量关于有效饱和度函
数变量；
10.步骤4：设置“固体力学”并设置各边界条件；增加“理查兹方程”并划分网格；
11.步骤5：设置水体“流体属性”“基本属性”各参数并划分网格；
12.步骤6：添加“研究”并计算“研究1”；
13.步骤7：对水库边界再次进行“布尔操作和分割”；
14.步骤8：设定“理查兹方程”参数；
15.步骤9：基于“研究1”，使用“辅助扫描”添加水库注水功能；添加并设定“研究2”时间单位、时间步；
16.步骤10：添加并设定“固体力学”并设置各参数条件并划分网格；
17.步骤11：基于“研究2”，添加并设定“研究3”；
18.步骤12：结果计算后数据成图。
19.本发明的特征还在于，
20.步骤1具体按照以下步骤实施：
21.步骤1.1：打开comsol软件，在“新建”点击“模型向导”，接下来在“选择空间维度”点击“二维”；
22.步骤1.2：接上一步，进入到“选择物理场”界面，选择“结构力学”模块，在该模块树中选择并添加物理场接口——“固体力学(solid)”；
23.步骤1.3：接上一步，接着点击“研究”，在“选择研究”点击“一般研究”，并在该板块树中选择“稳态”；后点击“完成”；
24.步骤1.4：接上一步，为构建水库、土体几何模型进入下一界面；在模型开发器中选择“几何”，右键“几何”选择“矩形”，该矩形a作为土体结构；在“几何设置”的“大小和形状”中分别设置“宽度”为a，“高度”为b，在“设置”窗口点击“构建选定对象”；再绘制矩形b(水库结构的框架)，右键“几何”选择“矩形”在“几何设置”的“大小和形状”中分别设置“宽度”为c，“高度”为d(上述的a，b为土体模拟尺寸，c，d为水库框架模拟尺寸)；在“图形”窗口栏选择“编辑”将矩形b移动到矩形a中部，在菜单栏选择“几何”中“布尔操作和分割”下的“差集”；然后在“差集”中将“要添加的对象”选为矩形a并激活，将“要减去的对象”选为矩形b并激活(其余选择为默认状态)，点击“构建选定对象”；再回到“几何”右键“矩形”绘制“宽度”为c，“高度”为d的矩形b2和b大小一致，点击绘制；在“图形”栏选择“编辑”将b2一定到矩形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。
25.步骤2具体按照以下步骤实施：
26.步骤2.1：在模型开发器中，右键“材料”选择“空材料”修改为“土”，将“几何实体层”“域”设定为“1”并激活；再从材料属性中找到“粘聚力”、“内摩擦角”增加到该“材料属性明细”，依次输入土的“杨氏模量”、“泊松比”、“密度”、“粘聚力”、“内摩擦角”；右键“材料”选
择“空材料”修改为“混凝土”，将“几何实体层”“域”设定为“2”并激活；再从材料属性中找到“单轴抗拉强度”“单轴抗压强度”“双轴抗压强度”增加到该“材料属性明细”，依次输入混凝土的“杨氏模量”、“泊松比”、“密度”、“单轴抗拉强度”、“单轴抗压强度”、“双轴抗压强度”(根据实际数据输入参数)；
27.步骤2.2：在模型开发器中，“几何”树的最后的“形成装配体”并选择全部区域。
28.步骤3具体按照以下步骤实施：
29.步骤3.1：在模型开发器中，右键单击“定义”在“组件耦合”下选择“广义拉伸”，在“图形”界面“几何实体层”选择“边界”，将矩形d结构的上下边界作为“源选择”，并将“目标映射”“x，y”均修改为大写；将“源”中“源坐标系”修改为“材料(x,y,z)”；再次右键单击“定义”在“组件耦合”下选择“广义拉伸”，在“图形”界面“几何实体层”选择“边界”，将矩形d结构的左右边界作为“源选择”，并将“目标映射”“x，y”均修改为大写；将“源”中“源坐标系”修改为“材料(x,y,z)”；
30.步骤3.2：在“模型开发器”下“组件”“定义”右键选择“变量”，设置侧边和底边的水压力函数，即侧边函数为随着水库高度变化的压力函数——load1＝h*9.8*(h1<y<h2)，并设定水库底部随水位高度变化的竖直向压力，load2＝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
‑
h1)*(h1<y<h2)”和water_h＝0(上式中h1是水库底边注水起始位置，h2是水库注水上限，water_h为水库储水高度)。
31.步骤4具体按照以下步骤实施：
32.步骤4.1：“模型开发器”下“固体力学”，“线弹性材料”右键选择“土壤塑性”、“混凝土”；其中，“土壤塑性”选择“域1”部分，“混凝土”选择“域2”部分，“屈服准则”采用“drucker
‑
prager”并选择“匹配莫尔
‑
库伦准则”，内聚力和内摩擦角均选自“来自材料”；
33.步骤4.2：在菜单栏“物理场”“域”中选择“重力”其作用域为全部域；“边界”中选择“固定约束”作用边界为域1的下边界；并选择“对称”作用边界为域2的左右边界；在“边界”中选择“指定位移”，依次作为x、y方向的指定位移；添加第一组“指定位移”，边界选择域2的左右外边界，在“设置”中“指定位移”勾选“在x方向指定”并写入前面“定义”中的广义拉伸的算子名称“genext1(x)”；添加第二组“指定位移”边界选择域2的上下外边界，在“设置”中“指定位移”勾选“在y方向指定”并写入前面“定义”中的广义拉伸的算子名称“genext1(y)”；
34.步骤4.3：在菜单栏选择“物理场”中的“边界荷载”，添加第一组“边界荷载”选择水库内部的沿高度方向的内壁，并在“力”中添加前面“定义”中的“load
1”，添加第二组“边界荷载”选择水库内部的底部内壁，并在“力”中添加前面“定义”中的“load
2”；
35.步骤4.4：在菜单栏中选择“物理场”“添加物理场”，在“流体流动”接口选择“多孔介质和地下水流”中的“理查兹方程”。
36.步骤5具体按照以下步骤实施：
37.步骤5.1：在“理查兹方程”下，“设置”中，选定域1、域2；在“流体属性”下，将密度设
为“1000kg/m
3”；在“基本属性”设置“多孔材料”为“土体”材料，“饱和液体体积分数”、残余液体体积分数，根据该地区地质勘察报告获取；“渗透率模型”选定“水力传导率”(各向同性)；根据地质勘察报告获；并在“储水模型”中选择“线性储水”，设定“流体的可压缩性”为“4.4e
‑
10[1/pa]”，“基体有效压缩率”为“10e
‑
6[1/pa]”；“滞留模型”选定“van genuchten”，“本构关系常数”的“α、n、l”根据勘察报告获取；
[0038]
步骤5.2：定位到“模型开发器”的“初始值”；定位到“初始值”的“设置”，选择全部区域；选定“压力头”设定“h
p
＝
‑
0.1m”(该数值可以根据实际调整，
‑
0.1m为举例)；定位到“模型开发器”下的“理查兹方程”右键选择“压力头”。
[0039]
步骤6具体按照以下步骤实施：
[0040]
基于上述“压力头”设定；回到“几何”设置；计算“固体力学”下的“研究1”；定位到菜单栏的“结果”，选定“二维绘图组”，右键单击选择“表面”，定位到“设置”，选定“数据”的“数据集”为“固体力学”，点击“模型开发器”的“表面”，定位到“设置”中的“表达式”，输入“solid.pm”，点击“绘制”；得到立柱在受周围土体重力、自身重力的弹性受力下的应力云图，获得立柱底部受力分布图，大概量测得到立柱底端受力从坡段到平缓过渡的水平投影长度，将该长度作为前“理查兹方程”立柱渗流区域的入渗端；
[0041]
步骤7具体按照以下步骤实施：
[0042]
再定位回“模型开发器”的“几何”部分，在菜单栏选择“布尔操作和分割”下的“分割边”；回到图形区域，基于前面选定的渗流投影区长度，分割出对应长度的渗流端，分别获得各立柱下的渗流边(例如立柱下渗流区、水库两角点处渗流区)；
[0043]
步骤8具体按照以下步骤实施：
[0044]
定位回“理查兹方程”，右键单击“模型开发器”下的“理查兹方程”选择“压力头”，定位到“设置”部分的“边界选择”选择上一步设置好的渗流边，在“压力头”处输入“h
p0
＝0.2m”(0.2m为举例，可根据实测压力头输入该值)。
[0045]
步骤9具体按照以下步骤实施：
[0046]
步骤9.1：定位回“固体力学”，在对应计算“研究1”；定位到“设置”，在“物理场和变量选择”只勾选“固体力学”；勾选“研究扩展”下的“辅助扫描”；选定“water_h”设定“参数值列表”即“range(h1,1,h2)”(h1是水库底边注水起始位置，h2是水库注水上限)；单击“计算”；
[0047]
步骤9.2：定位到菜单栏：单击“研究”选择“添加研究”，选定“一般研究”下的“瞬态”；定位到“设置”，选择“时间单位”为“d”，设置“时间步”为“range(0,3,90)”(本例按照三天为一个计算周期，计算三个月的渗流场)；在“物理场和变量选择”只勾选“理查兹方程”；在“因变量值”下的“设置”选择“用户控制”，“方法”选择“解”，“研究”选择“研究1”，“参数值water_h”选择“自动”，其余默认；单击“计算”。
[0048]
步骤10具体按照以下步骤实施：
[0049]
步骤10.1：定位回菜单栏的“物理场”，选择“添加物理场”，在“结构力学”下选择“固体力学”；添加该物理场，并重命名为“固体力学(二)”；定位到“模型开发器”的“固体力学(二)”，右键“线弹性材料”，选定“土壤塑性”、“混凝土”、“外部应力”；定位到“线弹性材料”，将“域1、2”选中作为域选择，“杨氏模量”选用“用户定义”输入“e”，其余默认；定位到“土壤塑性”将区域选择为“域1”，“屈服准则”选定为“drucker
‑
prager准则”，并勾选“匹配莫尔
‑
库伦准则”，“内聚力”选定“用户定义”填入“c”，“内摩擦角”选用“用户定义”填入
其余默认；定位到“混凝土”选定“域2”，“混凝土准则”选用“bresler
‑
pister准则”，“单轴抗拉强度”、“单轴抗压强度”、“双轴抗压强度”均选择“来自材料”(根据实际材料设定)；定位到“外部应力”，选择“域1”，“应力输入”选定“孔隙压力”，“绝对压力”选定“用户定义”并输入“p*(p>0)”，其余默认；
[0050]
步骤10.2：定位回“模型开发器”，右键“固体力学(二)”，选定“指定位移”，并同样设定“对称”、“重力”，同前述“固体力学”一样设置步骤。
[0051]
步骤11具体按照以下步骤实施：
[0052]
定位回菜单栏，点击“研究”并“添加研究”，选定“一般研究”下的“稳态”；定位到“设置”的“物理场和变量选择”，勾选“修改研究步骤的模型配置”，只禁用“固体力学”；定位到“因变量值”在“求解变量的初始值”下“设置”选定“用户控制”，“方法”选定“解”，研究选定“研究2”(即理查兹物理场对应的研究)，“解”选择“解2”，“时间(d)”选择“30d”(举例为30d)；点击“计算”；
[0053]
以上计算基于相同网格划分，如图2所示。
[0054]
步骤12具体按照以下步骤实施：
[0055]
按照上述方案，后处理绘制土体的压力、杨氏模量、内摩擦角、杨氏模量、粘聚力的表面分布图；结构的压力表面图，结构的位移变化曲线。
[0056]
本发明的有益效果是：本发明在comsol软件平台，建立简化的二维平面内土
‑
结构联结物理模型，实现水库不同储水高度下结构的力学响应；建立土体中采用粘聚力、杨氏模量、内摩擦角分别关于有效饱和度的函数作为固体力学的相应参数即考虑增湿作用下土
‑
结构相互作用计算，对结构在“破坏
‑
安全情况”即水库基底破坏发生渗漏的情况下进行渗流模拟，实现了结构破坏、土体渗流后水分场和力学场的耦合作用下结构构件的力学响应。该发明提出的该种计算方法实现地下水库未渗漏阶段、渗漏阶段的数值计算，为黄土湿陷性地区的非饱和渗流工程问题提供模拟预估，进一步为黄土地区的地下调蓄水库结构设计优化提供理论指导和技术支持。
附图说明
[0057]
图1是本发明模型构建的水库
‑
土体几何模型图；
[0058]
图2是本模型计算前土
‑
结构相互联结网格划分图；
[0059]
图3是本模型经计算得到水库结构受上覆土体荷载和注水情况压力表面分布图；
[0060]
图4是本发明实施案例中土体模拟渗流发生30天后土的有效饱和度随渗流时间发生改变的表面分布云图；
[0061]
图5是本发明实施案例中土体模拟渗流发生60天后土的有效饱和度随渗流时间发生改变的表面分布云图；
[0062]
图6是本发明实施案例中土体模拟渗流发生90天后土的有效饱和度随渗流时间发生改变的表面分布云图；
[0063]
图7是本发明实施案例中土体模拟渗流发生30天后土的粘聚力随渗流时间发生改变的表面分布云图；
[0064]
图8是本发明实施案例中土体模拟渗流发生60天后土的粘聚力随渗流时间发生改变的表面分布云图；
[0065]
图9是本发明实施案例中土体模拟渗流发生90天后土的粘聚力随渗流时间发生改变的表面分布云图；
[0066]
图10是本发明实施案例中土体模拟渗流发生30天后土的杨氏模量随渗流时间发生改变的表面分布云图；
[0067]
图11是本发明实施案例中土体模拟渗流发生60天后土的杨氏模量随渗流时间发生改变的表面分布云图；
[0068]
图12是本发明实施案例中土体模拟渗流发生90天后土的杨氏模量随渗流时间发生改变的表面分布云图；
[0069]
图13是本发明实施案例中土体模拟渗流发生30天后土的内摩擦角随渗流时间发生改变的表面分布云图；
[0070]
图14是本发明实施案例中土体模拟渗流发生60天后土的内摩擦角随渗流时间发生改变的表面分布云图；
[0071]
图15是本发明实施案例中土体模拟渗流发生90天后土的内摩擦角随渗流时间发生改变的表面分布云图；
[0072]
图16是本发明实施案例中水库结构发生渗漏30天得到的压力表面云图；
[0073]
图17是本发明实施案例中水库结构发生渗漏30天得到的结构四周位移量曲线分布结构图；
[0074]
图18是本发明实施案例中水库结构发生渗漏30天得到的结构四周压力分布结构图。
具体实施方式
[0075]
下面结合附图说明和具体实施方式对本发明进行详细说明。
[0076]
本发明提供一种地下水库增湿作用下土
‑
结构相互作用计算方法，
[0077]
步骤1：构建调蓄水库混凝土几何结构、土体几何结构模型。(如图1所示，内外依次划分为域2和域1，内部均匀分布三根支柱)
[0078]
步骤1.1：打开comsol软件，在“新建”点击“模型向导”并“选择空间维度”点击“二维”；
[0079]
步骤1.2:接上一步，进入到“选择物理场”界面，选择“结构力学”模块，在该模块树中选择并添加物理场接口——“固体力学(solid)”；
[0080]
步骤1.3：接上一步，接着点击“研究”，在“选择研究”点击“一般研究”，并在该板块树中选择“稳态”；后点击“完成”。
[0081]
步骤1.4：接上一步，为构建水库、土体几何模型进入下一界面。在模型开发器中选择“几何”，右键“几何”选择“矩形”，该矩形a作为土体结构。在“几何设置”的“大小和形状”中分别设置“宽度”为40.4m，“高度”为24m，在“设置”窗口点击“构建选定对象”。再绘制矩形b(水库结构的框架)，右键“几何”选择“矩形”在“几何设置”的“大小和形状”中分别设置“宽度”为16m，“高度”为12m。在“图形”窗口栏选择“编辑”将矩形b移动到矩形a中部，在菜单栏选择“几何”中“布尔操作和分割”下的“差集”。然后在“差集”中将“要添加的对象”选为矩形a并激活，将“要减去的对象”选为矩形b并激活(其余选择为默认状态)，点击“构建选定对象”。再回到“几何”右键“矩形”绘制“宽度”为12m，“高度”为16m的矩形b2和b大小一致，点击
绘制。在“图形”栏选择“编辑”将b2一定到矩形a的空白位置处(其余选项默认)。采取同上述步骤相同的原理，再绘制宽14.4m，高10.4m的矩形c，然后将矩形c移动到矩形b的中部，再使用“差集”，用矩形b减去矩形c(取消“保留内部对象”)，并点击“构建选定对象”得到水库内无立柱框架结构d。添加立柱，仍旧右键“几何”选择“矩形”，依次设置若根个宽为0.8m，高度为10.4m的立柱e，并再增设同样规格的立柱两根f和g。均匀分布这些立柱在矩形d结构内部，再利用“布尔操作和分割”下的“并集”，将水库框架d和立柱e、f、g取并集(取消“保留内部边界”)，点击“构建选定对象”。以上得到周围土体——域1、水库结构——域2。
[0082]
步骤2：设置土体、水库混凝土材料各参数
[0083]
步骤2.1：在模型开发器中，右键“材料”选择“空材料”修改为“土”，将“几何实体层”“域”设定为“1”并激活。再从材料属性中找到“粘聚力”、“内摩擦角”增加到该“材料属性明细”，依次输入土的“杨氏模量”取20e6pa、“泊松比”取0.3、“密度”取1800kg/m3，“粘聚力”取35e3pa、“内摩擦角”取26.4[deg]rad；右键“材料”选择“空材料”修改为“混凝土”，将“几何实体层”“域”设定为“2”并激活。再从材料属性中找到混凝土的“杨氏模量”、“泊松比”、“密度”、“单轴抗拉强度”、“单轴抗压强度”、“双轴抗压强度”增加到该“材料属性明细”。依次输入“单轴抗拉强度”取1.39e6pa，“单轴抗压强度”取13.8e6pa,“双轴抗压强度”取21e6pa。
[0084]
步骤2.2：在模型开发器中，“几何”树的最后的“形成装配体”并选择全部区域
[0085]
步骤3：定义“广义拉伸”、压力函数变量、土强度指标和模量关于有效饱和度函数变量
[0086]
步骤3.1：在模型开发器中，右键单击“定义”在“组件耦合”下选择“广义拉伸”，在“图形”界面“几何实体层”选择“边界”，将矩形d结构的上下边界作为“源选择”，并将“目标映射”“x，y”均修改为大写。将“源”中“源坐标系”修改为“材料(x,y,z)”；再次右键单击“定义”在“组件耦合”下选择“广义拉伸”，在“图形”界面“几何实体层”选择“边界”，将矩形d结构的左右边界作为“源选择”，并将“目标映射”“x，y”均修改为大写。将“源”中“源坐标系”修改为“材料(x,y,z)”；
[0087]
步骤3.2：在“模型开发器”下“组件”“定义”右键选择“变量”，设置侧边和底边的水压力函数，即侧边函数为随着水库高度变化的压力函数——load1＝h*9.8*(2.8<y<h2)，并设定水库底部随水位高度变化的竖直向压力，load2＝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”将该函数名称定义为eint；在“组件”的“定义”下的“变量”中写入“e＝eint(dl.se)”；在“全局定义”右键“变量”输入“h＝(y
‑
2.8)
*
(2.8<y<12.8)”。
[0088]
步骤4：设置“固体力学”并设置各边界条件；增加“理查兹方程”并划分网格。如图2，水库结构和土体均采用“较细化”网格划分；
[0089]
步骤4.1：“模型开发器”下“固体力学”，“线弹性材料”右键选择“土壤塑性”、“混凝土”。其中，“土壤塑性”选择“域1”部分，“混凝土”选择“域2”部分，“屈服准则”采用“drucker
‑
prager”并选择“匹配莫尔
‑
库伦准则”，内聚力和内摩擦角均选自“来自材料”。
[0090]
步骤4.2：在菜单栏“物理场”“域”中选择“重力”其作用域为全部域；“边界”中选择“固定约束”作用边界为域1的下边界；并选择“对称”作用边界为域2的左右边界；在“边界”中选择“指定位移”，依次作为x、y方向的指定位移。添加第一组“指定位移”，边界选择域2的左右外边界，在“设置”中“指定位移”勾选“在x方向指定”并写入前面“定义”中的广义拉伸的算子名称“genext1(x)”；添加第二组“指定位移”边界选择域2的上下外边界，在“设置”中“指定位移”勾选“在y方向指定”并写入前面“定义”中的广义拉伸的算子名称“genext1(y)”。
[0091]
步骤4.3：在菜单栏选择“物理场”中的“边界荷载”，添加第一组“边界荷载”选择水库内部的沿高度方向的内壁，并在“力”中添加前面“定义”中的“load
1”，添加第二组“边界荷载”选择水库内部的底部内壁，并在“力”中添加前面“定义”中的“load
2”。
[0092]
步骤4.4：在菜单栏中选择“物理场”“添加物理场”，在“流体流动”接口选择“多孔介质和地下水流”中的“理查兹方程”。
[0093]
步骤5：设定水体“流体属性”“基本属性”各参数
[0094]
步骤5.1：在“理查兹方程”下，“设置”中，选定域1、域2。在“流体属性”下，将密度设为“1000kg/m
3”；在“基本属性”设置“多孔材料”为“土体”材料，“饱和液体体积分数”取0.497、“残余液体体积分数”取0.144。“渗透率模型”选定“水力传导率”(各向同性)，水力传导率取1.5e
‑
7m/s。并在“储水模型”中选择“线性储水”，设定“流体的可压缩性”为“4.4e
‑
10[1/pa]”，“基体有效压缩率”为“10e
‑
6[1/pa]”。“滞留模型”选定“van genuchten”，“本构关系常数”的“α、n、l”分别为取22.21[1/m]、1.127、0.1126。
[0095]
步骤5.2：定位到“模型开发器”的“初始值”。定位到“初始值”的“设置”，选择全部区域。选定“压力头”设定“h
p
＝
‑
0.1m”。定位到“模型开发器”下的“理查兹方程”右键选择“压力头”。
[0096]
步骤6：添加“研究”并计算“研究1”[0097]
步骤6.1：基于上述“压力头”设定。回到“几何”设置。计算“固体力学”下的“研究1”。定位到菜单栏的“结果”，选定“二维绘图组”，右键单击选择“表面”，定位到“设置”，选定“数据”的“数据集”为“固体力学”，点击“模型开发器”的“表面”，定位到“设置”中的“表达式”，输入“solid.pm”，点击“绘制”。得到立柱在受周围土体重力、自身重力的弹性受力下的应力云图，获得立柱底部受力分布图，大概量测得到立柱底端受力从坡段到平缓过渡的水平投影长度，将该长度作为前“理查兹方程”立柱渗流区域的入渗端。
[0098]
步骤:7：对水库边界再次进行“布尔操作和分割”[0099]
步骤7.1：再定位回“模型开发器”的“几何”部分，在菜单栏选择“布尔操作和分割”下的“分割边”。回到图形区域，基于前面选定的渗流投影区长度，分割出对应长度的渗流端，分别获得各立柱下的渗流边(例如立柱下渗流区、水库两角点处渗流区)。
[0100]
步骤8：设定“理查兹方程”参数
[0101]
步骤8.1：定位回“理查兹方程”，右键单击“模型开发器”下的“理查兹方程”选择“压力头”，定位到“设置”部分的“边界选择”选择上一步设置好的渗流边，在“压力头”处输入“h
p0
＝0.2m”。
[0102]
步骤9：基于“研究1”，使用“辅助扫描”添加水库注水功能；添加并设定“研究2”时间单位、时间步
[0103]
步骤9.1：定位回“固体力学”，在对应计算“研究1”。定位到“设置”，在“物理场和变量选择”只勾选“固体力学”。勾选“研究扩展”下的“辅助扫描”。选定“water_h”设定“参数值列表”即“range(2.8,1,12.8)”(2.8m是水库底边注水起始位置，12.8m是水库注水上限)。单击“计算”。如图3所示，水库上下侧结构的压力在内部立柱位置左右出现最大值，立柱体沿纵向出现剪切应力。
[0104]
步骤9.2：定位到菜单栏。单击“研究”选择“添加研究”，选定“一般研究”下的“瞬态”，记为“研究2”。定位到“设置”，选择“时间单位”为“d”，设置“时间步”为“range(0,3,90)”。在“物理场和变量选择”只勾选“理查兹方程”。在“因变量值”下的“设置”选择“用户控制”，“方法”选择“解”，“研究”选择“研究1”，“参数值water_h”选择“自动”，其余默认。单击“计算”。
[0105]
如图4—图15所示，分别表示水库发生渗漏30天、60天、90天，水库底部的土的有效饱和度、粘聚力、杨氏模量、内摩擦角分布状态云图。明显发现，随着渗流时间的持续，土体的有效饱和度逐渐增大，而粘聚力、杨氏模量、内摩擦角随渗流持续时间逐渐减小；渗流发生的同一个时间，随着土层高程从大到小变化，有效饱和度从大到小变化而粘聚力、杨氏模量、内摩擦角反而由小到大变化。
[0106]
步骤10：添加并设定“固体力学”并设置各参数条件并划分网格
[0107]
步骤10.1：定位回菜单栏的“物理场”，选择“添加物理场”，在“结构力学”下选择“固体力学”。添加该物理场，并重命名为“固体力学(二)”。定位到“模型开发器”的“固体力学(二)”，右键“线弹性材料”，选定“土壤塑性”、“混凝土”、“外部应力”。定位到“线弹性材料”，将“域1、2”选中作为域选择，“杨氏模量”选用“用户定义”输入“e”，其余默认。定位到“土壤塑性”将区域选择为“域1”，“屈服准则”选定为“drucker
‑
prager准则”，并勾选“匹配莫尔
‑
库伦准则”，“内聚力”选定“用户定义”填入“c”，“内摩擦角”选用“用户定义”填入其余默认。定位到“混凝土”选定“域2”，“混凝土准则”选用“bresler
‑
pister准则”，“单轴抗拉强度”、“单轴抗压强度”、“双轴抗压强度”均选择“来自材料”。定位到“外部应力”，选择“域1”，“应力输入”选定“孔隙压力”，“绝对压力”选定“用户定义”并输入“p
*
(p>0)”，其余默认。
[0108]
步骤10.2：定位回“模型开发器”，右键“固体力学(二)”，选定“指定位移”，并同样设定“对称”、“重力”，同前述“固体力学”一样设置步骤。
[0109]
步骤11：基于“研究2”，添加并设定“研究3”：
[0110]
步骤11.1：定位回菜单栏，点击“研究”并“添加研究”，选定“一般研究”下的“稳态”。定位到“设置”的“物理场和变量选择”，勾选“修改研究步骤的模型配置”，只禁用“固体力学”。定位到“因变量值”在“求解变量的初始值”下“设置”选定“用户控制”，“方法”选定“解”，研究选定“研究2”(即理查兹物理场对应的研究)，“解”选择“解2”，“时间(d)”选择“30d”(举例为30d)。点击“计算”。
[0111]
如图16所示，表示水库发生渗漏后，增湿作用下土
‑
结构相互作用发生物理场耦合的结构压力云图即水库上下侧和立柱所受压力较未渗流时有明显变化。
[0112]
步骤12：结果计算后数据成图
[0113]
步骤12.1：按照上述方案，图3—图16分别示意后处理绘制土体的压力、杨氏模量、内摩擦角、杨氏模量、粘聚力的表面分布图；图17和图18示意结构的压力表面图，结构的位
移变化曲线，直观体现水库发生渗流后水库结构各边界位移大小以及压力大小。
[0114]
图17所示，水库上下侧结构被土体向内挤入，左右两侧结构发生反弯点上部向外扩出，下侧向内挤入的情况，整体位移变形左右对称；图18所示，水库上下两侧结构受到的压力大小按内部立柱位置变化依次大小交错，左右两侧压力上部小下部大。