Merge pull request #13331 from marcosvanella/master

FDS Source : Add CFACEs cartesian advection velocity for vorticity term in mome…

Source/ccib.f90

@@ -394,7 +394,12 @@ SUBROUTINE CUTFACE_VELOCITIES(NM,UU,VV,WW,CUTFACES)
 REAL(EB), POINTER, DIMENSION(:,:,:), INTENT(INOUT) :: UU,VV,WW
 LOGICAL, INTENT(IN) :: CUTFACES
 TYPE(CC_CUTFACE_TYPE), POINTER :: CF=>NULL()
-INTEGER :: ICF,I,J,K,X1AXIS
+TYPE(CFACE_TYPE), POINTER :: CFA=>NULL()
+TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
+TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
+INTEGER :: ICF,ICFA,JCF,I,J,K,X1AXIS
+REAL(EB):: AREA,VELN(3),PREDFCT
+
 
 CUTFACES_IF : IF (CUTFACES) THEN ! USE CUT_FACE(ICF)%VEL_CF
    DO ICF=1,MESHES(NM)%N_CUTFACE_MESH+MESHES(NM)%N_GCCUTFACE_MESH
@@ -410,6 +415,29 @@ SUBROUTINE CUTFACE_VELOCITIES(NM,UU,VV,WW,CUTFACES)
       END SELECT
    ENDDO
 
+   PREDFCT = 0._EB; IF(PREDICTOR) PREDFCT = 1._EB
+   ! CFACEs, set velocity in underlaying solid cartesian faces to be used in VELOCITY_FLUX:
+   DO ICF=1,MESHES(NM)%N_CUTFACE_MESH
+      CF => CUT_FACE(ICF); IF(CF%STATUS/=CC_INBOUNDARY) CYCLE
+      I = CF%IJK(IAXIS); J = CF%IJK(JAXIS); K = CF%IJK(KAXIS)
+      ! Area Average velocity for boundary CFACEs:
+      AREA = 0._EB; VELN(IAXIS:KAXIS) = 0._EB
+      DO JCF=1,CF%NFACE
+         ICFA=CF%CFACE_INDEX(JCF); IF(ICFA<1) CYCLE
+         CFA => CFACE(ICFA)
+         BC => BOUNDARY_COORD(CFA%BC_INDEX)
+         B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
+         AREA = AREA+CFA%AREA
+         VELN(IAXIS:KAXIS) = VELN(IAXIS:KAXIS) - &
+                             (PREDFCT*B1%U_NORMAL+(1._EB-PREDFCT)*B1%U_NORMAL_S)*CFA%AREA*BC%NVEC(IAXIS:KAXIS)
+      ENDDO
+      VELN(IAXIS:KAXIS) = VELN(IAXIS:KAXIS)/(AREA+TWO_EPSILON_EB)
+      ! Distribute into Solid cartesian faces:
+      WHERE(FCVAR(I-1:I,J,K,CC_FGSC,IAXIS)==CC_SOLID) UU(I-1:I,J,K) = VELN(IAXIS)
+      WHERE(FCVAR(I,J-1:J,K,CC_FGSC,JAXIS)==CC_SOLID) VV(I,J-1:J,K) = VELN(JAXIS)
+      WHERE(FCVAR(I,J,K-1:K,CC_FGSC,KAXIS)==CC_SOLID) WW(I,J,K-1:K) = VELN(KAXIS)
+   ENDDO
+
 ELSE CUTFACES_IF ! USE CUT_FACE(ICF)%VEL_CRT
    DO ICF=1,MESHES(NM)%N_CUTFACE_MESH+MESHES(NM)%N_GCCUTFACE_MESH
       CF => CUT_FACE(ICF); IF(CF%STATUS/=CC_GASPHASE) CYCLE
@@ -424,6 +452,14 @@ SUBROUTINE CUTFACE_VELOCITIES(NM,UU,VV,WW,CUTFACES)
       END SELECT
    ENDDO
 
+   DO ICF=1,MESHES(NM)%N_CUTFACE_MESH
+      CF => CUT_FACE(ICF); IF(CF%STATUS/=CC_INBOUNDARY) CYCLE
+      I = CF%IJK(IAXIS); J = CF%IJK(JAXIS); K = CF%IJK(KAXIS)
+      WHERE(FCVAR(I-1:I,J,K,CC_FGSC,IAXIS)==CC_SOLID) UU(I-1:I,J,K) = 0._EB
+      WHERE(FCVAR(I,J-1:J,K,CC_FGSC,JAXIS)==CC_SOLID) VV(I,J-1:J,K) = 0._EB
+      WHERE(FCVAR(I,J,K-1:K,CC_FGSC,KAXIS)==CC_SOLID) WW(I,J,K-1:K) = 0._EB
+   ENDDO
+
 ENDIF CUTFACES_IF
 RETURN
 END SUBROUTINE CUTFACE_VELOCITIES