Cálculos para o aterramento de uma subestação

projeto_aterramento2.m

clc
clearvars
close all

% resistividade aparente do solo [ohm*m]
rhoa = 448;

% resistividade da camada superior do solo [ohm*m]
rhoi = 472;

% resistividade do material de acabamento 
% da superfície da área da subestação [ohm*m]
rhos = 3000;

% corrente máxima de curto-circuito fase-terra [A]
Icft = 871;

% tempo de duração da corrente de curto-circuito fase-terra [s]
Tf = 0.5;

Vc = 15e3;
Ift_max = 55e3;
Lbrita = 0.15;

%% Tabelas

tabelaA30.Dh = [1/2 3/4 1];
tabelaA30.De = [2 3 4 5 9 12];
tabelaA30.A = [.2292 .1528 .1149 .0917 .0509 .0382; ...
               .2443 .1629 .1222 .0977 .0543 .0407; ...
               .2563 .1709 .1282 .1025 .0570 .0427];

tabelaA24.Dh = [1/2 3/4 1];
tabelaA24.De = [2 3 4 5 9 12];
tabelaA24.A = [.1898 .1266 .0949 .0759 .0422 .0316; ...
               .2028 .1352 .1014 .0811 .0450 .0338; ...
               .2132 .1421 .1066 .0853 .0474 .0355];

tabelaB.Nh = [4 9 16 25 36 49];
tabelaB.B = [2.7071 5.8917 8.5545 11.4371 14.0650 16.8933];


%% Número de condutores principais e de junção

Cm = 57; % comprimento
Lm = 41; % largura

% espaçamentos
Dl = Cm/17;
Dc = Lm/12;

% quantidades
Ncp = Cm/Dl + 1
Ncj = Lm/Dc + 1

%% Comprimento dos condutores da malha de terra [m]

Lcm = 1.05 * (Cm * Ncj + Lm * Ncp)

%% Coeficientes de ajuste Km

H = 0.5; % profundidade considerada da malha de terra
Dca = 14.33e-3; % diâmetro do condutor (magia da tabela 3.33, q nao existe)

% para condutores principais
Kmp = log(Dl^2/(4*pi*(Ncp-1)*H*Dca)) / (2*pi)

% para condutores de junção
Kmj = log(Dc^2/(4*pi*(Ncj-1)*H*Dca)) / (2*pi)

%% Coeficientes de ajuste Ks

% para condutores principais
Ksp = 1/pi * (1/(2*H) + 1/(Dl + H) + ...
    log(0.655 * (Ncp - 1) - 0.328) / Dl^2)

% para condutores de junção
Ksj = 1 / pi * (0.5 / H + 1 / (Dc + H) + ...
    log(0.655 * (Ncj - 1) - 0.328) / Dc^2)

%% Coeficientes de ajuste Ki

% para condutores principais
Kip = 0.65 + 0.172*Ncp

% para condutores de junção
Kij = 0.65 + 0.172*Ncj

%% Comprimento mínimo do condutor da malha [m]

% Adotar o maior produto Km*Ki para um direção considerada
if Kmp*Kip > Kmj*Kij
    Km = Kmp;
    Ki = Kip;
else
    Km = Kmj;
    Ki = Kij;
end

Lc = Km * Ki * rhoa * Icft * sqrt(Tf) / (0.116 + 0.174 * rhos)

% condição: Lcm >= Lc
if Lc > Lcm
    Lc, Lcm
    'Lc > Lcm'
    return
end

%% Tensão máxima de passo [V]

Epa = (116 + 0.7*rhos)/sqrt(Tf)

%% Tensão de passo existente na periferia da malha [V]

% Adotar o maior produto Ks*Ki para um direção considerada
if Ksp*Kip > Ksj*Kij
    Ks = Ksp;
    Ki = Kip;
else
    Ks = Ksj;
    Ki = Kij;
end

Eper = Ks * Ki * rhoi * Icft / Lcm

% condição: Epa >= Eper
if Eper > Epa
    Eper, Epa
    'Eper > Epa'
    return
end

%% Tensão máxima de toque [V]

Etm = (116 + 0.174*rhos)/sqrt(Tf)

%% Tensão de toque existente [V]

Ete = Km * Ki * rhoi * Icft / Lcm

% condição: Etm >= Ete
if Ete > Etm
    Ete, Etm
    'Ete > Etm'
    return
end

%% Corrente máxima de choque [mA]

Ich = 116/sqrt(Tf)

%% Corrente de choque existente devido à tensão de passo (sem brita) [mA]
% na periferia da malha

Ipmsb = 1000 * Eper / (1000 + 6 * rhoi)

% condição: Ipmsb < Ich
if Ipmsb > Ich
    disp('  Ipmsb     Ich')
    disp([Ipmsb, Ich])
    disp('Ipmsb > Ich: Deve-se usar brita')
end

%% Corrente de choque existente devido à tensão de passo (com brita) [mA]
% na periferia da malha

Ipmcb = 1000 * Eper / (1000 + 6 * (rhoi + rhos))

% condição: Ipmcb < Ich
if Ipmcb > Ich
    disp('  Ipmcb     Ich')
    disp([Ipmcb, Ich])
    disp('Ipmcb > Ich')
    return
end

%% Corrente de choque existente devido à tensão de toque existente (sem brita) [mA]

Itmsb = 1000 * Ete / (1000 + 1.5 * rhoi)

% condição: Itmsb < Ich
if Itmsb > Ich
    disp('  Itmsb     Ich')
    disp([Itmsb, Ich])
    disp('Itmsb > Ich: Deve-se usar brita')
end

%% Corrente de choque existente devido à tensão de toque existente (com brita) [mA]

Itmcb = 1000 * Ete / (1000 + 1.5 * (rhoi + rhos))

% condição: Itmcb < Ich
if Itmcb > Ich
    Itmcb, Ich
    disp('Itmcb > Ich')
    return
end

%% Corrente mínima de acionamento de relé de terra [A]

% Resistência considerada do corpo humano
Rch = 1000;

Ia = (Rch + 1.5 * rhos) * 9 * Lcm / (1000 * Km * Ki * rhoi)

%% Potenciais da região externa à malha [V]

% distância da cerca à periferia da malha
X = 5;
n = 2:Ncp-1;

% original
Kcx2 = log((H^2+X^2)*(H^2+(Dl+X)^2)/(H*Dca*(H^2+Dl^2)))/(2*pi) + ...
    1/pi*log(prod(1+X./(n*Dl)))

X1 = 4;
Kcx1 = log((H^2+X1^2)*(H^2+(Dl+X1)^2)/(H*Dca*(H^2+Dl^2)))/(2*pi) + ...
    1/pi*log(prod(1+X1./(n*Dl)))

% Kcx2 = log((H^2+X^2)*(H^2+(Dl+X)^2)/(H*Dca*(H+Dl)))/(2*pi) + ...
%     1/pi*log(prod(1+X./(n*Dl)))
% 
% X1 = 4;
% Kcx1 = log((H^2+X1^2)*(H^2+(Dl+X1)^2)/(H*Dca*(H+Dl)))/(2*pi) + ...
%     1/pi*log(prod(1+X1./(n*Dl)))

% original
% dEc = (Kcx2 - Kcx1) * rhoi * Icft / Lc
dEc = (Kcx2 - Kcx1) * rhoi * Icft / Lcm

% condição: dEc < Ete
if dEc > Ete
    dEc, Ete
    disp('dEc > Ete: A cerca necessita de aterramento.')
end

%% Resistência da malha de terra [ohm]

Amt = Cm*Lm;
R = sqrt(Amt/pi)

Rmc = rhoa/(4*R) + rhoa/Lcm

if (Vc <= 38e3) || (Vc >= 15e3)
    % condição: Rmc <= 10
    if Rmc >= 10
        disp('Rmc > 10')
    end
elseif Vc >= 69e3
    % condição: Rmc <= 5
    if Rmc >= 5
        disp('Rmc > 5')
    end
end

%% Resistência de aterramento de um eletrodo vertical [ohm]

Dh = 3/4;
Lh = 3;

Rel = rhoa/(2*pi*Lh)*log(400*Lh/(2.54*Dh))

%% Coeficiente de redução da resistência de um eletrodo vertical

Nh = 12;

if Lh == 3.0
    tabelaA = tabelaA30;
elseif Lh == 2.4
    tabelaA = tabelaA24;
end

% índice correspondente ao diametro do eletrodo
[~, iA] = min(abs(Dh-tabelaA.Dh));

% distância máxiima entre os eltrodos verticais
Dev = max(Dl,Dc)*4;

if Dev > max(tabelaA.De)
    jA = length(tabelaA.De);
    [~, jA] = min(abs(Dev-tabelaA.De));
    A = tabelaA.A(iA,jA)
else
    A = interp1(tabelaA.De,tabelaA.A(iA,:),Dev)
end

B = interp1(tabelaB.Nh,tabelaB.B,Nh)

Kh = (1 + A * B) / Nh

%% Resistência de aterramento do conjunto de eletrodos verticais

Rne = Kh * Rel

%% Resistência mútua dos cabos e eletrodos verticais

K = Cm/Lm
Lth = Nh * Lh

K1 = 1.14125 - 0.0425 * K
K2 = 5.49 - 0.1443 * K
S = Cm*Lm;

Rmu = rhoa / (pi * Lcm) * (log(2 * Lcm / Lth) + K1 * Lcm / sqrt(S) - K2 + 1)

%% Resistência total da malha

Rtm = (Rmc * Rne - Rmu^2) / (Rmc + Rne - 2 * Rmu)

red = (Rmc-Rtm)/Rmc*100;
disp(['Redução da resistência da malha de terra foi de ', num2str(red), '%'])