Why is my thermal lattice boltzmann model ignoring the source term? Attempting to model one-phase stefan problem (phase change) using an iterative enthalpy method.

Why is my LBM code ignoring the source term?

I am trying to solve the one-phase stefan problem of freezing using a D1Q2 LBM with an iterative enthalpy method. My code does not produce the expected results when compared to the analytical solution. It apparently solves the heat conduction equation only and ignores the liquid fraction (Fl) dependent source term in the collision step.

Any advice would be appreciated. The code is below.

%Thermal LBM D1Q2 lattice
%Liquid to solid one-phase phase change in 1D. Liquid initially at melting temperture. At t=0 liquid comes into contact
%with cold wall at T_c < T_m.
clear all;
%all variables in lattice units (LU)
dx=1; nx=98; Lx=nx*dx;
dt=1; nt=10^4; Lt=nt*dt;
x(1)=0; x(2:nx+1)=dx/2:dx:Lx-dx/2; x(nx+2)=Lx;
dt=1; nt=10^4; t=0:dt:nt;
alpha=0.5; %ice
L=1; Cps=1; St=1;
T_c=0; T_m=1;
tau_g=alpha+1/2;  %relaxation time
%enthalpy
Hs=Cps*T_m; Hl=Cps*T_m+L;
% Initial macroscopic temperature distribution and liquid fraction are given:
T=T_m*ones(nx+2,1); T(1)=T_c; %temperature
Fl=ones(nx+2,1); Fl(1)=0; %liquid fraction
H=Hl+zeros(nx+2,1); H(1)=Hs; %enthalpy
%initialize for tracking position of phase-change boundary
phi=zeros(1);
%analytical Neumann solution for parameter lambda
lambda=ste_nondim(St); %calls function. For these conditions lambda = 0.62
%Initialize distribution function
n_a=3; %for D2Q5 %n_a=9 for D2Q9
g=zeros(nx+2, n_a); %distribution function
%D1Q2 weights
w0=[0 1/2 1/2]; %weight factors for D2Q5 topology
for i=1:nx+2
    for a=1:n_a
        g(i,a)=w0(a)*T(i);
    end
end
g_eq=g;
%LBM loop
for steps=1:nt
    Fl_old=Fl;
    T_old=T;
    time=t(steps)
    %Collision
    epsilon=1e-8; error=1;
    while error>epsilon
        Fl_old_iter=Fl;
        T_old_iter=T;
        for i=1:nx+2
            if(Fl(i)<1) %only apply collision in parts of the domain containing the active phase
                for a=1:n_a
                    g_eq(i,a)=w0(a)*T(i);
                    g(i,a)=g(i,a)-(1/tau_g)*(g(i,a)-g_eq(i,a))-w0(a)*(1/St)*(Fl(i)-Fl_old(i));
                end
            else %bounceback at boundary
                boundary=i; %position of phase interface (where Fl<1)
                break
            end
        end
        %bounceback at phase interface
        g(boundary,2)=g(boundary,3);
        %boundary conditions
        g(1,2)=T_c-g(1,3); %constant temperature at wall
        g(nx+2,2)=g(nx+1,2); %adiabatic 
        g(nx+2,3)=g(nx+1,3); %adiabatic 
        %streaming g
        for i=nx+2:-1:2
            g(i,2)=g(i-1,2); %east to west
        end
        for i=2:nx+1
            g(i,3)=g(i+1,3); %west to east (but not i=1)
        end
        %obtain macroscopic variables
        %temperature
        for i=1:nx+2
            T_sum=0;
            for a=1:n_a
                T_sum=T_sum+g(i,a);
            end
            T(i)=T_sum;
        end
        %enthalpy
        for i=1:nx+2
            H(i)=Cps*T(i)+Fl(i)*L;
        end
        %liquid fraction
        for i=2:nx+2
            if(H(i)<Hs) %solid
                Fl(i)=0;
            elseif(H(i)>=Hs && H(i)<=Hs+L) %mushy
                Fl(i)=(H(i)-Hs)/(Hl-Hs);
            elseif(H(i)>Hs+L) %liquid
                Fl(i)=1;
            end
        end
        %check for convergence
        error_Fl=abs(max((Fl-Fl_old_iter))./Fl);
        error_T=abs(max((T-T_old_iter))./T);
        error=min([error_Fl, error_T]);
    end 
    %track freezing front position
    if any(Fl < 1) %if any has frozen yet
        position=max(find(Fl<=0.5));
        phi=[phi, position];
    end
end

Answer 1

[Abhishek Chakram]

Hi Gwendolyn Brouwer,

I see that you are trying to identify why your model is ignoring the source term. Here are few suggestions and potential corrections to the MATLAB code that you provided:

• The D1Q2 lattice model has only two velocities, so you should define only two weights. The weights should sum up to 1. For a D1Q2 lattice, the typical weights are “w0 = ½” and “w1 = ½”.
• Since you are using a D1Q2 lattice, you should initialize “g” with only two directions.
• The source term in the collision step seems to be subtracted from the distribution function. This might not be the correct approach if you are trying to model a source term as an addition to the equilibrium distribution. Ensure the sign and implementation of the source term are consistent with the physical model.
• The convergence check might exit the loop prematurely if “error_Fl” or “error_T” is zero. You should add a small number (like “eps”) to the denominator to prevent division by zero.
• The boundary conditions should be applied after the streaming step to ensure they are not overwritten.
• The way you track the freezing front position using “find” may not work correctly if the liquid fraction “Fl” does not have a value exactly equal to 0.5. Consider using a threshold value or interpolation to find the front position.

Here is an updated snippet of the code reflecting some of the above points:

% D1Q2 weights
w0 = [1/2 1/2]; % weight factors for D1Q2 topology

% Initialize distribution function
n_a = 2; % for D1Q2
g = zeros(nx+2, n_a); % distribution function

for i = 1:nx+2
    for a = 1:n_a
        g(i, a) = w0(a) * T(i);
    end
end

% rest of the code

% Collision
epsilon = 1e-8; error = 1;
while error > epsilon
    % rest of the collision code

    % Streaming g
    for i = nx+2:-1:2
        g(i, 1) = g(i-1, 1); % streaming for direction 1
    end
    for i = 1:nx+1
        g(i, 2) = g(i+1, 2); % streaming for direction 2
    end

    % Apply boundary conditions after streaming
    % boundary conditions code

    % rest of the LBM loop

    % Convergence check (add eps to prevent division by zero)
    error_Fl = abs(max((Fl - Fl_old_iter)) ./ (Fl + eps));
    error_T = abs(max((T - T_old_iter)) ./ (T + eps));
    error = max([error_Fl, error_T]);
end

I hope it helps!

Best Regards,

Abhishek Chakram