% 参数设置
L = 100;% 土柱长度，单位：cm
W = 50;% 土柱宽度，单位：cm
t_end = 360024365;% 模拟时间，单位：s
dt = 36002430;% 时间步长，单位：s
dz = 1;% 空间步长，单位：cm
z = 0:dz:L;% 空间网格
m = length(z);% 空间网格数
theta_f = zeros(m,1);% 孔隙中水含量初始值
theta_m = 0.1*ones(m,1);%土壤基质中水含量初始值
T = 5*ones(m,1);% 土柱温度初始值
S_f = zeros(m,1);% 孔隙中水源项初始值
S_m = zeros(m,1);% 土壤基质中水源项初始值
K_f = 1e-6;% 孔隙中渗透系数
K_m = 1e-6;% 土壤基质中水渗透系数
w = 0.3;% 孔隙率
rho_c = 2352;% 土壤比热容
alpha = 1.28;% 土壤热扩散系数
H = 10;% 上边界水头，单位：cm
T_H = 10;% 上边界温度，单位：°C
B = 0;% 左边界水头，单位：cm
T_B = 1;% 左边界温度，单位：°C
D = 0;% 右边界水头，单位：cm
T_D = 1;% 右边界温度，单位：°C
theta_0 = 0.1;% 初始含水率
% 计算边界温度
T_L = zeros(m,1); T_L(1) = T_B;% 左边界温度
T_L(m) = T_D;% 右边界温度
for i = 2:m-1
T_L(i) = (2*alpha/dz^2/rho_c*T(i-1) + 2*alpha/dz^2/rho_c*T(i+1) + (1-w)*rho_c*H/alpha) / (2*alpha/dz^2*rho_c+ (1-w)*rho_c/alpha); end
% 创建矩阵A和B
A = zeros(3*m,3*m); A(1,1) = 1;% 上边界
A(m,m) = 1;% 下边界
B = zeros(3*m,1);
for i = 2:m-1 % 孔隙中水方程系数矩阵
A(i,i-1) = -K_f/dz^2;
A(i,i) = 2*K_f/dz^2 + 1/dt;
A(i,i+1) = -K_f/dz^2; B(i) = S_f(i);
% 土壤基质水方程系数矩阵
A(i+m,i-1+m) = -K_m/dz^2;
A(i+m,i+m) = 2*K_m/dz^2 + 1/dt;
A(i+m,i+1+m) = -K_m/dz^2;
B(i+m) = S_m(i+m);
% 孔隙中水源项和基质水源项
S_f(i) = 1/dt - theta_f(i)/dt*(theta_f(i)-theta_0)/dz - (1-theta_f(i))/dt*(theta_f(i)-theta_m(i))/dz ...
+ K_f*(H-theta_f(i)/w)/dz;
S_m(i+m) = 1/dt - theta_m(i+m)/dt*(theta_m(i+m)-theta_f(i))/((1-w)*dz) ...
- (1-theta_m(i+m))/dt*(theta_m(i+m)-theta_0)/dz ...
+ K_m*(theta_f(i+m)/w-theta_m(i+m))/(1-w)/dz;
% 温度方程系数矩阵
A(i+2*m,i-1+2*m) = -alpha/dz^2/rho_c;
A(i+2*m,i+2*m) = 2*alpha/dz^2/rho_c + 1/dt;
A(i+2*m,i+1+2*m) = -alpha/dz^2/rho_c;
% 温度源项
B(i+2*m) = T(i)/dt*(1/dz^2/rho_c*(theta_f(i)-theta_m(i)) ...
- 2/dz^2*(alpha*(T(i+1)-T(i))-alpha*(T(i)-T(i-1))));
end
% 边界条件
A(1,1) = 1;
A(m,m) = 1;
B(1) = T_H;
B(m) = T_D;
% 求解方程组
sol = A\B;
theta_f = sol(1:m);
theta_m = sol(m+1:2*m);
T = sol(2*m+1:3*m);
% 绘图
plot(z,theta_f,'r',z,theta_m,'b',z,T,'k')
xlabel('Depth (cm)')
ylabel('Value')
title('Results')
错误信息
索引超出数组元素的数目(101)。
出错 test22 (line 46)
B(i+m) = S_m(i+m);
- Badges
- Science topic
- Similar topics
- Condensed Matter Physics
- Zero