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

Question

Joydeep Bagchi 2023-9-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2025792-which-boundary-condition-is-correct-or-what-will-be-the-correct-alternative

编辑： Joydeep Bagchi 2024-2-27

采纳的回答： Torsten

在 MATLAB Online 中打开

% 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

% 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

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2023-9-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2025792-which-boundary-condition-is-correct-or-what-will-be-the-correct-alternative#answer_1318607

编辑：Torsten 2023-9-26

在 MATLAB Online 中打开

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 个评论
显示 7更早的评论隐藏 7更早的评论

Torsten 2023-9-26

编辑：Torsten 2023-9-26

在 MATLAB Online 中打开

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 2023-9-26

在 MATLAB Online 中打开

% 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

% 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

hold off

Joydeep Bagchi 2023-9-27

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 2023-9-27

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.

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

9 个评论
显示 7更早的评论隐藏 7更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

9 个评论 显示 7更早的评论隐藏 7更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

WeChat

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

9 个评论
显示 7更早的评论隐藏 7更早的评论