about to introduce phase change

demo.m

 close all; clear all; clc;
 %% PAREMETERS
 x_left = 0.0;
-x_right = 1.0;
-rho = 1.0;
-c_p = 1.0;
-k   = 1.0;
-T_0 = 0.0;
-
-dt = 0.01;
-num_els = 100;
-num_nodes = num_els + 1;
-num_time_steps = 5;
+x_right = 0.1;
+%A = 0.01^2;
+global rho c_p k diffusivity lamma L T_m T_l T_s T_0;
+rho = 4510;
+c_p = 520;
+k   = 16;
+diffusivity = k / rho / c_p;
+lamma = 0.388150542167233;
+L = 325;
+T_m = 1670;
+T_l = 2000;
+T_s = 1500;
+T_0 = T_s;
+time = 0.0;
+S   = 4;
 
-wait_time_plots = 1;%seconds
+dt = 1;
+num_els = 1000;
+num_nodes = num_els + 1;
+num_time_steps = 2;
+wait_time_plots = 0.5;%seconds
 
 %% PAREMETERS
 Xs = linspace(x_left,x_right,num_nodes);
+node2dof = [-1, 1:(num_nodes-2), -1];
+% node2dof = 1:(num_nodes);
+mask = find(node2dof>0);
+forced_nodes = find(node2dof<0);
+num_dofs = sum(node2dof>0);
 
 %% INITIAL CONDITION
-f0  = @(x) sin(16*x).^2 + (0.5 +0.5*cos(9*x));
+f0  = @(x) T_0*ones(size(Xs))';
+fl   = @(tem) 1/2*(tanh(S*2/(T_l - T_s)*(tem - (T_s + T_l)/2.0))+1.0);
+flp  = @(tem) S/(T_l - T_s)*(1 - (tanh(S*2/(T_l - T_s)*(tem - (T_s + T_l)/2.0))).^2);
+flpp = @(tem) -4*S^2/(T_l - T_s)^2*(tanh(S*2/(T_l - T_s)*(tem - (T_s + T_l)/2.0)).*sech(S*2/(T_l - T_s)*(tem - (T_s + T_l)/2.0)).^2);
 u_np1 = f0(Xs');
 
-ylims = [min(u_np1), max(u_np1)];
+ylims = [T_s-100, T_l+100];
 
 fig = figure;
 plot(Xs, u_np1);
@@ -32,27 +49,67 @@ hold on
 %% TIME LOOP
 for time_iter=1:num_time_steps
   u_n = u_np1;
+  time = time + dt;
+  u_exact = exact_sol(Xs,time);
+  dirichlet_vals = nan([num_nodes,1]);
+  dirichlet_vals(1) = T_l;
+  dirichlet_vals(end) = T_s;
   % ASSEMBLY
-  A = sparse(num_nodes,num_nodes);
-  M = sparse(num_nodes,num_nodes);
-  b = zeros([num_nodes,1]);
+  A = sparse(num_dofs,num_dofs);
+  M = sparse(num_dofs,num_dofs);
+  b = zeros([num_dofs,1]);
   for ielem=1:num_els
     inodes_global = [ielem, ielem+1];
     xr = Xs(ielem+1);
     xl = Xs(ielem);
     h  = xr - xl;
     M_loc = h / 6.0 * [2.0, 1.0; 1.0, 2.0];
-    M_loc = M_loc * rho * c_p;
+    M_loc = M_loc * rho * c_p / dt;
     A_loc = 1 / h * [1.0, -1.0; -1.0, 1.0];
     A_loc = A_loc * k;
-    M(inodes_global,inodes_global) = M(inodes_global,inodes_global) + M_loc;
-    A(inodes_global,inodes_global) = A(inodes_global,inodes_global) + A_loc;
+    for inode=1:2
+      inodeg = inodes_global(inode);
+      idof = node2dof(inodeg);
+      if (idof<0)
+        continue;
+      end
+      for jnode=1:2
+        jnodeg = inodes_global(jnode);
+        jdof = node2dof(jnodeg);
+        if (jdof>0)
+          M(idof,jdof) = M(idof,jdof) + M_loc(inode,jnode);
+          A(idof,jdof) = A(idof,jdof) + A_loc(inode,jnode);
+        else
+          dval = dirichlet_vals(jnodeg);
+          b(idof) = b(idof) - A_loc(inode,jnode)*dval;
+          b(idof) = b(idof) - M_loc(inode,jnode)*u_n(jnodeg);
+        end
+      end
+    end
   end
   % SOLVE
-  b = b + (1 / dt) * M * u_n;
-  u_np1 = mldivide((M / dt + A), b);
+  b = b + M * u_n(mask);
+  stiffness_matrix = M + A;
+  u_np1(mask) = mldivide(stiffness_matrix, b);
+  u_np1(forced_nodes) = dirichlet_vals(forced_nodes);
   clf(fig);
-  plot(Xs, u_np1);
+  plot(Xs, u_np1,'DisplayName','uh');
+  hold on
+  plot(Xs, u_exact,'DisplayName','exact');
   ylim(ylims);
+  legend;
   pause(wait_time_plots)
 end
+
+function [sol] = exact_sol(x,t)
+  global rho c_p k diffusivity lamma L T_m T_l T_s T_0;
+  interface_pos = locate_interface(t);
+  sol = zeros(flip(size(x)));
+  sol(x<=interface_pos) = T_l - (T_l - T_m)*erf(x(x<=interface_pos)/(2*sqrt(diffusivity*t)))/erf(lamma);
+  sol(x>interface_pos) = T_s + (T_m - T_s)*erfc(x(x>interface_pos)/(2*sqrt(diffusivity*t)))/erfc(lamma);
+end
+
+function [interface_pos] = locate_interface(t)
+  global rho c_p k diffusivity lamma L T_m T_l T_s T_0;
+  interface_pos = 2*lamma*sqrt(diffusivity*t);
+end