PDE solve with numerical way

Hello, I need some help with building a code for numerical solve.

The pictures above is about the partial differential equation (PDE) and boundary conditions for the dimensionless temperature (θ_n1, θ_n2, θ_w) along the z and y directions. I have modified the terms (z, y1, y2) to del_y1, del_y2, and del_z. My goal is to compare the numerical solutions with the approximate ones by LMTD method, and assess the accuracy.

clear all;
clc;
clear;

m_0=1;
k=1;

y1 = 1;
y2 = 1;
z_L = 1;

del_y1 = 1e-3;
del_y2 = 1e-3;
del_z = 1e-3;

%% Set up matrix
m = z_L/del_z+1; % x-direc
n1=y1/del_y1+1; % y1-direc
n2=y2/del_y2+1; % y2-direc

theta_n1 = ones(m,n1); %m*n1 with all '1' element
theta_n2 = ones(m,n2); %m*n2 with all '1' element
theta_w = ones(m,1); %m*1 array with all '1' element


%%Initial theta(BC)
theta_1_in = 0; % Inlet_n1 : 0
theta_2_in = 1; % Inlet_n1 : 1

theta_n1 = theta_n1*0.5; % theta_n1 : 0.5
theta_n2 = theta_n2*0.5; % theta_n2 : 0.5
theta_w = theta_w*0.5; % theta_w : 0.5


%%Tolerance
errormax=1e-5;
error=1;
errors = [];
iteration=0;

while (error>errormax)
    theta_n1_old=theta_n1;
    theta_n2_old=theta_n2;
    theta_w_old= theta_w;

%% Inlet + wall initiate
    theta_n1(1,1:n1)=theta_1_in; 
    theta_n2(m,1:n2)=theta_2_in;

    theta_w(1) = 0; % zeta = 0 -> 0
    theta_w(m) = 1; % zeta = zeta_L -> 1
%% B.Cs
%% Top & Bottom
    for M=2:m-1
        % condition in contact line
        theta_w(M)= (k /(1+k)) *  theta_n2(M,n2-1) + (1/1+k) * theta_n1(M,n1-1);
        theta_n1(M,n1)= theta_w(M);
        theta_n2(M,n2)= theta_w(M);
        % condition in outer line
        theta_n1(M,1)= theta_n1(M,2);
        theta_n2(M,1)= theta_n1(M,2);
    end

%% Main + Output
    for N=2:n1-1
        C1 = (1/(del_y1*del_y1) - 0.75*(1-(N*del_y1)^2)/del_z);
        for M=2:m-1
            theta_n1(M,N)= (1/C1)*((theta_n1(M,N+1)+theta_n1(M,N-1))/(2*del_y1*del_y1) - (0.75*(1-(del_y1*N)^2)*theta_n1(M+1,N))/del_z);
        end
        theta_n1(m,N) = ((theta_n1(m-1,N+1)+theta_n1(m-1,N-1))/(2*del_y1*del_y1)-theta_n1(m-1,N)*C1)*del_z/(0.75*(1-(del_y1*N)^2));
    end

    for N=2:n2-1
        C2 = (1/(del_y2*del_y2) + 0.75*m_0*(1-(N*del_y2)^2)/del_z);
        for M=2:m-1
            theta_n2(M,N)= (1/C2)*((theta_n2(M,N+1)+theta_n2(M,N-1))/(2*del_y2*del_y2) + (0.75*m_0*(1-(del_y2*N)^2)*theta_n2(M+1,N))/del_z);
        end
        theta_n2(m,N) = (-(theta_n2(m-1,N+1)+theta_n2(m-1,N-1))/(2*del_y2*del_y2)+theta_n2(m-1,N)*C2)*del_z/(0.75*m_0*(1-(del_y2*N)^2));
    end
    theta_n1(m,1)= theta_n1(m,2);
    theta_n2(m,1)= theta_n2(m,2);
    if iteration ==5
        contourf(theta_n1);
    end

%% iteration
    error=max(max(sqrt(((theta_n1-theta_n1_old).^2 + (theta_n2-theta_n2_old).^2 + (theta_w-theta_w_old).^2)/3)));
    iteration=iteration +1;
    errors(iteration) = error;
end

The solutions do not converge; instead, they diverge. Simulation results(theta_n1, theta_n2, theta_w) should be in the range of 0 to 1. I have checked the code multiple times, but I am still struggling with this problems.

Are there any issues with the code? Or is there a mathematical flaw inherent in the difference equation?

It would be helpful if someone helps me to solve this problem.

Thank you.

Ref 1)

Vera, M., & Linan, A. (2010). Laminar counterflow parallel-plate heat exchangers: exact and approximate solutions. International Journal of Heat and Mass Transfer, 53(21-22), 4885-4898.

Ref 2)

How I change PDE into difference equation.

+)

I think there may be some errors in the '%% Main + Output' section.

Theta_n1 and Theta_n2 share the same axis (z), and the given boundary conditions are at z = 0 and z = 1. If everything goes well, I will upload the modified version.

++)

I found some mistakes with codes, but it still doesn't work.

Here are the codes that I've changed.

clear all;
clc;
clear;

m_0=1;
k=1;

y1 = 1;
y2 = 1;
z_L = 1;

del_y1 = 1e-3;
del_y2 = 1e-3;
del_z = 1e-3;

%% Set up matrix
m = z_L/del_z+1; % x-direc
n1=y1/del_y1+1; % y1-direc
n2=y2/del_y2+1; % y2-direc

theta_n1 = ones(m,n1); %m*n1 with all '1' element
theta_n2 = ones(m,n2); %m*n2 with all '1' element
theta_w = ones(m,1); %m*1 array with all '1' element

%%Initial theta(BC)
theta_1_in = 0; % Inlet_n1 : 0
theta_2_in = 1; % Inlet_n1 : 1

theta_n1 = theta_n1*0.5; % 초기 설정_n1 : 0.5
theta_n2 = theta_n2*0.5; % 초기 설정_n2 : 0.5
theta_w = theta_w*0.5; % 초기 설정_w : 0.5

%%Tolerance
errormax=1e-5;
error=1;
errors = [];
iteration=0;

while (error>errormax)
    theta_n1_old=theta_n1;
    theta_n2_old=theta_n2;
    theta_w_old= theta_w;

%% Inlet + wall initiate
    theta_n1(1,:)=theta_1_in; 
    theta_n2(m,:)=theta_2_in;

    theta_w(1) = 0; % zeta = 0 -> 0
    theta_w(m) = 1; % zeta = zeta_L -> 1

    theta_n1(m,n1)= 1; 
    theta_n2(1,n2)= 0;
    
%% B.Cs

%% Main + Output
    for N=2:n1-1
        a = del_y1*del_y1;
        b = 0.75*(1-(N*del_y1)^2);
        c = del_z;
        for M=2:m
            theta_n1(M,N) = (c/(2*a*b))*(theta_n1(M-1,N+1)+theta_n1(M-1,N-1))-((c/(a*b))-1)*(theta_n1(M-1,N));
        end        
    end

    for N=2:n2-1
        a = del_y2*del_y2;
        b = -0.75*m_0*(1-(N*del_y2)^2);
        c = del_z;
        for M=1:m-1
            theta_n2(m-M,N) = -(c/(2*a*b))*(theta_n2(m-M+1,N+1)+theta_n2(m-M+1,N-1))+((c/(a*b))+1)*(theta_n2(m-M+1,N));
        end
    end
%% Top & Bottom
    for M=2:m-1
        % condition in contact line
        theta_w(M)= (k /(1+k)) *  theta_n2(M,n2-1) + (1/1+k) * theta_n1(M,n1-1);
        theta_n1(M,n1)= theta_w(M);
        theta_n2(M,n2)= theta_w(M);
        % condition in outer line
        theta_n1(M,1)= theta_n1(M,2);
        theta_n2(m-M+1,1)= theta_n2(m-M+1,2);
    end
    theta_n1(m,1)= theta_n1(m,2);
    theta_n2(1,1)= theta_n2(1,2);

%% iteration
    error=max(max(abs(theta_n1-theta_n1_old)));
    iteration=iteration +1;
    errors(iteration) = error;
end

[Torsten]

Are theta_w(xi) and nu(xi) externally given or are they to be computed due to a coupling of the temperature fields at the position of the thermally thin wall ?

[JUN-HO PARK]

They are not given. They should be computed with coupling conditions.

[Torsten]

So the main task is to determine the temperature profile for theta on y1 = y2 = 1 that goes from theta = 0 at xi = 0 up to theta = 1 at xi = xi_L .

Assume a temperature profile theta at the interface.

Then you can integrate both partial differential equations for theta1 and theta2 separately using MATLAB's "pdepe".

Usually, the heat flux condition

dtheta1/dy1 = -k*dtheta2/dy2 

at the interface will be violated for the two solutions theta1(xi,y1) and theta2(xi,y2).

Define the error equations as

dtheta1/dy1 + k*dtheta2/dy2 = 0

with the temperatures theta at the interface as unknowns and use "fsolve" to adjust them such that the heat flux condition holds.

[JUN-HO PARK]

I'll try it with "pdepe" ftn. Thank you for helping me.

[Torsten]

As said, you will have to couple "fsolve" with "pdepe".

The unknowns for "fsolve" are the temperatures at the interface. Their number depends on the number of xi points you choose in "pdepe". xi corresponds to t in "pdepe", yi corresponds to x.

Within the function for fsolve where you define your equations dtheta1/dy1 + k*dtheta2/dy2 = 0, you call "pdepe" twice: once for both counterflow equations.

In the definition of the boundary condition for "pdepe" at yi = 1, you should use the "interp1" function to interpolate the assumed temperature profile at the interface to the value of xi needed by "pdepe" during the integration.

To evaluate dtheta1/dy1 and dtheta2/dy2 for the error equation dtheta1/dy1 + k*dtheta2/dy2 = 0, use "pdeval":

https://uk.mathworks.com/help/matlab/ref/pdeval.html