C------------------------------------------------------------------------
SUBROUTINE Dyf4 ! Dla programu MultiFlower Wersja 1.0
C-------------------------------------------------------------------------
C.....Oblicza tylko cztery pierwsze składowe strumienia dyfuzyjnago
c na wybranym boku l, t, r, b - pierwsza składowa dyfuzyja jest równa 0
c składowa druga i trzecia to składowe napręzeń lepkich i turbulentnych,
c sładowa czwarta to dyfuzyjny strumień energii
c
c Przypadki(typ bloku l,r,t,b)
c
c 1) Regularna objętość skończona,
c 2) Bok jest scianką,
c 3) Bok jest wlotem / wylotem
c 4) Bok styka się z pasem (lentą)
C.....
c.....Przypadek 1. (regularna objętość skończona)
c.....1.1 Odszukanie wartosci Q w osmiu punktach
c wektor Q na granicach l,t,r,b komorki
c Q = [R,u,v,p,c_k] = [q1,q2,q3,q4,q5]
c c_k(ilosc skladnikow) - 1
c..... Oznaczenia: (R) gęstość; (u,v) prędkoęci w kierunku x,y;
c..... (p)ciśnienie; (c){w opracowaniu c:=c_k:=Y_k}liczba skladników
c.....1.2a Utworzenie Q_N i Q_S dla scianki w kierunku ksi
ij = (i-1) * ny + j + ijbaz ! ( i , j )
ij_e = i * ny + j + ijbaz ! ( i+1 , j )
ij_ne = (i-1) * ny + j + ijbaz ! ( i+1 , j+1 )
ij_n = (i-1) * ny + j + ijbaz + 1 ! ( i , j+1 )
ij_wn = (i-2) * ny + j + ijbaz + 1 ! ( i-1 , j+1 )
ij_w = (i-2) * ny + j + ijbaz ! ( i-1 , j )
ij_sw = (i-2) * ny + j + ijbaz - 1 ! ( i-1 , j-1 )
ij_s = (i-1) * ny + j + ijbaz - 1 ! ( i , j-1 )
ij_es = i * ny + j + ijbaz - 1 ! ( i+1 , j-1 )
c.....dla kierunku ksi
Ksi_N = 0.25 * ( ksi(ij_n) + ksi(ij_ne) + ksi(ij_e) + ksi(ij) )
Ksi_S = 0.25 * ( ksi(ij_e) + ksi(ij_es) + ksi(ij_s) + ksi(ij) )
Ksi_W = ksi(ij)
Ksi_E = ksi(ij_e)
Eta_N = 0.25 * ( eta(ij_n) + eta(ij_ne) + eta(ij_e) + eta(ij) )
Eta_S = 0.25 * ( eta(ij_e) + eta(ij_es) + eta(ij_s) + eta(ij) )
Eta_W = eta(ij)
Eta_E = eta(ij_e)
Ksi_NE = Ksi_E - Ksi_N
Ksi_SE = Ksi_E - Ksi_S
Ksi_NW = Ksi_N - Ksi_W
c.... Utworzenie przyrostów dla kierunku ksi
DeltaKsi_E = Ksi_NE - Ksi_SE
DeltaKsi_N = Ksi_NW - Ksi_NE
DeltaEta_E = Eta_NE - Eta_SE
DeltaEta_N = Eta_NW - Eta_NE
c.....Rekonstrukcja wartości U_N U_S
Q1_N = 0.25 * (q1(ii_ne) + q1(ii_e) + q1(ii_n) + q1(ii) )
Q1_S = 0.25 * (q1(ii_es) + q1(ii_e) + q1(ii_s) + q1(ii) )
Q1_E = q1(ij_e)
Q1_W = q1(ij)
Q2_N = 0.25 * (q2(ii_ne) + q2(ii_e) + q2(ii_n) + q2(ii) )
Q2_S = 0.25 * (q2(ii_es) + q2(ii_e) + q2(ii_s) + q2(ii) )
Q2_E = q2(ij_e)
Q2_W = q2(ij)
Q3_N = 0.25 * (q3(ii_ne) + q3(ii_e) + q3(ii_n) + q3(ii) )
Q3_S = 0.25 * (q3(ii_es) + q3(ii_e) + q3(ii_s) + q3(ii) )
Q3_E = q3(ij_e)
Q3_W = q3(ij)
Q4_N = 0.25 * (q4(ii_ne) + q4(ii_e) + q4(ii_n) + q4(ii) )
Q4_S = 0.25 * (q4(ii_es) + q4(ii_e) + q4(ii_s) + q4(ii) )
Q4_E = q4(ij_e)
Q4_W = q4(ij)
C..... 1.3 Utworzenie pochodnych (dQ/dKsi) i (dQ/dEta) w punkcie (1+1/2, j)
c Aksi - powierzchnia ścianki(dla 2d - długość boku) prostpadła do
c kierunku ksi
c Aksi1 = Aksi^(-1)
Aksi = 0.5 * ( eta(ij_ne) - eta(ij_es) )
Aksi1 = 1/Aksi
dq2dKsi = Aksi1 * ( Q2_E * DeltaEta_E + Q2_N * DeltaEta_N +
& Q2_W * DeltaEta_W + Q2_S * DeltaEta_S )
dq2dEta = -Aksi1 * ( Q2_E * DeltaKsi_E + Q2_N * DeltaKsi_N +
& Q2_W * DeltaKsi_W + Q2_S * DeltaKsi_S )
dq3dKsi = Aksi1 * ( Q3_E * DeltaEta_E + Q3_N * DeltaEta_N +
& Q3_W * DeltaEta_W + Q3_S * DeltaEta_S )
dq3dEta = -Aksi1 * ( Q3_E * DeltaKsi_E + Q3_N * DeltaKsi_N +
& Q3_W * DeltaKsi_W + Q3_S * DeltaKsi_S )
dq4dKsi = Aksi1 * ( Q4_E * DeltaEta_E + Q4_N * DeltaEta_N +
& Q4_W * DeltaEta_W + Q4_S * DeltaEta_S )
dq4dEta = -Aksi1 * ( Q4_E * DeltaKsi_E + Q4_N * DeltaKsi_N +
& Q4_W * DeltaKsi_W + Q4_S * DeltaKsi_S )
dq4dKsi = Aksi1 * ( Q4_E * DeltaEta_E + Q4_N * DeltaEta_N +
& Q4_W * DeltaEta_W + Q4_S * DeltaEta_S )
dq4dEta = -Aksi1 * ( Q4_E * DeltaKsi_E + Q4_N * DeltaKsi_N +
& Q4_W * DeltaKsi_W + Q4_S * DeltaKsi_S )
c..... koniec 1.2a
c.....1.2b Utworzenie Q_N i Q_S dla scianki w kierunku eta
c dla kierunku eta podobnie jak lda kierunku ksi...
c..... koniec 1.2a
c.....1.4 Określenie efektywnej turbulencji(Procedura Turbulencja)
c........ oraz współczynnika przewodności,
c mi_l lepkość laminarna według formuły Stutherlanda
c T0 tempetatura spoczynkowa na wylocie
mi_l = T**3/2 * (( 1.0 + 110.4 * T0 ) / ( T + 110.4 * T0 ))
mi_ef = mi_l + mi_t
C.....1.5 Utworzenie TauXX, TauXY, TauYY
TauXX = mi_ef * ( 4 * ( dKsidX * dUdKsi + dNidX * dUdEta ) - 2 *
$( dKsidY * dVdKsi + dEtadY * dVdEta ) )
TauXY = mi_ef * ( dKsiDY * dUdKsi + dEtadY * dUdEta + dKsidX *
& dVdKsi + dEtadX * dVdEta )
TauYY = mi_ef/3 * ( -2 * ( dKsidX * dUdKsi + dEtadX * dUdEta ) +
$ 4 * ( dKsidY * dVdKsi + dEtadY * dVdEta ) )
C.....1.6 Utworzenie składowych strumieni ciepła q_x, q_y
Pom=Mi_ef/Pr_ef*(1.0 / (Kappa-1.0) )
q_x = Pom * ( dKsidX * dTdKsi + dEtadX * dTdEta )
q_y = Pom * ( dKsiY * dTdKsi + dEtadY * dTdEta )
C.....1.7 Zbudowanie E^ni lub F^ni
c E_ev - czwarta składowe E_ni w układzie kartezjańskim
c F_ev - czwarta składowe F_ni w układzie kartezjańskim
C Zbudowanie E^ni =[E1_ni ,E2_ni ,E3_ni ,E4_ni]' w krzywoliniowym
c układzie
U = dKsidX * q2 + dKsidY * q3 ! sładowe kontrawariantne
V = dEtadX * q2 + dEtadY * q3 ! prędkości
E_ev = q2 * TauXX + q3 * TauXY + q_x
F_ev = q2 * TauXY + q3 * TauYY + q_y
E1_ni = 0.0
E2_ni = J1 * (dKsidX * TauXX + dKsidY * TauXY )
E3_ni = J1 * (dKsidX * TauXY + dKsidY * TauYY )
E4_ni = J1 * (dKsidX * E_ev + dKsidY * E_ev )
C Zbudowanie F^ni = [F1_ni ,F2_ni ,F3_ni ,F4_ni]' w krzywoliniowym
c układzie
F1_ni = 0.0
F2_ni = J1 * (dEtadX * TauXX + dEtadY * TauXY )
F3_ni = J1 * (dEtadX * TauXY + dEtadY * TauYY )
F4_ni = J1 * (dEtadX * F_ev + dEtadY * F_ev )
end ! SUBROUTINE Dyf4 dla programu MultiFlower