Local variables to be added
The following local variables need to be defined for the examples in this section:
integer iel, ielpdc, ikpdc
integer ilelt, nlelt
integer izone
double precision alpha, cosalp, sinalp, vit, ck1, ck2
double precision, dimension(:,:), pointer :: cvara_vel
integer, allocatable, dimension(:) :: lstelt
Initialization and finalization
The following initialization block needs to be added for the following examples:
At the end of the subroutine, it is recommended to deallocate the work array:
In theory Fortran 95 deallocates locally-allocated arrays automatically, but deallocating arrays in a symetric manner to their allocation is good pratice, and avoids using a different logic for C and Fortran.
Map field array
Body
Calculation and identification of the number of cells with imposed head loss term
if (iappel.eq.1.or.iappel.eq.2) then
2 calls:
iappel
= 1: Calculation of the number of cells where a head loss term is imposed: ncepdp
. Called once at the beginning of the calculation.
iappel
= 2: Identification of the cells where a head loss term is imposed: array icepdc(ncepdc)
. Called once at the beginning of the calculation.
- Note
- Do not use ckupdc in this section (it is defined with
iappel = 3)
- Use icepdc in this section only with
(iappel = 2)
To be completed by the user: cell selection
Head loss examples
Example 1: No head loss (default)
Example 2: Head losses defined by coordinates for zone
(4 <= x; 2 <= y <= 8)
No head losses else.
if (1.eq.0) then
izone = 0
ielpdc = 0
call getcel(
'X <= 6.0 and X >= 4.0 and Y >= 2.0 and'// &
'Y <= 8.0',nlelt,lstelt)
Generic subsection
For iappel = 1
, define ncepdp
, the number of cells with head losses. This is valid for both examples above.
izone = izone + 1
do ilelt = 1, nlelt
iel = lstelt(ilelt)
izcpdc(iel) = izone
ielpdc = ielpdc + 1
if (iappel.eq.2) icepdc(ielpdc) = iel
enddo
endif
Defining the number of cells with head losses
if (iappel.eq.1) then
ncepdp = ielpdc
endif
Computing the head loss coefficient values
Third call, at each time step
iappel = 3:
ckupdc:
compute head loss coefficients in the calculation coordinates, organized in order k11, k22, k33, k12, k13, k23
Note:
- make sure diagonal coefficients are positive. The calculation may crash if this is not the case, and no further check will be done
Example 1: head losses in direction x
Diagonal tensor : Example of head losses in direction x
if (.true.) return
do ielpdc = 1, ncepdp
iel=icepdc(ielpdc)
vit = sqrt(cvara_vel(1,iel)**2 + cvara_vel(2,iel)**2 + cvara_vel(3,iel)**2)
ckupdc(ielpdc,1) = 10.d0*vit
ckupdc(ielpdc,2) = 0.d0*vit
ckupdc(ielpdc,3) = 0.d0*vit
enddo
Example 2: alpha = 45 degres
3x3 tensor: Example of head losses at alpha = 45 degres x,y direction x
resists by ck1
and y
by ck2
ck2 = 0
represents vanes as follows: in coordinate system x
, y
Orthogonal reference frame sketch
if (.true.) return
ck1 = 10.d0
ck2 = 0.d0
Filling the cells of the head loss matrix
do ielpdc = 1, ncepdp
iel = icepdc(ielpdc)
vit = sqrt(cvara_vel(1,iel)**2 + cvara_vel(2,iel)**2 + cvara_vel(3,iel)**2)
ckupdc(ielpdc,1) = (cosalp**2*ck1 + sinalp**2*ck2)*vit
ckupdc(ielpdc,2) = (sinalp**2*ck1 + cosalp**2*ck2)*vit
ckupdc(ielpdc,3) = 0.d0
ckupdc(ielpdc,4) = cosalp*sinalp*(-ck1+ck2)*vit
ckupdc(ielpdc,5) = 0.d0
ckupdc(ielpdc,6) = 0.d0
enddo
endif