% 定义常数和参数 dt = 0.1;% 时间步长 dx = 0.1;% 空间步长 L = 1;% 空间长度 最大温度 = 100;% 最大模拟时间 Nt = 最大/分;% 时间步数 Nx = L/dx;% 空间步数 RHO = 1;% 密度 C = 1;% 热容 λ = 1;% 热导率 L = 1;% 潜热 rho_l = 1;% 液体密度 rho_w = 1;% 水密度 D = 1;% 扩散系数 k = 1;% 热对流系数 % 初始化温度和液相温度 T = 零（Nx+1， Nt+1）;T（：，1） = 0;% 初始温度为0 theta_l = 零（Nx+1， Nt+1）;theta_l（：，1） = 0;% 初始液相温度为0 % 迭代求解 对于 n = 1：Nt % 求解温度方程 对于 i = 2：Nx T（i，n+1） = T（i，n） + dt/rho/C/dx^2 * lambda * （T（i+1，n） - 2 T（i，n） + T（i-1，n）） ... + dt L rho_l/rho/C * （theta_l（i，n+1） - theta_l（i，n））; 结束 % 求解液相温度方程 对于 i = 2：Nx theta_u = T（i，n）;% 上层温度即为该位置温度 theta_z = T（i，n） - theta_l（i，n）;% 上下层温度差 theta_l（i，n+1） = theta_l（i，n） + dt/rho_w/rho_l/dx^2 * D * （theta_l（i+1，n） - 2theta_l（i，n） + theta_l（i-1，n）） ... + 分rho_w * k * theta_z;结束 结束 % 绘制温度随时间和位置的变化 [x， t] = meshgrid（0：dx：L， 0：dt：Tmax）;数字;冲浪（x， t， t'）;xlabel（'位置'）;ylabel（'时间'）;zlabel（'温度'）;title（'温度随时间和位置的变化'）;% 绘制液相温度随时间和位置的变化 数字;冲浪（x， t， theta_l'）;xlabel（'位置'）;ylabel（'时间'）;zlabel（'液相温度'）;title（'液相温度随时间和位置的变化'）;为以上代码添加并应用边界条件的代码

边界条件可以根据具体问题来选择，以下是两种常见的边界条件：

1. Dirichlet 边界条件：在边界处给定温度值

对于左右边界，可以将其温度固定为0，即：

T(1,n) = 0;
T(Nx+1,n) = 0;

对于上下边界，可以将其温度固定为一个常数，比如最大温度，即：

T(:,1) = 最大温度;
T(:,Nt+1) = 最大温度;

对于液相温度theta_l，可以将其边界条件设置为与温度T的边界条件相同，即：

theta_l(1,n) = 0;
theta_l(Nx+1,n) = 0;
theta_l(:,1) = 最大温度;
theta_l(:,Nt+1) = 最大温度;

2. Neumann 边界条件：在边界处给定热流量（温度梯度）

对于左右边界，可以将其热流量固定为0，即：

T(1,n+1) = T(2,n+1);
T(Nx+1,n+1) = T(Nx,n+1);

对于上下边界，可以将其热流量固定为0，即：

T(:,1) = T(:,2);
T(:,Nt+1) = T(:,Nt);

对于液相温度theta_l，同样可以将其边界条件设置为与温度T的边界条件相同，即：

theta_l(1,n+1) = theta_l(2,n+1);
theta_l(Nx+1,n+1) = theta_l(Nx,n+1);
theta_l(:,1) = 最大温度;
theta_l(:,Nt+1) = 最大温度;

需要注意的是，在迭代求解过程中，对于液相温度方程，由于theta_l的边界条件不一定与T的边界条件相同，因此需要根据具体问题来确定液相温度的边界条件。