_node=(0:(N-1))*delta_x;
x_s(2:N)=x_node(2:N)-delta_x/2;
x_n(1:(N-1))=x_node(1:(N-1))+delta_x/2;
d_lens(1:N)=0;
%存储初始变量
T1(1:N)=2;P1(1:N)=0;
nn=1;tao=0;lens_num=1;Nnn=1;new_lens=0;
Ttao(Nnn)=tao;X_NODE(:,Nnn)=x_node';T(:,Nnn)=T1';P(:,Nnn)=P1';
tic
n_stop=24*1*3600/3;
taoc=0;
while nn<n_stop
nn=nn+1
tao=tao+delta_tao;
P0=P1;T0=T1;
%定义新的 T1,P1,也就是 n+1 时刻边界条件
T_warm=2;T_cold=-4;P_warm=0;
%计算一个时间步,因为用全隐式,要迭代两次,计算完后再考虑透镜体的形成
for sb=1:1
%定义各结点的导热及导湿系数,并计算控制容积界面上的数值;给出每个控制容积的等效热容
clear lamtaii;clear kii;clear Cvii;
lamtaii=lamta(P1,T1);kii=k_water(P1,T1);
ii=2:(N-1);
lamta_s(ii)=(delta_x+d_lens(ii-1))./(delta_x/2./lamtaii(ii)+delta_x/2./lamtaii(ii-1)+d_lens(ii-1)/lamta_i);
lamta_n(ii)=(delta_x+d_lens(ii))./(delta_x/2./lamtaii(ii)+delta_x/2./lamtaii(ii+1)+d_lens(ii)/lamta_i);
k_s(ii)=2./(1./kii(ii)+1./kii(ii-1));
k_n(ii)=2./(1./kii(ii)+1./kii(ii+1));
Cvii(ii)=(Cv(P1(ii),T1(ii))*delta_x+C_i*rou_i*d_lens(ii))./(delta_x+d_lens(ii));
clear D;
%给出系数矩阵,用 5 列表示,第 3 列表示对角元素
D=zeros(2*N1+N2+1,5);
%node 1,未定义的自动为 0
D(1,3)=1;d(1)=T_warm;
D(2,3)=1;d(2)=P_warm;
%node2~N1
ii=2:N1;
a11=rou_i*delta_v*T_a*F_diff(P1(ii),T1(ii));
a12=Cvii(ii)+L_f*rou_i*F_diff(P1(ii),T1(ii));
a21=(rou_w-rou_i)*delta_v*T_a*F_diff(P1(ii),T1(ii))/L_f/rou_w;
a22=(rou_w-rou_i)*F_diff(P1(ii),T1(ii))/rou_w;
D(2*ii-1,1)=-lamta_s(ii)'*delta_tao/delta_x^2;
D(2*ii-1,3)=a12'+(lamta_s(ii)'+lamta_n(ii)')*delta_tao/delta_x^2;
D(2*ii-1,4)=a11';
109D(2*ii-1,5)=-lamta_n(ii)'*delta_tao/delta_x^2;
d(2*ii-1)=a11.*P0(ii)+a12.*T0(ii);
D(2*ii,1)=-k_s(ii)'*delta_tao/delta_x^2;
D(2*ii,2)=a22';
D(2*ii,3)=a21'+(k_n(ii)'+k_s(ii)')*delta_tao/delta_x^2;
D(2*ii,5)=-k_n(ii)'*delta_tao/delta_x^2;
d(2*ii)=a21.*P0(ii)+a22.*T0(ii);
if N2==1
D(2*N1+1,3)=1;d(2*N1+1)=T_cold;
D(2*N1+2,2)=L_f/v_w/T_a;D(2*N1+2,3)=-1;d(2*N1+2)=0;%若 N2=1,则结束
else
% N1+1 node
D(2*N1+1,1)=-lamta_s(N1+1)/delta_x;
D(2*N1+1,2)=-L_f*rou_w*k_s(N1+1)/delta_x;
D(2*N1+1,3)=lamta_s(N1+1)/delta_x+lamta_n(N1+1)/(x_node(N1+2)-x_node(N1+1))+Cvii(N1+1)/delt
a_tao*(x_n(N1+1)-x_s(N1+1));
D(2*N1+1,4)=L_f*rou_w*k_s(N1+1)/delta_x;
D(2*N1+1,5)=-lamta_n(N1+1)/(x_node(N1+2)-x_node(N1+1));
d(2*N1+1)=Cvii(N1+1)*(x_n(N1+1)-x_s(N1+1))*T0(N1+1)/delta_tao;
D(2*N1+2,2)=L_f/v_w/T_a;D(2*N1+2,3)=-1;d(2*N1+2)=0;
if N2==2
D(2*N1+N2+1,3)=1;d(2*N1+N2+1)=T_cold;
elseif N2==3
D(2*N1+3,1)=-lamta_s(N1+2)/(x_node(N1+2)-x_node(N1+1));
D(2*N1+3,3)=(x_n(N1+2)-x_s(N1+2))*Cvii(N1+2)/delta_tao+lamta_s(N1+2)/(x_node(N1+2)-x_node(
N1+1))+lamta_n(N1+2)/(x_node(N1+3)-x_node(N1+2));
D(2*N1+3,4)=-lamta_n(N1+2)/(x_node(N1+3)-x_node(N1+2));
d(2*N1+3)=(x_n(N1+2)-x_s(N1+2))*Cvii(N1+2)*T0(N1+2)/delta_tao;
D(2*N1+N2+1,3)=1;d(2*N1+N2+1)=T_cold;
else
D(2*N1+3,1)=-lamta_s(N1+2)/(x_node(N1+2)-x_node(N1+1));
D(2*N1+3,3)=(x_n(N1+2)-x_s(N1+2))*Cvii(N1+2)/delta_tao+lamta_s(N1+2)/(x_node(N1+2)-x_node(
N1+1))+lamta_n(N1+2)/(x_node(N1+3)-x_node(N1+2));
D(2*N1+3,4)=-lamta_n(N1+2)/(x_node(N1+3)-x_node(N1+2));
d(2*N1+3)=(x_n(N1+2)-x_s(N1+2))*Cvii(N1+2)*T0(N1+2)/delta_tao;
ii=(N1+3):(N1+N2-1);
D(N1+ii+1,2)=-lamta_s(ii)'./(x_node(ii)'-x_node(ii-1)');
D(N1+ii+1,3)=(x_n(ii)'-x_s(ii)').*Cvii(ii)'/delta_tao+lamta_s(ii)'./(x_node(ii)'-x_node(ii-1)')+lamta_n(ii)'
110/(x_node(ii+1)'-x_node(ii)');
D(N1+ii+1,4)=-lamta_n(ii)'./(x_node(ii+1)'-x_node(ii)');
d(N1+ii+1)=(x_n(ii)'-x_s(ii)').*Cvii(ii)'.*T0(ii)'/delta_tao;
%N1+N2 node
D(2*N1+N2+1,3)=1;d(2*N1+N2+1)=T_cold;
end
end
clear AA;clear BB;clear CC;
%五对角阵的求解算法
AA(1)=-D(1,5)/D(1,3);BB(1)=-D(1,4)/D(1,3);CC(1)=d(1)/D(1,3);
AA(2)=-D(2,5)/(D(2,3)+BB(1)*D(2,2));BB(2)=(-D(2,4)-AA(1)*D(2,2))/(D(2,3)+BB(1)*D(2,2));CC(2)=
(d(2)-CC(1)*D(2,2))/(D(2,3)+BB(1)*D(2,2));
for ii=3:(2*N1+N2+1)
SB=-D(ii,2)-BB(ii-2)*D(ii,1);
AA(ii)=-D(ii,5)/(D(ii,3)+AA(ii-2)*D(ii,1)-SB*BB(ii-1));
BB(ii)=(-D(ii,4)+SB*AA(ii-1))/(D(ii,3)+AA(ii-2)*D(ii,1)-SB*BB(ii-1));
CC(ii)=(d(ii)-CC(ii-2)*D(ii,1)+SB*CC(ii-1))/(D(ii,3)+AA(ii-2)*D(ii,1)-SB*BB(ii-1));
end
fai(2*N1+N2+1)=CC(2*N1+N2+1);
fai(2*N1+N2)=BB(2*N1+N2)*fai(2*N1+N2+1)+CC(2*N1+N2);
for ii=1:(2*N1+N2-1)
jj=2*N1+N2-ii;
fai(jj)=AA(jj)*fai(jj+2)+BB(jj)*fai(jj+1)+CC(jj);
end
%赋值
ii=1:(N1+1);
T1(ii)=fai(2*ii-1);
P1(ii)=fai(2*ii);
if N2>=2
ii=(N1+2):(N1+N2);
T1(ii)=fai(ii+N1+1);
end
end
%计算活动透镜体生长速度,即冻胀速度,改变节点及界面坐标
V_H(nn-1)=rou_w/rou_i*k_n(N1)*(P1(N1)-P1(N1+1))/delta_x;
delta_s=V_H(nn-1)*delta_tao;
if N2>=2
x_node((N1+2):N)=x_node((N1+2):N)+delta_s;
x_s((N1+2):N)=x_s((N1+2):N)+delta_s;
x_n((N1+1):(N-1))=x_n((N1+1):(N-1))+delta_s;
d_lens(N1+1)=d_lens(N1+1)+delta_s;
end
%判断有无新透镜体形成,若有则确定其位置,改变活动透镜体底端压力
111分离冰冻胀模型计算 Matlab 程序(以 Konrad 恒温冻结试验为例)
主程序
format long;
%基本固定常物性参数
global L_f;L_f=333.6*10^3;global T_a;T_a=273.15;global g;g=9.8;%相变潜热%热力学温度基准值%
重力加速度
global rou_w;rou_w=1000;global rou_i;rou_i=900;global rou_d;rou_d=1670;%液态水密度%固态冰密
度%土体干密度
global v_w;v_w=1/rou_w;global v_i;v_i=1/rou_i;global delta_v;delta_v=v_i-v_w;%液态水比容%固态
冰比容%比容差
global lamta_i;lamta_i=2.32;global lamta_w;lamta_w=0.58;global lamta_soil;lamta_soil=1.95;%冰的导
热系数%水的导热系数%土骨架的导热系数
global C_w;C_w=4.18*10^3;global C_i;C_i=2.09*10^3;global C_s;C_s=2.2*10^3;%液态水的比热%
固态冰的比热%土骨架的比热
global sigma_iw;sigma_iw=7.3*10^(-2)/2.2; %冰水界面的张
力
%与未冻水膜及土颗粒有关参数
global R_cl;R_cl=2.5*10^(-6);global R_si;R_si=40*10^(-6);global R_sa;R_sa=100*10^(-6);
global m_cl;m_cl=0.29;global m_si;m_si=0.61;global m_sa;m_sa=0.1;
global K_av;K_av=-(m_cl/R_cl+m_si/R_si+m_sa/R_sa);
%饱和质量含水量和饱和体积含水量
global theta_sat;theta_sat=0.3698;global w_tot;w_tot=rou_w*theta_sat/rou_d;
%黏土、粉土、沙土的颗粒尺寸界限常数(mm)
d_cl=0.002;d_si=0.026;d_sa=0.12;
%平均颗粒尺寸及几何方差偏移
d_g=exp(m_cl*log(d_cl)+m_si*log(d_si)+m_sa*log(d_sa));
sigma_g=exp((m_cl*(log(d_cl))^2+m_si*(log(d_si))^2+m_sa*(log(d_sa))^2-d_g^2)^0.5);
%空气进入势及指数参数 b
b=d_g^(-0.5)+0.2*sigma_g;
global k_sat;k_sat=4*10^(-5)*(0.5/(1-theta_sat))^(1.3*b)*exp(-6.88*m_cl-3.63*m_si-0.025)/rou_w/g;
%饱水导湿系数!!!!!!!!!!!!!!!!!
%未冻水含量曲线
global A;A=0.07;
global B;B=0.33;
global T_f;T_f=-(w_tot/A)^(1/-B);
%临界分离压力
global P_c;P_c=0.025*10^6;
%网格参数
long=0.1;N=1001;delta_x=long/(N-1);delta_tao=3;
N1=N-1;N2=1;
108lear P_sep;
frozen_node=find(theta_u(P1,T1)<theta_sat);
front_node=min(frozen_node);
f_node(nn)=front_node;
if front_node<(N1+1)
P_sep=v_w/v_i*P1(front_node:N1)-sigma_iw*K_av-L_f*T1(front_node:N1)/v_i/T_a;
mm=max(P_sep);
if mm>P_c
lens_posi_pre=find(P_sep>P_c);
N1=min(lens_posi_pre)+front_node-2;N2=N-N1;
lens_index1(lens_num,:)=[nn,front_node,N1+1]
lens_index2(lens_num,:)=[P1(N1+1),T1(N1+1),P_sep(lens_posi_pre(1))]
lens_num=lens_num+1;
lens_time(N1+1)=tao;P1(N1+1)=L_f/v_w/T_a*T1(N1);
end
end
lin_depth(nn)=x_node(N)-x_node(min(find(T1<=0)));
%选择存储相关变量
if new_lens | mod(nn,10)==1
Nnn=Nnn+1;
Ttao(Nnn)=tao;X_NODE(:,Nnn)=x_node';
T(:,Nnn)=T1';P(1:(N1+1),Nnn)=P1(1:(N1+1))';
new_lens=0;
end
end
toc
辅助程序(work 目录下被调用函数)
function zanswer=Cv(z_p,z_t)
%formula from dongtu wulixue
global w_tot;global rou_d;global rou_w
global C_s;global C_w;global C_i;
zanswer=(C_s+(w_tot-theta_u(z_p,z_t)*rou_w/rou_d)*C_i+theta_u(z_p,z_t)*rou_w/rou_d*C_w)*rou_d
;
function zanswer=F_diff(z_p,z_t)
global rou_d;global rou_w;
global A;global B;
global delta_v;global T_a;global L_f;
global theta_sat;
z_pt=z_t+delta_v*T_a*z_p/L_f;
zanswer_pre=rou_d/rou_w*A*B*(-z_pt).^(-B-1);
kk1=find(theta_u(z_p,z_t)<theta_sat); zanswer(kk1)=zanswer_pre(kk1);
kk2=find(theta_u(z_p,z_t)==theta_sat);zanswer(kk2)=0;
function zanswer=k_water(z_p,z_t)
global k_sat;global theta_sat;
112gama=9; %指数参数,决定导湿系数衰减快慢的!!!!!!!!!!!!
zanswer=k_sat*(theta_u(z_p,z_t)/theta_sat).^gama;
