How to use an interpolated value as index in a for loop? Stefan problem

Hi,

I am attempting to solve a 1D stefan problem using the improved explicit enthalpy method algoritm from Voller and Cross 1980, which can be found: https://doi.org/10.1016/0017-9310(81)90062-4.

The aim is to be able to track the moving phase change boundary between solid and liquid. As the boundary doesn't necessarily move one spatial step with each time step, they propose using linear interpolation to find the time the boundary is on a node.

The time that the boundary is on the node is ti = (n+X)*dt, where n is the current time step, X is the linearly interpolated value (X<0) and dt is the time step.

However the next step in the method is calculating the temperatures at the nodes when at time ti.

The temperature at time ti for node i will be the melting temperature of the material and they provide an equation for calculating the temperatures at the neighbouring nodes at time ti, which I have attached. Note: in the eqn j is the time step (n in my example) and k is the spatial step (i in my example).

But I can't find a way to assign those values, as the temperatures are occuring at a time between two indices, n and n+1, i.e. n*dt < ti < (n+1)*dt.

When trying to set the value of temperature at that time using indexing, matlab throws an error saying 'Indicies must be positive integers'. As for example, the first time the boundary is at a node is t= 15.09 hours and matlab doesn't allow you to input u(15.09, i)=T_m.

The rest of the code runs fine but my problem is with the last if statement, I am unsure how I first can assign the calulated temperatures to the temperature vector and then use those values when calculating the enthalpy in the next time step.

% Script to solve the simple Stefan problem using an explicit enthalpy
% method as shown by Voller and Cross - 1980
clear

% set grid sizes, and create matrices etc
xmax = 10;
dx=0.125; % 12.5 cm
x = (0:dx:xmax);
nx=length(x);

tmax = (86400*50);
dt = 3600; % 1 hour = 3600seconds
t=(0:dt:tmax);
nt = length(t);
u = zeros(nt, nx);
H = zeros(nt, nx);

% Load properties of block 
k=2; % Thermal conductivity W/m degC
cp=2.5e6; % Specific heat J/kg deg C
rho=1; % Density kg/m3
T_m = 0; % Melting Temp deg c
Lf = 100e6; % J/kg - Latent heat
pcm_alph = (k/(rho*cp));
p = ((pcm_alph * dt)/(dx)^2); % Fourier number to check stability of explicit method

% Initial conditions
T_0 = 2; % degC
u(1,:) = T_0; % Initial temperature across the block
if u(1,:) <= T_m
    H(1,:) = cp*(u(1,:)); 
elseif u(1,:) > T_m
    H(1,:) = ((cp*(u(1,:))) + Lf);
end

for n=1:nt-1
    for i=2:nx-1
        % LHS Boundary
        u(n+1,1)= -10; 
        if u(n+1,1) <= T_m
            H(n+1,1) = cp*(u(n+1,1));
        elseif u(n+1,1) > T_m
            H(n+1,1) = ((cp*(u(n+1,1))) + Lf);
        end
        % RHS Boundary 
        u(n+1,nx) = u(n+1,nx-1);
        if u(n+1,nx) <= T_m
            H(n+1,nx) = cp*(u(n+1,nx));
        elseif u(n+1,nx) > T_m
            H(n+1,nx) = ((cp*(u(n+1,nx))) + Lf);
        end
        
        % Enthalpy for interior nodes
        H(n+1,i) = H(n,i) + ((k/rho)*(dt/(dx^2)))*(u(n,i-1) - (2*u(n,i)) + u(n,i+1));
        
        if  H(n+1,i) < (cp*T_m)
            u(n+1,i) = (H(n+1,i)/cp);
        elseif (cp*T_m) <= H(n+1,i) && H(n+1,i) <= ((cp*T_m) + Lf)
            u(n+1,i) = T_m;
        elseif H(n+1,i)> ((cp*T_m) + Lf)
            u(n+1,i) = ((H(n+1,i)-Lf)/cp);
        end
        
        %%%%%%%%%%%%%%% Problem section %%%%%%%%%%%%%%%
        if H(n,i)> ((cp * T_m) + (Lf/2)) && H(n+1,i)< ((cp * T_m) + (Lf/2))
            loc(n)= i*dx;                                               % x location of node i, where interface falls
            X(n) = (((Lf/2) + (cp * T_m)) - H(n,i))/(H(n+1,i) - H(n,i)); % Linear interp to find time at node i
            inter_t(n) = (n+X(n))*dt;                          % Time when interface is at node i
            inter_hrs(n) = round((inter_t(n)/3600),2);          % Convert seconds to hours
   
            % X<1 and can't index to non integer
            u(n+X,i)= T_m;
            u(n+X,i-1) = (X*(u((n+1),(i-1)) - u(n,i-1)) + u(n,i-1));
            u(n+X,i+1) = (X*(u((n+1),(i+1)) - u(n,i+1)) + u(n,i+1)); 
        end
    end
end

I would be grateful for any suggestions.

Thanks

Answer 1

[Shivam]

Hi,

I understand you are getting an error while assigning the calculated temperature to the temperature vector and calculating enthalpy using this vector.

The problem you're encountering is due to the nature of MATLAB's matrix indexing. MATLAB only accepts integer indices, so you can't directly assign a value to a non-integer index like "n+X."

You can work around this by storing the interpolated temperatures in a separate array and then using these values to calculate the enthalpy in the next step.

I hope it helps.

Thanks,

Shivam