diff --git a/Source/ccib.f90 b/Source/ccib.f90 index 02067cc5704..a55eacb4427 100644 --- a/Source/ccib.f90 +++ b/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