产品
人工地下冻结(AGF, Artificial Ground Freezing)最早出现于俄国，用于金矿开采，后由德国工程师用于煤矿矿井建设并获得专利，技术渐趋成熟，现在已广泛应用于地铁、深基坑、矿井建设等工程中。人工地下冻结法的原理是将人工冻结管埋置于土体内，利用冻结管内循环的冷媒剂将土体中的热量逐渐带走从而形成强度、高密封性的冻结土壁，形成的冻结土壁通常构成封闭的环形，能够起到承受荷载和密封防水的作用。AGF法的适应性强、安全 可靠、无污染，但和其他方法相比，一般造价较高。
AGF法进行土体冻结，可以适用于存在潜在地下水的地层，为保证地下开挖施工的安全，先开凿冻结孔，下放冻结管并将冻结管与负责冷量循环的冻结站相连，待冻结站工作后每根冻结管周围会逐渐形成冻结圆柱，相邻冻结圆柱半径逐渐增大到相互连接后就形成了闭合的冻结墙(frozen wall)，冻结墙外侧的潜在地下水被封闭，同时冻结土体的强度比原土高很多，都有利于保证后续冻结墙内向下开挖土体施工的安全。
FLAC3D软件作为有限差分软件，除广泛应用于岩土体的力学问题分析外，还可以用于冻结温度场的分析。本文后续演示了采用人工地下冻结法进行煤矿井筒开挖的分析代码，是我在硕士阶段编写的，虽然代码还有不少值得完善的地方。但里面的包含了诸如坐标数据的读入、自动截图和数据文件的自动取名保存等等。 希望对大家学习FLAC3D有帮助。
冻结方式为竖向冻结，冷媒是盐水。
代码：
; ----------------------------------------
; Description: 竖井温度场分析
; Author: Qiao Cheng
; ----------------------------------------
new
config thermal
; --------------------
; 定义几何模型参数
; --------------------
def Geocons ; geometrical constants
eps = 1.e-3
meps = -1 * eps
;z coordinate of top
z1 = 1.
z1m = z1 - eps
z1p = z1 + eps
;
out_b = 21. ;模型最外侧边界
mout_b = - out_b
out_b_m = out_b - eps ;for boundary condition define
out_b_p = out_b + eps
;
tempbc = 16.
tempbcm = tempbc - eps
;
c0rad = 12.
mc0rad = - c0rad
;
; 冻结管布置圈径(半径)/m
c1rad = 15.3/2.
mc1rad = - c1rad
c1radm = c1rad - eps
c1radp = c1rad + eps
;
outer_rad = 10.4/2. ; equalant to original 10.4
mouter_rad = -1*outer_rad
outer_radm = outer_rad - eps
outer_radp = outer_rad + eps
;
inner_rad = 10.4/2. -0.5 ; out side of inner wall
minner_rad = -1*inner_rad
inner_radm = inner_rad - eps
inner_radp = inner_rad + eps
;
excrad = 7.6/2.
mexcrad = (-1.0)*excrad
excradp = excrad + eps
;
_nz = 1 ;devided along z direction
_zgpn = _nz + 1 ;gridpoints along z direction
_zslen = z1 / _nz ;沿z轴一个节段的长度
;
P_times = 60
M_times = 40
;
contiflag = 1.
have_frozen = 0.
thawed = 0.
end
Geocons
; ----------------------------------------
; 冻结管根数及温度变量声明
def PipeCal
_qtrp1num = C1N / 2
_doup1num = 2 * C1N
_c1angle = 2 * pi / C1N
;
_yanshui = -10.
_initemp = 13.
end
;设置每圈的冻结孔数
set C1N 36
PipeCal
; ----------------------------------------
def setup_io
IO_READ = 0
IO_WRITE = 1
IO_FISH = 0
IO_ASCII = 1
filename1 = 'c1.dat'
positemp = 'p_temp.dat'
maintemp = 'm_temp.dat'
end
setup_io
; ----------------------------------------
; Read the brine pipe coordinates from file
def ReadCoords
;冻结孔坐标(考虑偏斜)
array _x1(C1N), _y1(C1N), _z1(5)
array eatit(1)
status = open(filename1, IO_READ, IO_ASCII)
if status = 6 then
oo = out('文件'+filename1+'不存在')
endif
status = read(_x1, C1N)
status = read(eatit,1)
status = read(_y1, C1N)
if status # C1N then
oo = out('数据文件行数不正确')
endif
status = close
end
set C1N 10
ReadCoords
; ----------------------------------------
; Read temperature of positive freezing stage from file
def ReadPositiveTemp
array _positivetemp(P_times)
status = open(positemp, IO_READ, IO_ASCII)
if status = 6 then
oo = out('文件'+positemp+'不存在')
endif
status = read(_positivetemp, P_times)
if status # P_times then
oo = out('数据文件行数不正确')
endif
status = close
end
ReadPositiveTemp
; ----------------------------------------
def ReadMaintainTemp
array _mainttemp(M_times)
status = open(maintemp, IO_READ, IO_ASCII)
if status = 6 then
oo = out('文件'+maintemp+'不存在')
endif
status = read(_mainttemp, M_times)
if status # M_times then
oo = out('数据文件行数不正确')
endif
status = close
end
ReadMaintainTemp
; ----------------------------------------
def Calparas
c_k = 55.56e6 ; bulk modulus
c_g = 18.52e6 ; shear modulus
;
_cond = 1.45 ;导热系数，未冻结[W/(m*℃)]
_cv = 1.3e3 ;比热，未冻结 [J/kg/℃]
;
_cond0 = 1.55
_cv0 = 1.2e3
;
_cond1 = 1.6
_cv1 = 1.9e3
;
_cond2 = 1.65
_cv2 = 1.0e3
;
_cond3 = 1.7
_cv3 = 0.95e3
;
_cond4 = 1.75
_cv4 = 0.9e3
;
_cond5 = 1.8
_cv5 = 0.85e3
end
Calparas
; ----------------------------------------
def CalculateSteps
_positive_days = 150
_totaldays = 270
_1day_steps = int(24*3600/_dtval)
_positive_steps = _positive_days * _1day_steps
_one_positive_stage_steps = int(_positive_steps / P_times)
_maintainsteps = (_totaldays - 150)*_1day_steps
_one_maint_stage_steps = int(_maintainsteps / M_times)
_8day_age = 8*86400
end
set _dtval 360.
CalculateSteps
; print _1day_steps
; print _one_positive_stage_steps
; pause
; ----------------------------------------
; setThawProps
; ---------------------
def SetThawProps
p_z = zone_head
loop while p_z # null
if z_model(p_z) # 'null' then
if z_group(p_z) # 'model_board' then
if z_group(p_z) # 'inner_concrete' then
if z_group(p_z) # 'outer_concrete' then
if z_temp(p_z) > -10. then
z_prop(p_z,'conductivity') = _cond3
z_prop(p_z,'spec_heat') = _cv3 + 0.3e3
z_prop(p_z,'thexp') = 1.5e-5
if z_temp(p_z) > -2. then
z_prop(p_z,'conductivity') = _cond1
z_prop(p_z,'spec_heat') = _cv0 + 0.2e3
; _cv1 is too large, lead to too slowly thraw
z_prop(p_z,'thexp') = 1.e-5
if z_temp(p_z) > 0. then
z_group(p_z) = 'unfrozen'
thawed = 1.
z_prop(p_z,'conductivity') = _cond0
z_prop(p_z,'spec_heat') = _cv0 + 0.4e3
z_prop(p_z,'thexp') = 0.
endif
endif
endif
endif
endif
endif
endif
p_z = z_next(p_z)
end_loop
end
; --------------------
gen zone cyl p1 excrad 0 0 p2 0 0 z1 p3 0 mexcrad 0 size 10 _nz 54 rat 0.8 1 1 & group excavatedsoil gen zone cshell p1 inner_rad 0 0 p2 0 0 z1 p3 0 minner_rad 0 size 4 _nz 54 & rat 1 1 1 dim excrad excrad excrad excrad group inner_concrete range group excavatedsoil not gen zone cshell p1 outer_rad 0 0 p2 0 0 z1 p3 0 mouter_rad 0 size 4 _nz 54 & rat 1 1 1 dim inner_rad inner_rad inner_rad inner_rad group outer_concrete range group excavatedsoil not group inner_concrete not gen zone cshell p1 c1rad 0 0 p2 0 0 z1 p3 0 mc1rad 0 size 16 _nz 54 & rat 0.98 1 1 dim outer_rad outer_rad outer_rad outer_rad ; pause gen zone cshell p1 c0rad 0 0 p2 0 0 z1 p3 0 mc0rad 0 size 22 _nz 54 & rat 1.02 1 1 dim c1rad c1rad c1rad c1rad
gen zone cshell p1 out_b 0 0 p2 0 0 z1 p3 0 mout_b 0 size 26 _nz 54 & rat 1.05 1 1 dim c0rad c0rad c0rad c0rad ;gen zone reflect normal (1,0,0) origin (0,0,0) ;gen zone reflect normal (0,1,0) origin (0,0,0) ; attach face range cyl end1 0 0 -1 end2 0 0 2 rad c1radm not ; pause expgrid geomodel.Flac3D ; impgrid geomodel.Flac3D ; ; pause plot create view00 plot add surf plot set rot 90 0 0 plot show ; pause range name _out cyl end1 0 0 0 end2 0 0 z1p rad out_b_m not ; --- mechanical model --- model mo prop bulk c_k shear c_g coh 30e3 fric 20 ; --- thermal model --- model th_iso prop cond _cond spec_heat _cv thexp 0. ini density 2.e3 ini t _initemp ; ; 力学边界条件 fix z range z 0.99 1.01 fix z range z -0.01 0.01 fix x range x -0.01 0.01 fix y range y -0.01 0.01 fix x y range _out ; ; assign initial stress state ; set grav 0 0 -10 ; ini szz -6.e6 ;grad 0 0 10000 ; ini sxx -6.e6 ;grad 0 0 10000 ; ini syy -6.e6 ;grad 0 0 10000 ; range name _out2 cyl end1 0 0 -1 end2 0 0 2 rad tempbcm not fix t _initemp range _out2 ; plot create view01 plot add con temp outline on plot set mag 1.95 plot set rot 90 0 0 plot set center 6. -4.5 0.5 plot show ; --- settings --- hist id = 1 Monitor1t hist id = 2 Monitor2t hist id = 3 Monitor3t hist id = 4 Monitor4t hist id = 5 Monitor5t ; plot set background white set plot jpg set plot jpg size (1024,768) quality 2 Positive_Freezing ; ; save jijidongjie.sav ; plot show 'Base' plot set back white plot his 1 skip 4 2 skip 4 3 skip 4 4 skip 4 5 skip 4 & xlabel '天数/d' ylabel '温度/℃' alias '温度-天数曲线' plot hard 'Base' file jijidongjie.jpg his write 1 2 3 4 5 skip 4 file jijidongjie.dat ;------------------------------ ; Excavation ;------------------------------ model th_null range group excavatedsoil model th_null range group inner_concrete model th_null range group outer_concrete ; ; ini temp 18. range group model_board ; cast outer concrete model hydrat th_hyd_concrete1 range group outer_concrete prop dens 2000 bulk 1e3 shear .7e3 range group outer_concrete prop cond 2. thex 1e-5 spec_heat 0.2 range group outer_concrete prop qmax 1e5 cement 330 b_const -1.114 t1_const 7.2e4 range group outer_concrete prop gas_const 8.314 e_activate 33.5e3 cte_alpha 0.20 cte_tension 2.0 range group outer_concrete prop cte_bulk 0.98e3 cte_shear 0.50e3 c_const 0.4 a_const 0.6 range group outer_concrete prop dalpha_min 1e-4 cte_young 1e3 range group outer_concrete ; ; ini temp 18. range group outer_concrete ; 外壁的内侧 range name innerofoutercrete cyl end1 0 0 0 end2 0 0 1. rad inner_radp ; fix t 10. range innerofoutercrete ; apply flux=50. range innerofoutercrete ; model hydrat th_hyd_concrete1 range group inner_concrete prop dens 2000 bulk 1e3 shear .7e3 range group inner_concrete prop cond 2. thex 1e-5 spec_heat 0.2 range group inner_concrete prop qmax 1e5 cement 330 b_const -1.114 t1_const 7.2e4 range group inner_concrete prop gas_const 8.314 e_activate 33.5e3 cte_alpha 0.20 cte_tension 2.0 range group inner_concrete prop cte_bulk 0.98e3 cte_shear 0.50e3 c_const 0.4 a_const 0.6 range group inner_concrete prop dalpha_min 1e-4 cte_young 1e3 range group inner_concrete ; ini temp 18. range group inner_concrete ; range name inner_of_model_board cyl end1 0 0 -1. end2 0 0 2. rad 3.69 ; fix t _initemp range inner_of_model_board ; ; fix t _initemp range group model_board ; fix x y range inner_of_model_board ; ; ; model th_iso range group poam_board ; prop cond 0.1 thexp 1.e-8 spec_heat 1700 range group poam_board ; ; ini temp _initemp range group poam_board ; set geom_rep 600 ; def monitor_inner _ptgp01 = gp_near(4.175, 0., 0.) monitor_inner = gp_temp(_ptgp01) end def monitor_outer _ptgp02 = gp_near(4.8, 0., 0.) monitor_outer = gp_temp(_ptgp02) end ; hist id = 6 monitor_inner hist id = 7 monitor_outer ; plot show 'view01' plot set mag 4.0 ; ; set thermal implicit off Maintain_Freezing ; plot show 'Base' plot set background white plot his 1 skip 4 2 skip 4 3 skip 4 4 skip 4 5 skip 4 6 skip 4 7 skip 4 & xlabel '天数/d' ylabel '温度/℃' alias '温度-天数曲线' plot hard 'Base' file final.jpg
his write 1 2 3 4 5 6 7 skip 4 file final.dat ; save final.sav ; return ; Shut down computer by the command below ; system ("shutdown -s -t 15")
; 计算代码结束
图2 冻结温度场分布云图
代码在计算过程中要读入的数据文件共有3个，
文件1为：p_temp.dat
以下为该文件的内容：
-16. -16.5 -17. -17.5 -18. -18.5 -19. -19. -20. -20.5 -21. -21.5 -22. -22.5 -23. -23.5 -24. -24.5 -25. -25.5 -26. -26.5 -27. -27. -27.5 -27.5 -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -28. -30. -30. -30. -30. -30. -30. -30.
另外文件2是：m_temp.dat
内容如下：
-28. -28. -27. -27. -26. -26. -25. -25. -24. -24. -23. -23. -22. -22. -21. -21. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20. -20.
文件3读入冻结管位置坐标文件：c1.dat
内容如下：
7.6500e+000 7.5338e+000 7.1886e+000 6.6251e+000 5.8602e+000 4.9173e+000 3.8250e+000 2.6165e+000 1.3284e+000 4.6841e-016 ; 0.0000e+000 -1.3284e+000 -2.6165e+000 -3.8250e+000 -4.9173e+000 -5.8602e+000 -6.6251e+000 -7.1886e+000 -7.5338e+000 -7.6500e+000
免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删
155-2731-8020
