Implementing infinity like boundary condition for 1D diffusion equation solved with implicit finite difference method

Viewed 638

I am using following MATLAB code for implementing 1D diffusion equation along a rod with implicit finite difference method. `

xsize   =   10;                          % Model size, m
xnum    =   10;                          % Number of nodes
xstp    =   xsize/(xnum-1);              % Grid step
tnum    =   504;                         % number of timesteps
kappa   =   833.33;                      % Thermal diffusivity, m^2/s
dt     =   300;                          % Timestep
x       =   0:xstp:xsize;                %Creating vector for nodal point positions

tlbc = sin(linspace(0.1,2.9,tnum));      % left boundary condition

%Define initial temperature profile
tback   =   0;                           % background temperature, K

for i=1:1:xnum
    % Background profile  
    t0imp(i)    =   tback;               % profile for implicit solving
end
% Time cycle
timesum=0; % Elapsed time
for t=1:1:tnum

    % Matrix of coefficients initialization for implicit solving
    L       =       sparse(xnum,xnum);
    % Vector of right part initialization for implicit solving
    R       =       zeros(xnum,1);

    % Implicit solving of 1D temperature equation: dT/dt=kappa*d2T/dx2
    % Composing matrix of coefficients L()
    % and vector (column) of right parts R()
    % First point: T=tback
    L(1,1)  =       1;
    R(1,1)  =       tlbc(t);
    % Intermediate points
    for i=2:1:xnum-1
        % dT/dt=kappa*d2T/dx2
        % Tnew(i)/dt-kappa*(Tnew(i-1)-2*Tnew(i)+Tnew(i+1))/dx^2=Told(i)/dt

        L(i,i-1)    =   -kappa/xstp^2;
        L(i,i)      =   1/dt+2*kappa/xstp^2;
        L(i,i+1)    =   -kappa/xstp^2;
        R(i,1)      =   t0imp(i)/dt;
    end
    % Last point:T=tback
    L(xnum,xnum)    =   1;

    R(xnum,1)       =  tback;
    % Obtaining solution for implicit temperature profile
    t1imp           =   L\R;

end `

I want to change right boundary condition to infinite boundary condition so that presence of boundary at 10 m has no effect on heat flow and temperature flows in 10 m domain as if there was no boundary at 10 m.

1 Answers

If you have localized phenomenon that is far away from the boundary then I suggest you use a non-uniform grid for your x-space discretization. The mesh should have more points in the region of interest. You can then choose a large xsize but keep the number of mesh points reasonable.

Non-uniform grids are trickier to implement (not that bad for 1-d) but fortunately Matlab has pdepe that solves non-uniform equations of the type

\begin{equation} c\left(x , t,u, \frac{\partial u}{\partial x} \right) \frac{\partial u}{\partial t} = x^{-m} \frac{\partial}{\partial t} \left( x^m f\left(x , t,u, \frac{\partial u}{\partial x} \right)\right) + S\left(x , t,u, \frac{\partial u}{\partial x} \right) \end{equation}

with user defined meshes. You can read more about pdepe here

https://www.mathworks.com/help/matlab/ref/pdepe.html

They also solve the equation you are interested in.

Related