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

2 次查看（过去 30 天）
Terry Keeley
Terry Keeley 2022-4-5
回答： Shivam 2023-10-20
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

请先登录，再回答此问题。

回答（1 个）

Shivam
Shivam 2023-10-20
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

类别

在 Help Center 和 File Exchange 中查找有关 Gas Dynamics 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by