I am trying to solve a system of PDEs for a 1D flow ($Z$-axis) in time $t$ having several phases distinct by different boundary conditions of mixed type (their duration is fixed by defining a time for each phase). I am using FVM with TVD.
Problem description: my approach
Initially, I define initial conditions for t = 0, which I put in the resulting vector of unknown parameters. Next, I implement boundary conditions for the first phase for $Z = 0$ ($Z = 0.5$ for FVM) and $Z = N$ ($Z = N+0.5$), respectively, and execute the ODE solver ($\texttt{ode15s}$ as the set of PDEs is stiff).
I have noticed that during the time integration, the ode15s overwrites the boundary values for $Z = 0.5$, which I need to remain constant. Specifically, $\texttt{ode15s}$ returns a non-zero values for the time derivatives at $Z = 0.5$, thus "changing" the boundary values for each time layer.
Arguably, this is not correct as the $\texttt{ode15s}$ returns a non-zero values for the time derivatives on the $Z = 0.5$ boundary where they should be constant.
My solution: To "fix" this error, I have not found a better way than to zero (by "force") the time derivatives for $Z = 0.5$, after the $\texttt{ode15s}$ has been computed, to achieve constant values at this boundary throughout the simulation.
Principally, the algorithm should work by defining the time derivatives only once at the beginning and then changing only the boundary conditions, in my opinion.
Numerical scheme and implementation of ICs and BCs
- I use FVM with TVD (with ghost cells) and specify boundary conditions for the first and last finite volume walls ($Z = 0.5$ and $Z = N+0.5$, respectively).
- At the beginning of my algorithm, initial conditions are assigned to all central values of the finite volumes for time $t = 0$ for all $N$ volumes.
- At the beginning of each phase, the boundary conditions are again assigned to the first and last walls while zeroing of the time derivates (likely a mistake) is performed for the initial finite volume centre, i.e. $Z = 1$.
MATLAB example: my implementation I made a simple code document my problem and my implementation (the original code is "long"):
If I introduce ghost cells at the walls of the numerical grid, the function returns non-zero time derivatives at its beginning and end, with the $\texttt{ode15s}$ subsequently changing these values by time integration.
If I "hard" set the first and last time derivatives equal to 0, then these temperature values (first and last) are constant throughout the $\texttt{ode15s}$ time integration.
Questions
- What computation (implementation), different to mine, of the boundary conditions for the first time layer and for all time derivatives, to keep them constant all the time, would you use, rather than force these values to remain constant by introducing zero time derivatives ($\texttt{dTdt(1) = 0;}$ , $\texttt{dTdt(nz+1) = 0;}$)? (In other words, would you please have any recommendations on how to keep the boundary values constant in the ODE solver without overwriting the time derivatives within the run of the algorithm?).
% Simulation of a heat conduction
% 24012023
clear all
clc
% Input parameters
L = 1; % [m]
t = 1; % [s]
alfa = 0.1; % [degC/m]
Ta = 80; % [degC]
Tb = 20; % [degC]
% Mesh
nz = 20;
dz = L/nz; % [m]
% Initial conditions
T0 = Tb*ones(1,nz+1);
T0(1) = Ta; % Boundary condition/values for Z = 0
T0(nz+1) = Tb; % Boundary condition/values for Z = L
z = linspace(0,L,nz+1);
% ODE solver
opts = odeset('RelTol',1e-4,'AbsTol',1e-5);
% Ghost cell
[time,T_ODE] = ode15s(@heat_calc, [0 t], T0, opts, nz, alfa, dz, Ta, Tb, 0);
% dTdt = 0
[time2,T2_ODE] = ode15s(@heat_calc, [0 t], T0, opts, nz, alfa, dz, Ta, Tb, 1);
% Plotting results
subplot(2,1,1), plot(z,T_ODE),title('Ghost ODE')
subplot(2,1,2), plot(z,T2_ODE),title('dTdt=0 ODE')
%%
function dTdt = heat_calc(t, T, nz, alfa, dz, Ta, Tb, parametr)
% Add ghost cell
% Boundary condition for Z = 0 and Z = L, is such implementation correct?
Tg = [Ta; T; Tb];
dTdt = (alfa/dz^2).*(Tg(3:nz+3,1)...
- 2.*Tg(2:nz+2,1) + Tg(1:nz+1,1));
% For keeping Ta and Tb constant, is such implementation correct?
if parametr == 1
dTdt(1) = 0;
dTdt(nz+1) = 0;
end
end
In my opinion, the best way to ensure that the boundary conditions remain constant is to not make them dependent variables in the first place and only simulate the dynamics on the interior cells. The ODE function would be exactly the same except you would have to manually write out the cases for the leftmost and rightmost interior cells as their values would de generated by the finite volume stencil with the boundary condition values inside them. If you want to keep the exterior nodes in the state (this makes the bookkeeping easier on you, but you get these problem in exchange), then setting the derivatives to zero manually seems like the best course of action.
Edit: With minimal modifications to your code, something like this should work. I don't have MATLAB on this computer, so I haven't tested it. However, something like this should work. The important changes are in the grid and problem dimension. I changed the grid so that
nzis the number of internal nodes, then add two for the boundaries. The ODE function then only simulates the internal nodes. This is very important. Notice that the part of the ODE function that touches the boundaries still had to be manually specified, but it still matches the FV stencil. I then concatenate the output of the ODE to add the boundary conditions at both ends of the domain for plotting. The rows/columns might be slightly off here since I couldn't test but this should work.