% ------------------------------------------------------------------------
% Corso di Termofluidodinamica Computazionale (CFD)
% Homework No. 2, AA 2018/2019
% 
% Problem: 2D steady-state heat conduction
%                T_xx + T_yy = 0
%
% Domain and boundary conditions (BC): [0,L]x[-t/2,t/2] (rectangular fin)
%
%             y ^
%               |
%               |
%               |    Convection (h, Ti)
%           t/2 +--------------------------+
%               |                          |
%               |                          |
% Dirichlet (Tb)| -.-.-.-.- symm -.-.-.-.- | Convection (h, Ti) ------> x
%               |                          |
%               |  thermal conductivity k  |
%          -t/2 +--------------------------+
%               0    Convection (h, Ti)    L
%
% Solution methodology: semi-analytical
%
% T-Ti = sum (i=1 to Nq) [ Ai*exp(qi*x) + Bi*exp(-qi*x) ] * cos(qi*y)
%   Ai, Bi, qi = ?
%   The term [Ai*exp(qi*x)+Bi*exp(-qi*x)]*cos(qi*y) is called mode
%   Convective BCs are made valid for each mode (linearity)
%   Dirichlet  BC @ fin base is global, i.e., it is satisfied by modes sum
%
% Date: 22/3/2019
% Author: R. Zamolo
% ------------------------------------------------------------------------

% Dati
k  = 50    ; % W/(m K)
h  = 500   ; % W/(m^2 K)
t  = 0.100 ; % m
L  = 0.200 ; % m
Tb = 200   ; % �C
Ti = 25    ; % �C

% Number of modes (1 - 1000)
Nq = 100 ; % hint: try to investigate the effect of Nq on the results
%                  (heat flux for example) for different aspect ratios L/t

% Convective BC @ upper surface (y=t/2)
% ( BC @ lower surface y=-t/2 automatically satisfied for symmetry 
%   of solution along y: cos(qi*y) )
%                     -k * dT/dy = h * (T-Ti)
%
% =>  zi * tan( zi ) = b
%     where:
%     zi = qi * t / 2
%     b  = h * t / ( 2 * k )
%     =>  we have to find Nq roots zi
b = h * t / ( 2 * k ) ;
z = pi * ( 0 : Nq-1 )' ; % initial guess for zi (roots of tan(zi)=0)
z = z + atan( b ./ z ) ; % improvement of initial guess
z( 1 ) = atan( b ) ;     % improvement of first zero
% root-finding by Newton's method (step = -function/derivative)
dz = 1 ;
while max( abs( dz ./ z ) ) > eps % convergence to machine-zero
 tg = tan( z ) ;
 d  = tg + z .* ( 1 + tg .^ 2 ) ; % derivative
 dz = -( z .* tg - b ) ./ d  ;    % step = -function / derivative
 z = z + dz ;
end
q = z * 2 / t ; % qi = zi * 2 / t

% Dirichlet BC @ fin base (x=0), T = Tb
% =>  sum (i=1 to Nq) Ci * cos(qi*y) = DT
%     where:
%     Ci = Ai + Bi 
%     DT = Tb - Ti
DT = Tb - Ti ;
Ny = 4 * Nq ;                     % # of points for least squares
y = linspace( 0 , t/2 , Ny )' ;   % y vector
M = cos( y * q' ) ;               % matrix for least squares
C = M \ ( DT * ones( Ny , 1 ) ) ; % Ci solved by least squares

% Convective BC @ right surface (x=L)
%                     -k * dT/dx = h * (T-Ti)
%
% =>  alf * Ai + bet * Bi = 0
%     where:
%     alf = exp( qi*L) * ( 1 + k*qi/h )
%     bet = exp(-qi*L) * ( 1 - k*qi/h )
alf = exp(  q * L ) .* ( 1 + k * q / h ) ;
bet = exp( -q * L ) .* ( 1 - k * q / h ) ;

% So we can calculate Ai and Bi by solving a 2x2 linear system
%       Ai +       Bi = Ci <= we know Ci
% alf * Ai + bet * Bi = 0
%
% =>  Ai = bet * Ci / ( bet - alf )
% =>  Bi = alf * Ci / ( alf - bet )
A = bet .* C ./ ( bet - alf ) ;
B = alf .* C ./ ( alf - bet ) ;
%   but we must be careful with numerics:
%   - alf grows exponentially with Nq => overflow
%   - bet decreases exponentially with Nq
%   if alf approaches infinity (overflow) (bet approaches 0)
%   => Ai --> 0
%   => Bi --> Ci
A( isinf(alf) ) = 0 ;
B( isinf(alf) ) = C( isinf(alf) ) ;
% (modes associated to overflowed alf CAN NOT be neglected!!)

% Let's evaluate the unit heat flux qf [W/m]
%   qf = integral (y from -t/2 to t/2, x=0) (-k * dT/dx) * dy
%      = 2 * k * sum (i=1 to Nq) ( Bi - Ai ) * sin(qi*t/2)
qf = 2 * k * sum( ( B - A ) .* sin( q * t / 2 ) ) ;
fprintf( 'Semi-analytical heat flux = %3.2f W/m\n' , qf ) ;

% Let's evaluate the full temperature field
x = linspace(    0 , L   , round(L*500) )' ;
y = linspace( -t/2 , t/2 , round(t*500) ) ;
T = 0 * x * y ;
%   to avoid overflows, the solution is explicitly written in terms of
%   alf and bet and simplified when alf(i) approaches infinity (overflow)
for i = 1 : Nq
 if isfinite( alf(i) )
  T = T + C(i)*( exp(q(i)*(x-  L)) * ( 1 - k * q(i) / h ) - exp(q(i)*(L-x)) * ( 1 + k * q(i) / h ) ) * cos(q(i)*y) / (bet(i)-alf(i)) ;
 else
  T = T - C(i)*( exp(q(i)*(x-2*L)) * ( 1 - k * q(i) / h ) / ( 1 + k * q(i) / h ) - exp(q(i)*(-x)) ) * cos(q(i)*y) ;
 end
end
T = T + Ti ;

% 3D temperature plot
surf( x , y , T' ) ;