目录
一、模型向导
二、建模
1. 全局参数
2. 定义
（1）插值
（2）变量
（3）状态变量
（4）分段函数
3. 几何
4. 材料
5. 物理场
（1）固体力学
（2）达西定律
（3）多孔介质传热
6. 网格
三、计算
1. 求解器设置
2. 计算
损伤及属性效果：
一、模型向导
空间：【二维】
物理场：【多孔弹性 实体】+【多孔介质传热】
研究：【稳态】，向导完成后再添加一个【瞬态】
二、建模
1. 全局参数
| 名称 | 表达式 | 值 | 描述 |
|---|---|---|---|
| fc | 160[MPa] | 1.6E8 Pa | 单轴抗压强度 |
| phi | 40[deg] | 0.69813 rad | 内摩擦角 |
| poro_0 | 0.01 | 0.01 | 初始孔隙率 |
| poro_r | 1.00E-03 | 0.001 | 残余孔隙率 |
| mu | 0.2 | 0.2 | 泊松比 |
| k0 | 1e-19[m^2] | 1E-19 m² | 初始渗透率 |
| pinj | 30[MPa] | 3E7 Pa | 注入压力定值 |
| landa0 | 4[W/m/K] | 4 W/(m·K) | 初始固体热传导系数 |
| Cs0 | 950[J/kg/K] | 950 J/(kg·K) | 固体比热容 |
| h_tc0 | 1000[W/(m^2*K)] | 1000 W/(m²·K) | 初始传热系数 |
| Tinit | 373.15[K] | 373.15 K | 初始温度 |
| Tinj | 293.15[K] | 293.15 K | 注入温度 |
| h_expan_s | 6e-6[1/K] | 6E-6 1/K | 固体热膨胀系数 |
| h_expan_l | 1e-3[1/K] | 0.001 1/K | 用于算c2 |
| h_expan_m | 6.2e-6[1/K] | 6.2E-6 1/K | 用于算c2 |
2. 定义
（1）插值
【定义】→【函数】→【插值】
创建两个插值函数，一个是初始杨氏模量 E，另一个是初始抗拉强度 ft
【数据源】：文件
选择 MATLAB 随机生成的 E 和 ft 数据表后导入
MATLAB 参考代码：（适用于2019以下版本）
value = 16; %抗拉强度 ft 参考数值 %value = 35; %杨氏模量 E 参考数值 hang = 0.3; %m %宽度 lie = 0.3; %m %高度 n = 0.002; %取样间隔 pro = wblrnd(value,20,hang/n,lie/n); mat = []; A = linspace(0,hang,hang/n); %矩阵转换 mat = ones(hang/n,1)*A; B = mat(:); C = mat'; C = C(:); D = pro'; D = D(:); T = [B C D]; min(D) max(D) mean(D) xlswrite('weibull-ft.xlsx',T); %xlswrite('weibull-E.xlsx',T); [x1,y1]=meshgrid(A,A); surf(x1,y1,pro) shading interp杨氏模量 E 的单位为 GPa
抗拉强度 ft 单位为 MPa
（2）变量
| 名称 | 表达式 |
|---|---|
| F1 | solid.sp1-ft0 |
| F2 | -solid.sp3+(1+sin(phi))/(1-sin(phi))*solid.sp1-fc |
| E0 | int1(x,y) |
| et0 | ft0/E0*0.5 |
| ec0 | fc/E0 |
| poro | poro_r+(poro_0-poro_r)*exp(alpha*p) |
| alpha | 3*(1-2*mu)/E0/(poro_0-poro_r) |
| km | k0*(poro/poro_0)^3*exp(5*D) |
| Dt | if((F1>=0)&&(d(F1,TIME)>0),min(0.99,max(0,1-(et0/solid.ep1)^2)),0) |
| Dc | -if((F2>=0)&&(d(F2,TIME)>0),min(0.99,max(0,1-(ec0/solid.ep3)^2)),0) |
| DD | min(0.99,Dt+abs(Dc)) |
| ft0 | int2(x,y) |
| landa | landa0*exp(D*5) |
| Cs | Cs0 |
| h_tc | if(abs(D)>0.95,5*h_tc0,h_tc0) |
| biot | if(D<0.9,0.5,1) |
| c1 | biot |
| c2 | poro*h_expan_l+(1-poro)*h_expan_s-h_expan_m*(1-biot) |
| Qs | solid.K*h_expan_s*Td*d(solid.eth11+solid.eth22,TIME) |
| Td | solid.Tdiff |
（3）状态变量
先在【显示更多选项】中勾选上【方程贡献】
（comsol 6.4版本是勾选这个，其他版本可能不是，在解释里找）
【定义】→【方程贡献】→【状态变量】
表达式：if(D<=DD,DD,D)
（4）分段函数
【定义】→【函数】→【分段】
3. 几何
【正方形】：边长 0.3
【圆】：半径 0.015，xy 0.15
【布尔】→【差集】，正方形减去圆形
【形成联合体】
4. 材料
（1）空材料
缺失属性会在设置完物理场的边界条件后自动显示
可以先创建好边界条件后再来设置具体的属性值
（2）water
添加【比热率】属性，这里设置为 1.0066
5. 物理场
（1）固体力学
【辊支承】：左 + 下边界，约束法向位移
【边界荷载】：上 + 右 + 内部圆（3个分开）
竖向荷载小
右边界
内部边界
【线弹性材料】→【热膨胀】
体积参考温度 Tref 设置为 Tinit，温度用 ht
（2）达西定律
【流体】
mat2 指的是材料中的水
【初始值】
【压力】内部边界设置，压力设置为 0.1 MPa
【入口】内部边界设置
压力
【质量源】（耦合项）
（3）多孔介质传热
初始温度用 Tinit
【流体】
【基体】
【温度】四周边界设置
【初始值】和【温度】中的 T 设置为 Tinit
【热通量】内部边界设置
【热源】（耦合项）
6. 网格
【自由三角形】→【大小 1】：整个域
单元最大 0.003，最小 0.0004，最大增长率 1.05
【大小 2】：内部边界
单元最大、最小：0.001
【细分方法】：Delaunay，会根据几何形状和网格大小要求，自动把计算区域划分成质量较好的三角形网格，尽量避免生成又细又扁的三角形
三、计算
1. 求解器设置
【稳态】求解器：
禁用【固体力学】的【压力】、【达西定律】的【入口】、【多孔介质传热】的【热通量】
【瞬态】求解器：
设置【输出时步】、【相对容差】
禁用【达西定律】中的【压力】
配置【因变量值】
打开【显示默认求解器】
【时间步进】设置
【求解时显示】设置
【全耦合】设置
2. 计算
计算【稳态】
【应力】改为【损伤】，修改表达式
禁用【变形】，删了也行
【温度】的单位改成摄氏度 degC
【损伤】、【压力】、【温度】的数据集改为研究2
计算【瞬态】
计算过程中会展示损伤变化
【损伤】计算结果
【压力】
【温度】
【孔隙率】
