/////////////////////////////////Variavies obtidas nos sensores:////////////////////////

//Ta = temperatura ambiente (°C) - Sensor DHT 22
//UR = umidade relativa do ar (%) - Sensor DHT 22
//Patm = pressão atmosferica (Pa) - Sensor BMP 280
//Te = temperatura de entrado na placa 1 (°C) - SENSOR DSB1820
//Ts1 = temperatura de saída da placa 1 e entrada na placa 2 (°C) - SENSOR DSB1820
//Ts2 = temperatura de saída da placa 2 (°C) - SENSOR DSB1820
//Tw1 = temperatura superficie aquecedor 1 (°C) - SENSOR DSB1820
//Tw2 = temperatura superficie aquecedor 2 (°C) - SENSOR DSB1820
//Vel = velocidade do vento (m/s) - Anemometro Caseiro
//Ti1 = temperatura inferior aquecedor (°C) - SENSOR DSB1820
//Ti2 = temperatura inferior aquecedor 2 (°C) - SENSOR DSB1820
//vazao = vazão de agúa no aquecedor  (l/s) - Sensor Hall Medidor Fluxo Água 3/4 2-45l  -  Freq = 8.1 Q-3
//Rad = Radiação Solar (W/m²) - LDR 5mm -

////////////////////////Capacidades físicas da Ar:////////////////////////////////////

Kar = 0.03;                     //Condutiblidade Termica do ar (W/m°C)
Prar = 0.707;                   //N° de Prandl do ar (-)
Uar = 17.4 * pow(10,-6);        //Viscosidade Dinamica do ar (Pa.s)
Cpar = 1004.831565;             //Capacidade Calorifica do ar (J/kg°C)

///Caracteristicas Físicas do Aquecedor Solar

De = 1.71;                                        //Diametro externo aquecedor (m)
Di = 0.46;                                        //Diametro interno aquecedor (m)
Pi = 3.14159265359;                               //Número pi
Per = (De * Pi);                                  //Perimetro aquecedor (m)
Aes = ((De * De - Di *Di) * Pi * 0.25);           //Área de exposição solar (m²)
Lc = Aes / Per;                                   //Comprimento caracteristico (m)
di0 = 0.0195;                                     //Diametro interno do tubo (m)
de0 = 0.0225;                                     //Diametro externo do tubo (m)
L = 100;                                          //Comprimento Tubo (m)
At = de0*Pi*L                                     //Aréa total de Troca de calor (m²)
/////////////////////Massa Especifica do Ar Úmido  (kg/m³)////////////////////////////////////

Psat = (6.1078 * pow(10,((7.5 * (Ta)) / (Ta + 237.3)));                            //Pressão de saturação do ar
Pv = Psat * UR;                                                                   //Pressão de vapor do ar
Pd = Patm - Pv;                                                                   //Pressão dinamica
par = (Pd / (287.058 * (Ta + 273.15))) + (Pv / (461.495 * (273.15 + Ta)));        //Massa Especifica do Ar

///////////////Calculo do Coeficiente Convectivo placa 1 e 2:////////////////////////////////

//N° de Reynolds para Escoamento do Vento na Superficie do Aquecedor:

Re = (vel * par * Lc) / Uar;                                                                //N° de Reynolds
Gr1 = (9.81 * (1 / (Ta + 273.15) * (Tw1 - Ta) * pow(Lc,2)) / (Uar / par));                  //N° de Grashft  placa 1
Gr2 = (9.81 * (1 / (Ta + 273.15) * (Tw2 - Ta) * pow(Lc,2)) / (Uar / par));                  //N° de Grashft placa 2
Ri1 = Gr1 / pow(Re,2);                                                                      //N° de Richardson placa 1
Ri2 = Gr2 / pow(Re,2);                                                                      //N° de Richardson placa 2
Ra1 = Gr1 * Prar;                                                                           //N° de Rayleigh placa 1
Ra2 = Gr2 * Prar;                                                                           //N° de Rayleigh placa 2
Nuf1 = 1.72 * pow(Re,0.5) * pow(Prar,(1/3));                                                //Nussfelt forçado placa 1
Nuf2 = 1.72 * pow(Re,0.5) * pow(Prar,(1/3));                                                //Nussfelt forçado placa 2

///////Condicional para escolha do tipo de escoamento exterior:///////////////////////

////Nussfelt natural para placa 1 :

if (Ra1 >=104,Ra1 <= 109)
{
Nun1 = 0.59 * pow(Ra1,(1/4));
}
else if (Ra1 > 109,Ra1 <= 1012)
{
Nun1 = 0.13 * pow(Ra1,(1/3));
}
else
{
Nun1 = 0;
}
end if

////Nussfelt natural para placa 2 :

if (Ra2 >= 104,Ra2 <= 109)
{
Nun2 = 0.59 * pow(Ra2,(1/4));
}
else if (Ra2 > 109,Ra2 <= 1012)
{
Nun2 = 0.13 * pow(Ra2,(1/3));
}
else
{
Nun2 = 0;
}
end if

////Nussfelt forçada para placa 1

if (Ri1 <= 0.1)
{
Nu1 = Nuf1;
}
else if (Ri1 > 0.1,Ri1 <= 10)
{
Nu1 = Nun1 + Nuf1;
}
else
{
Nu1 = Nun1;
}
end if

////Nussfelt forçada para placa 2

if (Ri2 <= 0.1)
{
Nu2 = Nuf2;
}
else if (Ri2 > 0.1,Ri2 <= 10)
{
Nu2 = Nun2 + Nuf2;
}
else
{
Nu2 = Nun2;
}
end if

////Coeficiente Convectivo com o ar://////////////////////////////////////////////

hconv1 = Nu1 * Kar / Lc;           //Coeficiente convectivo externo placa 1 (W/m²°C)
hconv2 = Nu2 * Kar / Lc;           //Coeficiente convectivo externo placa 2 (W/m²°C)

////Quantidade de calor perdida por convecção com o ar:///////////////////////////

Qconv1 = hconv1 * Aes * (Tw1 - Ta);         //Calor perdido na placa 1 p/ ar (W)
Qconv2 = hconv1 * Aes * (Tw2 - Ta);        //Calor perdido na placa 2 p/ ar (W)

/////Fluxo de Calor por Radiação solar://////////////////////////////////////////

Qsolar = Rad * Aes * 0.9;                   //Calor que entra no sistema (W)

/////////////Perda de calor por irradiação://///////////////////////////////////
TaK = Ta + 273.15;
Tw1K = Tw1 + 273.15;
Tw21K = Tw2 + 273.15;
Qrad1 = 5.6697 * 10 ^ (-8) * Aes * 0.95 * ((pow(Tw1K,4)) - (pow((0.0552 * pow(TaK,1.5)),4)); //(W)
Qrad2 = 5.6697 * 10 ^ (-8) * Aes * 0.95 * ((pow(Tw2K,4)) - (pow((0.0552 * pow(TaK,1.5)),4)); //(W)

//////////////////////////Velocidade agua no interior do tubo://///////////////////////////////////////

Vagua = (4 * (vazao * 0.001) / (pow(0.0195,2) * Pi));//Velocidade de escoamento água nos tubos //(m/s)

//////////////////Temeperaturas acessorias para os calculos:///////////////////////

Twi1 = Tw1 - (((Qsolar * 0.35 * Log(0.0225 / 0.0195)) / (2 * Pi * 100 * 0.17)));     //Temperatura interna da parede tubo 1
Twi2 = Tw2 - (((Qsolar * 0.35 * Log(0.0225 / 0.0195)) / (2 * Pi * 100 * 0.17)));     //Temperatura interna da parede tubo 2
Tm1 = (Te1 + Ts1) * 0.5;                                                             //Temperatura média da água placa 1
Tm2 = (Ts1 + Ts2) * 0.5;                                                             //Temperatura média da água placa 2

////Curvas para cálculo dos elementos das caracteristicas físico-quimicas da água:////////////////////////////

//Massa Especifica da água (kg/m³) :

Rou1 = ((0.0035) * Tm1*Tm1) - 0.0722 * Tm1 + 1000.7;
Rou2 = ((0.0035) * Tm2*Tm2) - 0.0722 * Tm2 + 1000.7;

//Viscosidade dianamica da água ud1 (Pa.s):

ud1 = ((-2 * pow(10,- 9)) * pow(Tm1,3)) + ((6 * pow(10,- 7)) * pow(Tm1,2)) - (5 * pow(10,- 5) * Tm1) + 0.0017;
ud2 = ((-2 * pow(10,- 9)) * pow(Tm2,3)) + ((6 * pow(10,- 7)) * pow(Tm2,2)) - (5 * pow(10,- 5) * Tm2) + 0.0017;

//Condutibilidade termica da água(W/m°C) :

Ka1 = (1 * pow(10,(-8)) * pow((Tm1,3))) - 9 * pow(10,(-6)) * pow((Tm1,2)) + 0.002 * (Tm1) + 0.5604;
Ka2 = (1 * pow(10,(-8)) * pow((Tm2,3))) - 9 * pow(10,(-6)) * pow((Tm2,2)) + 0.002 * (Tm2) + 0.5604;

//Capacidade termica da água(J/kg°C):
Cp1 = ((-9 * pow(10,-5)) * pow(Tm1,3)) + (0.0291 * pow(Tm1,2)) - (1.9549 * Tm1 + 4210.2);
Cp2 = ((-9 * pow(10,-5)) * pow(Tm2,3)) + (0.0291 * pow(Tm2,2)) - (1.9549 * Tm2 + 4210.2);

//Numero de Prandl da água:

Prh1 = Cp1 * ud1 / Ka1;
Prh2 = Cp2 * ud2 / Ka2;

//Viscosidade dianamica da água a temperatura da parede do tubo:

ut1 = ((-2 * pow(10,- 9)) * pow(Twi1,3) + ((6 * pow(10,-7)) * pow(Twi1,2)) - (5 * pow(10,- 5 * Twi1)) + 0.0017;
ut2 = ((-2 * pow(10,- 9)) * pow(Twi2,3) + ((6 * pow(10,-7)) * pow(Twi2,2)) - (5 * pow(10,- 5 * Twi2)) + 0.0017;

//Número de Reynolds no interior do tubo:

rey1 = (Vagua * 0.0195 * Rou1) / (ud1);  //Tubo 1
rey2 = (Vagua * 0.0195 * Rou2) / (ud2);  //Tubo 2

//Definição do coeficiente de atrito por Colebrook placa 1 - equação implicita

ra = 0.001;        //Rugosidade aparente
Ed = ra / di0;     //Rugosidade Relativa


F01 = 0.3164 * pow(rey1,-0.25);
float erro1 = 1;
while (erro1>= 0.0001)
{
func1 = 1 / pow(F01,0.5) + 2 * Log(Ed / 3.72 + 2.51 / (rey1 * pow(F01,0.5))) / Log(10);

dfunc1 = -2.51 / (Log(10) * pow(F01,(3 / 2)) * rey1 * (Ed / 3.72 + 2.51 / (rey1 * pow(F01,0.5)))) - 1 / (2 * pow(F01,(3 / 2));

f1 = F01 - func1 / dfunc1;

erro1 = Abs(F01 - f1);

F01 = f1;

}
end while

// Definição do coeficiente de atrito por Colebrook placa 2 - equação implicita

F02 =0.3164 * pow(rey2,-0.25);
float erro1 = 1;
while (erro1>= 0.0001)
{
func2 = 1 / pow(F02,0.5) + 2 * Log(Ed / 3.72 + 2.51 / (rey1 * pow(F02,0.5))) / Log(10);

dfunc2 =  -2.51 / (Log(10) * pow(F02,(3 / 2)) * rey2 * (Ed / 3.72 + 2.51 / (rey2 * pow(F02,0.5)))) - 1 / (2 * pow(F02,(3 / 2));

f2 = F02 - func2 / dfunc2;

erro2 = Abs(F02 - f2);

F02 = f2;
}
end while

///////////////////////////Calculo do Nussfelt no interior do tubo 1:///////////////////////////////////

a1 = (f1 / 8);
b1 = (rey1 - 1000) * Prh1;
c1 = 12.7 * ((pow((f1 / 8),0.5)) * ((pow(Prh1,(2 / 3))) - 1));
d1 = (ud1 / ut1);
a2 = (f2 / 8);
b2 = (rey2 - 1000) * Prh2;
c2 = 12.7 * ((pow((f2 / 8),0.5)) * ((pow(Prh2,(2 / 3))) - 1));
d2 = (ud2 / ut2);

if (rey1 <= 2200)
{
Nus1 = 3.657;
}
else  (rey1 > 2200)
{
Nus1 = ((a1 * b1) / (c1)) * d1;
}
end if

if  (rey2 <= 2200)
{
Nus2 = 3.657;
}
else if (rey2 > 2200)
{
Nus2 = ((a2 * b2) / (c2)) * d2;
}
end if

//Nusfelt corrigido para serpentina:

Nus1c = 1# + 3.4 * (2 * 0.0195 / 1.71) * Nus1;
Nus2c = 1# + 3.4 * (2 * 0.0195 / 1.71) * Nus2;

//Coeficiente de troca de calor no interior dos tubos:

h1 = Nusc1 * Ka1 / 0.0195;
h2 = Nusc2 * Ka2 / 0.0195;

//Fluxo de calor no interior dos tubos:

Q1 = h1 * Aes * (Twi1 - Tm1);
Q2 = h2 * Aes * (Twi2 - Tm2);

//Fluxo de calor por condução:

lt = 0.006;    //expessura telha
kT = 0.35;     //coef de cond. telha

Q1cond = kT * Aes * (Ti1 - ta) / lt;
Q2cond = kT * Aes * (Ti2 - ta) / lt;

///Coeficiente de troca de calor por condução:

hcond = kT / lt;

///Calculo do calor util pela vazao massica:

Vmass1 = vazao * Rou1;    //vazao massica 1
Vmass2 = vazao * Rou2;    //vazao massica 2


///Calculo do Coeficiente Global metodo 1:

U1 = 1 / ((1 / hconv1) + (hcond) + ((Math.Log(0.0225 / 0.0195)) / (2 * Pi * 100 * 0.17)) + (1 / h1));
U2 = 1 / ((1 / hconv2) + (hcond) + ((Math.Log(0.0225 / 0.0195)) / (2 * Pi * 100 * 0.17)) + (1 / h2));


// Calculo do Rendimento total:////////////////

//// Calculo da Temperatura de saída alternando vazão://///////////////////////////////////////
float Vv = 0.5;
while(Vv ==45); // varia a vazão de 0.5l/s a 45l/s
{
//Iteração para temperatura, Tscalc = Temperatura de saída calculada, varia até que a última Tscalc calculada diminuida da anterior tenha uma diferença <=0.01.

 float Tscalc1 = Ts1;                              //define a temperatura de saida calculada inicial como a ultima medida pelo sensor
 int N1 = 0;                                      //define o indice da Tscal
 diferenca1 = abs(Tscalc1(N1) - Tscalc1(N1-1))    //calcula a diferença da temperatura calculada atual e a anterior
for (diferenca1 <= 0.01; N1++)                    //estipula a diferença minima para o calculo parar
{
Tmcalc1 = (Tscalc1 + Te1)/2;                                                                                    //calcula a temperatura media para o calculo da temperatura de saída
Cpc1 =((-9 * pow(10,-5)) * pow(Tmcalc1,3)) + (0.0291 * pow(Tmcalc1,2)) - (1.9549 * Tmcalc1 + 4210.2);           //calcula a capacidade especifica para o calculo da temp de saida
Rouc1 = ((0.0035) * pow(Tmcalc1,2) - 0.0722 * Tmcalc1 + 1000.7;                                                  //calcula a massa especifica calculo da temp de saida
Vmassc1 = (Vv * Rouc1)                                                                                              //calcula a vazao massica p calculo da temp de saida
TScalc1(N1) = (((Vmassc1 * Cpc1 * Te1) - (0.5 * U * At * ta) + U * At * Ta)) / (Vmassc1 * Cpc1 - U * At * 0.5)  //calcula a temp de saida
}
Tsc1 = Tscalc1;      //grava a temp de saida calculada

//Iteração para temperatura na placa 2:

float Tscalc2 = Ts2;
int N2 = 0;
diferenca2 = abs(Tscalc1(N2) - Tscalc1(N2-1));
for (diferenca2 <= 0.01; N2++)
{
Tmcalc2 = (Tscalc2 + Tsc1)/2;                                                                                    //calcula a temperatura media para o calculo da temperatura de saída
Cpc2 =((-9 * pow(10,-5)) * pow(Tmcalc2,3)) + (0.0291 * pow(Tmcalc2,2)) - (1.9549 * Tmcalc2 + 4210.2);           //calcula a capacidade especifica para o calculo da temp de saida
Rouc2 = ((0.0035) * pow(Tmcalc2,2) - 0.0722 * Tmcalc2 + 1000.7;                                                  //calcula a massa especifica calculo da temp de saida
Vmassc2 = (Vv * Rouc2)                                                                                              //calcula a vazao massica p calculo da temp de saida
TScalc2(N2) = (((Vmassc2 * Cpc2 * Tsc1) - (0.5 * U * At * ta) + U * At * Ta)) / (Vmassc2 * Cpc2 - U * At * 0.5)  //calcula a temp de saida
}
Tsc2 = Tscalc2;      //grava a temp de saida calculada
}

Qutil1 = Vmassc1 * Cpc1 * (Tsc1 - Te1); //Calor Útil 1
Qutil2 = Vmassc2 * Cpc2 * (Tsc2 - Tsc1); //Calor Útil 2

Rend = (Quti11 + Qutil2) / (Qrad * 2);

**Equação para escolher a Vv que retorna o maior Rendimento....