tcn_std.m

function q=tcn_std(l1,l2,T,r0,n)
%TCN_STD - Calculate the steady state discharge for a well located between a constant head and a no flow boundaries
%
% Syntax: q = tcn_std(  l1, l2, T, r0, n )
%  
%   l1  = Distance to the constant head boundary
%   l2  = Distance to the no flow head boundary
%   T  = Transmissivity
%   r0 = radius of the well
%   n  = optional argument to increase the number of image wells in the
%        computation of the solution. n should be infinite but for
%        practical purposes, n=10 is usually sufficient.
%
% Description:
%   The function computes both q and the drawdown map around the well by
%   applying the principle of superposition and adding a finite number of
%   imaginery wells.
%
%   The well is positionned in (0,0) on the maps.
%
% Exemple:
%   q=tcn_std(100,50,1e-2,0.1)
%
% See also: tcn_drw, tcn_dls



% Number of pairs of imaginery wells
if(nargin==4) % Default value if not given by user
  n=10;
end

% Domain size
l=l1+l2;
beta=l1/l;

% Build the vectors X,Y of the location of the imaginery wells
% and the vector S such that S=1 for a pumping
% well and S=-1 for an injection well
i=1:n;
y1=(i-1)*2*l+2*beta*l;
sy1=(-1).^i;
y2=-i*2*l+2*beta*l;
sy2=(-1).^(i+1);
y3=2*i*l;
y4=-y3;
yr=[0,y1,y2,y3,y4];
xr=0;
sr=[1,sy1,sy2,sy1,sy1];
k=1;
for i=1:size(xr,2)
    for j=1:size(yr,2)
        X(k)=xr(i);
        Y(k)=yr(j);
        S(k)=sr(j);
        k=k+1;
    end
end

% Map of the drawdown within the aquifer
nm=30;
dxm=(2*l)./nm;
dym=l./nm;
[xm,ym] = meshgrid(-l:dxm:l,-l2:dym:l1);
rm=sqrt(xm.^2+ym.^2);
sm=zeros(nm+1,nm+1);
s0=0;
s1=0;
for i=1:size(X,2)
    sm=sm+S(i)*log((X(i)-xm).^2+(Y(i)-ym).^2);
    s0=s0+S(i)*log((X(i)).^2+(Y(i)-r0).^2);
    s1=s1+S(i)*log((X(i)).^2+(Y(i)-l1).^2);
end
c=l1/(s1-s0);
q=4*pi*T*c;
hm=c.*(sm-s0);
hm(find(rm<=r0))=0;

% draw the map of the boundaries and locate the pumping well
figure(1)
clf
hold on
contourf(xm,ym,hm)
axis image
colorbar
title('Hydraulic heads')
xlabel('x')
ylabel('z')

% draw the real boundaries
% xl1=[-3*l,3*l]; yl1=[l1,l1];
% xl2=[-3*l,3*l]; yl2=[-l2,-l2];
% hold on
% plot(xl1,yl1,'b',xl2,yl2,'r')

% % Plot the pumping wells
% i=find(S>0);
% plot(X(i),Y(i),'+') 
% % Plot the injection wells
% i=find(S<0);
% plot(X(i),Y(i),'o')

% Surface plot of the heads
figure(2)
clf
surf(xm,ym,hm)
shading interp

% Plot of the norm of the velocity in logscale
% figure(3)
% clf
% [px,py] = gradient(hm,dxm,dym);
% v=T*sqrt(px.^2+py.^2);
% contourf(xm,ym,log(v))
% colorbar