Implementing symmetry boundary condition for the diffusion equation

[Atr cheema]

The following lines of codes implements 1D diffusion equation on 10 m long rod with fixed temperature at right boundary and right boundary temperature varying with time.

    xsize   =   10;               % Model size, m
xnum    =   10;                 % Number of nodes
xstp    =   xsize/(xnum-1);     % Grid step
tnum    =   504;                 % number of timesteps
kappa   =   833.33;            %     k/rhocp;            % 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
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
        %%% Eq 10.6
        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

My question is how can I implement a symmetry boundary condition. I mean the temperature calculated for 9th node is assigned to 10th node? Also I want to know how can I implement infinite boundary condition?

 

[NFuller]

I mean the temperature calculated for 9th node is assigned to 10th node?

Node10=Node9 ...?
I don't understand why you would want to do this though. Can you explain how this relates to the diffusion problem stated?

Also I want to know how can I implement infinite boundary condition?

For numerical solvers, you can't do this in general. I'm also not aware of any physical problem which would require an infinite boundary condition. If you want an infinite boundary for whatever reason, the best you can do numerically is assign a very large value at the needed boundary.

 

[Atr cheema]

Node10=Node9 ...?
I don't understand why you would want to do this though. Can you explain how this relates to the diffusion problem stated?

I want to implement a boundary condition where presence of boundary has no impact . I am using diffusion equation to model fluid flow in porous medium and want to correlate the numerical model with an analytical model where flux at right boundary is zero i.e. $\frac{\partial u}{\partial x} = 0$ . Consider u as fluid head instead of temperature. I think in this case, having zero flux at right boundary means that right boundary is having no impact on fluid head/heat flow in x direction.

 

[NFuller]

I think in this case, having zero flux at right boundary means that right boundary is having no impact on fluid head/heat flow in x direction.

Zero flux is equivalent to one end having an impermeable barrier, and would definitely have an impact on the solution. If your goal is to model the absence of a boundary or a boundary at infinity, then your boundary condition can be taken as the diffusion equation at the end point,
$$\frac{\partial u}{\partial t}=D\frac{\partial^{2}u}{\partial x^{2}}\bigg |_{x=10}$$
Depending on how your algorithm works, this can be done with a one-sided finite difference at the boundary.

 

[Atr cheema]

Depending on how your algorithm works, this can be done with a one-sided finite difference at the boundary.

Could you please explain this that in the given case, how can this infinity boundary condition be implemented here?