Which boundary condition is correct or what will be the correct alternative?

% Constants
Ly = 0.001;        % Thickness of the sheet (meters)
T = 15;             % Total simulation time (seconds)
Ny = 10;            % Number of spatial grid points
Nt = 300;           % Number of time steps
k1 = 0.15;           % Thermal conductivity of PVC W/m-K
rho = 1250;          % density in kg/m3
cp  = 1350 ;         % specific heat in J/kg-K
alpha = k1/(rho*cp);   % Thermal diffusivity (m^2/s)
T_ambient = 25;     % Ambient temperature (°C)
velocity_x = 0.1667; % Velocity in x-direction (m/s)

%contact with hot roller
k2 = 50;             %Thermal conductivity of steel roller W/m-K
h_roller = 500;      %convective heat transfer coefficient of roller W/m2-K
d_roller = 0.4;      %diameter of roller in m
contact_angle = 180; %in °
contact_length = pi*d_roller*contact_angle/360; %contact length with roller
T_roller = 100;      % Roller temperature

%Contact in air
distance_x = 0.5;   % Distance traveled in x-direction (meters)
h_air = 12;         % Convective heat transfer coefficient of air W/m2-K

%Contact in IR field
distance_ir = 0.3 ; % in m
radiative_flux = 50e3; %W/m2
absorption = 50; % in %
performance = 80; % in %
Net_radiative_intensity = radiative_flux*(absorption/100)*(performance/100);

%Contact with air again
distance_v = 1 ; % in m

%Naming the distance variables
y1 = contact_length;
y2 = contact_length + distance_x;
y3 = contact_length + distance_x+ distance_ir;
y4 = contact_length + distance_x+ distance_ir+distance_v;

% Discretization
dy = Ly / (Ny - 1);
dt = T / Nt;
y = linspace(0, Ly, Ny);
t = linspace(0, T, Nt);

% Initial temperature distribution
initial_temperature = 25;
distance_travelled = zeros(1,Nt); %initial distance travelled

% Initialize temperature matrix
u = zeros(Ny, Nt);
u(:, 1) = initial_temperature;

% Forward euler method
r = alpha * dt / (dy^2); % as 0 < r < 0.5

% Time-stepping loop (explicit method)
for n = 1:Nt - 1
    
    % Update position in x-direction
    distance_travelled(n+1) = distance_travelled(n) + velocity_x * dt;
    
    % Check if the total distance is within the contact length
    if distance_travelled(n) <= y1
        for i = 2:Ny - 1
            u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
        end
        
        % Apply boundary conditions for contact with the roller
        u(1, n + 1) = ((h_roller * T_roller * dy) + (k1 * u(2, n + 1))) / (k1 + h_roller * dy);
        u(Ny, n + 1) = ((h_air * T_ambient * dy) + (k1 * u(Ny - 1, n + 1))) / (k1 + h_air * dy);
        
        % For the next part in contact with air 
    elseif distance_travelled(n) <= y2
        for i = 2:Ny - 1
            u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
        end
        % Apply boundary conditions for contact with air
        u(1, n + 1) = ((h_air * T_ambient * dy) + (k1 * u(2, n + 1))) / (k1 + h_air * dy);
        u(Ny, n + 1) = ((h_air * T_ambient * dy) + (k1 * u(Ny - 1, n + 1))) / (k1 + h_air * dy);
        
        %In contact with IR field radiation
    elseif distance_travelled(n) <= y3
        for i = 2:Ny - 1
            u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
        end
        % Apply boundary conditions for contact with IR field
        u(1, n + 1) = ((h_air * T_ambient * dy) + (k1 * u(2, n + 1))) / (k1 + h_air * dy);
        u(Ny, n + 1) = u(Ny - 1, n + 1) + ((Net_radiative_intensity * dy) / k1);
        
        % Distance in contact with air again
    else
        for i = 2:Ny - 1
            u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
        end
        % Apply boundary conditions for contact with air
        u(1, n + 1) = ((h_air * T_ambient * dy) + (k1 * u(2, n + 1))) / (k1 + h_air * dy);
        u(Ny, n + 1) = ((h_air * T_ambient * dy) + (k1 * u(Ny - 1, n + 1))) / (k1 + h_air * dy);
        
        
        
    end    
    
    % Break the loop when the total distance is reached
    if distance_travelled(n + 1) >= y4
        break;
    end
end

% Plot the temperature distribution over time
figure;
for n = 1:min(Nt,numel(distance_travelled))
    plot(y, u(:, n));
    xlabel('Thickness (m)');
    ylabel('Temperature (°C)');
    title(['Time = ', num2str(t(n)), ' seconds']);
    axis([0 Ly 20 150]); % Adjust the axis limits if needed
    pause(0.1); % Pause to create animation effect
    
    % Break the plotting loop when the loop breaks
    if distance_travelled(n) >= y4
        break;
    end
end

which boundary condition is correct? The PVC material's 1st element is in contact with the roller whose temperature is maintained at 100°C so ideally the mode of heat transfer should be conduction. In that case, either dirichlet boundary condition or neumann boundary condition should be implemented

If i apply dirichlet boundary condition u(1,n+1) = T_roller so it is fixing the temperature to 100°C until it is contact with roller which will not be the case . In reality, the temperature of the first element will rise more and try to reach 100°C as compared to other elements but won't be 100°C

For Neumann boundary condition, I need the information of temperature gradient or flux but unfortunately, I don't have neither of that.

I am not sure the boundary condition which I have applied here is correct or not? Can someone please clarify this or help me in this issue?

Answer 1

[Torsten]

If you can estimate the heat transfer coefficient k between PVC material and roller, you could use

-lambda_PVC * dT_PVC/dn = k * (T_PVD - T_roller)

where d/dn is the spatial derivative in the direction of the outer normal. Look up "Robin boundary condition".

If k is large, this means that the PVC temperature follows the roller temperature very fast. In this case, it is justified to use the Dirichlet boundary condition for T_PVC.

But if the two materials are in direct contact, they should share the same temperature at this location, shouldn't they ?

[Joydeep Bagchi]

Thankyou for the reply

Could you please specify how to compute dT_PVC/dn because the temperatures are not changing if I do not give any boundary condition so if you implement the code (u(2,n+1) - u(1,n+1))/dy will give 0 as result

And if I use Dirichlet condition which will be u(1,n+1) = 100°C so it will stay at 100 °C for the time in contact with roller which is not possible , rather it will rapidly increase to 100 °C and it will decrease gradually not in a steep manner

[Torsten]

I don't understand your code. What is the partial differential equation you are trying to solve ? What does the distance_travelled array mean ?

[Joydeep Bagchi]

I am trying to solve diffusion equation and distance_travelled is that the element or the geometry is moving

basically, I have to simulate the transient temperature distribution of a web across the thickness so I assume it as a geometry consisting of various elements and it is exposed to different environment as it is written in comments.

[Torsten]

Ok.

Is the following plot as expected if the contact to the environment would not change ?

I think it's not advisable to change the boundary conditions so drastically as you do. For each change, you should use a "numerical transition zone".

% Constants
Ly = 0.001;        % Thickness of the sheet (meters)
T = 15;             % Total simulation time (seconds)
Ny = 30;            % Number of spatial grid points
Nt = 3000;           % Number of time steps
k1 = 0.15;           % Thermal conductivity of PVC W/m-K
rho = 1250;          % density in kg/m3
cp  = 1350 ;         % specific heat in J/kg-K
alpha = k1/(rho*cp);   % Thermal diffusivity (m^2/s)
T_ambient = 25;     % Ambient temperature (°C)
velocity_x = 0.1667; % Velocity in x-direction (m/s)
%contact with hot roller
k2 = 50;             %Thermal conductivity of steel roller W/m-K
h_roller = 500;      %convective heat transfer coefficient of roller W/m2-K
d_roller = 0.4;      %diameter of roller in m
contact_angle = 180; %in °
contact_length = pi*d_roller*contact_angle/360; %contact length with roller
T_roller = 100;      % Roller temperature
%Contact in air
distance_x = 0.5;   % Distance traveled in x-direction (meters)
h_air = 12;         % Convective heat transfer coefficient of air W/m2-K
%Contact in IR field
distance_ir = 0.3 ; % in m
radiative_flux = 50e3; %W/m2
absorption = 50; % in %
performance = 80; % in %
Net_radiative_intensity = radiative_flux*(absorption/100)*(performance/100);
%Contact with air again
distance_v = 1 ; % in m
%Naming the distance variables
y1 = contact_length;
y2 = contact_length + distance_x;
y3 = contact_length + distance_x+ distance_ir;
y4 = contact_length + distance_x+ distance_ir+distance_v;
% Discretization
dy = Ly / (Ny - 1);
dt = T / Nt;
y = linspace(0, Ly, Ny);
t = linspace(0, T, Nt);
% Initial temperature distribution
initial_temperature = 25;
distance_travelled = zeros(1,Nt); %initial distance travelled
% Initialize temperature matrix
u = zeros(Ny, Nt);
u(:, 1) = initial_temperature;
% Forward euler method
r = alpha * dt / (dy^2); % as 0 < r < 0.5
% Time-stepping loop (explicit method)
for n = 1:Nt - 1
    
    % Update position in x-direction
    %distance_travelled(n+1) = distance_travelled(n) + velocity_x * dt;
    
    % Check if the total distance is within the contact length
    %if distance_travelled(n) <= y1
        for i = 2:Ny - 1
            u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
        end
            % Apply boundary conditions for contact with the roller
        %u(1, n + 1) = ((h_roller * T_roller * dy) + (k1 * u(2, n))) / (k1 + h_roller * dy);
        %u(Ny, n + 1) = ((h_air * T_ambient * dy) + (k1 * u(Ny - 1, n))) / (k1 + h_air * dy);
        u(1, n + 1) = u(1,n) + dt*(alpha*(u(2,n)-u(1,n))/dy - h_roller/(rho*cp)*(u(1,n)-T_roller))/(dy/2);
        u(Ny, n + 1) = u(Ny,n) + dt*(-h_air/(rho*cp)*(u(Ny, n)-T_ambient)-alpha*(u(Ny,n)-u(Ny-1,n))/dy)/(dy/2);
end
plot(y,u(:,1:100:Nt))     

[Torsten]

% Constants
Ly = 0.001;        % Thickness of the sheet (meters)
T = 15;             % Total simulation time (seconds)
Ny = 10;            % Number of spatial grid points
Nt = 3000;           % Number of time steps
k1 = 0.15;           % Thermal conductivity of PVC W/m-K
rho = 1250;          % density in kg/m3
cp  = 1350 ;         % specific heat in J/kg-K
alpha = k1/(rho*cp);   % Thermal diffusivity (m^2/s)
T_ambient = 25;     % Ambient temperature (°C)
velocity_x = 0.1667; % Velocity in x-direction (m/s)

%contact with hot roller
k2 = 50;             %Thermal conductivity of steel roller W/m-K
h_roller = 500;      %convective heat transfer coefficient of roller W/m2-K
d_roller = 0.4;      %diameter of roller in m
contact_angle = 180; %in °
contact_length = pi*d_roller*contact_angle/360; %contact length with roller
T_roller = 100;      % Roller temperature

%Contact in air
distance_x = 0.5;   % Distance traveled in x-direction (meters)
h_air = 12;         % Convective heat transfer coefficient of air W/m2-K

%Contact in IR field
distance_ir = 0.3 ; % in m
radiative_flux = 50e3; %W/m2
absorption = 50; % in %
performance = 80; % in %
Net_radiative_intensity = radiative_flux*(absorption/100)*(performance/100);

%Contact with air again
distance_v = 1 ; % in m

%Naming the distance variables
y1 = contact_length;
y2 = contact_length + distance_x;
y3 = contact_length + distance_x+ distance_ir;
y4 = contact_length + distance_x+ distance_ir+distance_v;

% Discretization
dy = Ly / (Ny - 1);
dt = T / Nt;
y = linspace(0, Ly, Ny);
t = linspace(0, T, Nt);

% Initial temperature distribution
initial_temperature = 25;
distance_travelled = zeros(1,Nt); %initial distance travelled

% Initialize temperature matrix
u = zeros(Ny, Nt);
u(:, 1) = initial_temperature;

% Forward euler method
r = alpha * dt / (dy^2); % as 0 < r < 0.5

% Time-stepping loop (explicit method)
for n = 1:Nt - 1
    
    % Update position in x-direction
    distance_travelled(n+1) = distance_travelled(n) + velocity_x * dt;
    
    % Check if the total distance is within the contact length
    if distance_travelled(n) <= y1
        for i = 2:Ny - 1
            u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
        end
            % Apply boundary conditions for contact with the roller
        %u(1, n + 1) = ((h_roller * T_roller * dy) + (k1 * u(2, n))) / (k1 + h_roller * dy);
        %u(Ny, n + 1) = ((h_air * T_ambient * dy) + (k1 * u(Ny - 1, n))) / (k1 + h_air * dy);
        u(1, n + 1) = u(1,n) + dt*(alpha*(u(2,n)-u(1,n))/dy - h_roller/(rho*cp)*(u(1,n)-T_roller))/(dy/2);
        u(Ny, n + 1) = u(Ny,n) + dt*(-h_air/(rho*cp)*(u(Ny, n)-T_ambient)-alpha*(u(Ny,n)-u(Ny-1,n))/dy)/(dy/2);
        
        % For the next part in contact with air 
    elseif distance_travelled(n) <= y2
        for i = 2:Ny - 1
            u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
        end
        % Apply boundary conditions for contact with air
        %u(1, n + 1) = ((h_air * T_ambient * dy) + (k1 * u(2, n + 1))) / (k1 + h_air * dy);
        %u(Ny, n + 1) = ((h_air * T_ambient * dy) + (k1 * u(Ny - 1, n + 1))) / (k1 + h_air * dy);
        u(1, n + 1) = u(1,n) + dt*(alpha*(u(2,n)-u(1,n))/dy - h_air/(rho*cp)*(u(1,n)-T_ambient))/(dy/2);
        u(Ny, n + 1) = u(Ny,n) + dt*(-h_air/(rho*cp)*(u(Ny, n)-T_ambient)-alpha*(u(Ny,n)-u(Ny-1,n))/dy)/(dy/2);        
        %In contact with IR field radiation
    elseif distance_travelled(n) <= y3
        for i = 2:Ny - 1
            u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
        end
        % Apply boundary conditions for contact with IR field
        %u(1, n + 1) = ((h_air * T_ambient * dy) + (k1 * u(2, n + 1))) / (k1 + h_air * dy);
        %u(Ny, n + 1) = u(Ny - 1, n + 1) + ((Net_radiative_intensity * dy) / k1);
        u(1, n + 1) = u(1,n) + dt*(alpha*(u(2,n)-u(1,n))/dy - h_air/(rho*cp)*(u(1,n)-T_ambient))/(dy/2);
        u(Ny, n + 1) = u(Ny,n) + dt*(Net_radiative_intensity/(rho*cp)-alpha*(u(Ny,n)-u(Ny-1,n))/dy)/(dy/2);
        
        % Distance in contact with air again
    else
        for i = 2:Ny - 1
            u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
        end
        % Apply boundary conditions for contact with air
        %u(1, n + 1) = ((h_air * T_ambient * dy) + (k1 * u(2, n + 1))) / (k1 + h_air * dy);
        %u(Ny, n + 1) = ((h_air * T_ambient * dy) + (k1 * u(Ny - 1, n + 1))) / (k1 + h_air * dy);
        u(1, n + 1) = u(1,n) + dt*(alpha*(u(2,n)-u(1,n))/dy - h_air/(rho*cp)*(u(1,n)-T_ambient))/(dy/2);
        u(Ny, n + 1) = u(Ny,n) + dt*(-h_air/(rho*cp)*(u(Ny, n)-T_ambient)-alpha*(u(Ny,n)-u(Ny-1,n))/dy)/(dy/2);          
        
        
    end    
    
    % Break the loop when the total distance is reached
    if distance_travelled(n + 1) >= y4
        break;
    end
end

% Plot the temperature distribution over time
hold on
for n = 1:100:min(Nt,numel(distance_travelled))
    if distance_travelled(n) <= y1
      plot(y, u(:, n),'k');
    elseif distance_travelled(n) <= y2
      plot(y, u(:, n),'b');
    elseif distance_travelled(n) <= y3
      plot(y, u(:, n),'r');   
    elseif distance_travelled(n) <= y4    
      plot(y, u(:, n),'y');          
    end
    xlabel('Thickness (m)');
    ylabel('Temperature (°C)');
    title(['Time = ', num2str(t(n)), ' seconds']);
    axis([0 Ly 20 150]); % Adjust the axis limits if needed
    pause(0.1); % Pause to create animation effect
    
    % Break the plotting loop when the loop breaks
    if distance_travelled(n) >= y4
        break;
    end
end
hold off

[Joydeep Bagchi]

Thankyou for the reply, but unfortunately the plot is not the expected one because we don't need the transition zone however I saw some changes which I will implement in the code

But I think the main question is still there

However I found one link as a reference to the master thesis where the boundary condition for the roller how I have used is mentioned here

pg 31

https://shareok.org/bitstream/handle/11244/45277/Lu_okstate_0664D_14232.pdf?sequence=1

[Torsten]

Thankyou for the reply, but unfortunately the plot is not the expected one because we don't need the transition zone

There is no transition zone used in the last code I posted.

But I think the main question is still there

Which question ? You mean the boundary condition you have to use on both ends of the web ?

However I found one link as a reference to the master thesis where the boundary condition for the roller how I have used is mentioned here

Your boundary condition terms use values for u at time t(n+1) on the right-hand sides. Thus you have a mixture of explicit and implicit Euler for which I'm not able to predict whether it will work out correctly.

[Joydeep Bagchi]

Once agian thanks Torsten for the reply.

There is no transition zone used in the last code I posted.

Ok, but I need a continuous transition not like the temperature distribution at each time step rather a continuous temperature distribution through the whole time span.

Which question ? You mean the boundary condition you have to use on both ends of the web ?

The link of the pdf which I have posted has the answer in it so I have implemented that condition which is right.

Your boundary condition terms use values for u at time t(n+1) on the right-hand sides. Thus you have a mixture of explicit and implicit Euler for which I'm not able to predict whether it will work out correctly.

I need to give the term for gradient thatswhy I have written like that which is working fine for me.

Thanks a lot for your effort and time.

[Torsten]

Ok, but I need a continuous transition not like the temperature distribution at each time step rather a continuous temperature distribution through the whole time span.

The temperature distribution over the web will be continuous if you choose more curves to plot.

The link of the pdf which I have posted has the answer in it so I have implemented that condition which is right.

I cannot find a boundary condition for the roller on page 31 of the document.

I need to give the term for gradient thatswhy I have written like that which is working fine for me.

It's your code, but your advancement in time is not Explicit Euler.