Which boundary condition is correct or what will be the correct alternative?
7 views (last 30 days)
Show older comments
% 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?
0 Comments
Accepted Answer
Torsten
on 26 Sep 2023
Edited: Torsten
on 26 Sep 2023
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 ?
9 Comments
Torsten
on 27 Sep 2023
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.
More Answers (0)
See Also
Categories
Find more on Heat Transfer in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!