基于FLAC3D軟件的人工凍結溫度場分析
人工地下凍結(AGF, Artificial Ground Freezing)最早出現于俄國，用于金礦開采，后由德國工程師用于煤礦礦井建設并獲得專利，技術漸趨成熟，現在已廣泛應用于地鐵、深基坑、礦井建設等工程中。人工地下凍結法的原理是將人工凍結管埋置于土體內，利用凍結管內循環的冷媒劑將土體中的熱量逐漸帶走從而形成強度、高密封性的凍結土壁，形成的凍結土壁通常構成封閉的環形，能夠起到承受荷載和密封防水的作用。AGF法的適應性強、安全 可靠、無污染，但和其他方法相比，一般造價較高。
AGF法進行土體凍結，可以適用于存在潛在地下水的地層，為保證地下開挖施工的安全，先開鑿凍結孔，下放凍結管并將凍結管與負責冷量循環的凍結站相連，待凍結站工作后每根凍結管周圍會逐漸形成凍結圓柱，相鄰凍結圓柱半徑逐漸增大到相互連接后就形成了閉合的凍結墻(frozen wall)，凍結墻外側的潛在地下水被封閉，同時凍結土體的強度比原土高很多，都有利于保證后續凍結墻內向下開挖土體施工的安全。
FLAC3D軟件作為有限差分軟件，除廣泛應用于巖土體的力學問題分析外，還可以用于凍結溫度場的分析。本文后續演示了采用人工地下凍結法進行煤礦井筒開挖的分析代碼，是我在碩士階段編寫的，雖然代碼還有不少值得完善的地方。但里面的包含了諸如坐標數據的讀入、自動截圖和數據文件的自動取名保存等等。 希望對大家學習FLAC3D有幫助。
凍結方式為豎向凍結，冷媒是鹽水。
圖1 人工凍結法凍結壁形成示意圖
代碼：
; ----------------------------------------
; 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
工程師必備
- 項目客服
- 培訓客服
- 平臺客服
TOP
