Pmma

Source/func.f90

@@ -1731,6 +1731,9 @@ SUBROUTINE PACK_BOUNDARY_COORD(NM,IC,RC,OS,BC_INDEX,UNPACK_IT,COUNT_ONLY)
 IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%IIG,UNPACK_IT)
 IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%JJG,UNPACK_IT)
 IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%KKG,UNPACK_IT)
+IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%IIS,UNPACK_IT)
+IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%JJS,UNPACK_IT)
+IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%KKS,UNPACK_IT)
 IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%IOR,UNPACK_IT)
 
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),BC%X,UNPACK_IT)
@@ -1957,7 +1960,7 @@ END SUBROUTINE INITIALIZE_BOUNDARY_ONE_D
 
 SUBROUTINE INITIALIZE_BOUNDARY_PROP1(NM,B1_INDEX,SURF_INDEX)
 
-USE GLOBAL_CONSTANTS, ONLY: SIGMA,TMPA,RHOA,TMPA4,N_TRACKED_SPECIES, RADIATION
+USE GLOBAL_CONSTANTS, ONLY: SIGMA,TMPA,RHOA,TMPA4,N_TRACKED_SPECIES,RADIATION,MU_AIR_0,RSUM0
 INTEGER, INTENT(IN) :: NM,B1_INDEX,SURF_INDEX
 TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
 TYPE(SURFACE_TYPE), POINTER :: SF
@@ -1977,8 +1980,13 @@ SUBROUTINE INITIALIZE_BOUNDARY_PROP1(NM,B1_INDEX,SURF_INDEX)
 B1%TMP_F           = SF%TMP_FRONT
 B1%TMP_G           = TMPA
 B1%TMP_B           = SF%TMP_BACK
+B1%TMP_FS          = TMPA
+B1%RHO_FS          = RHOA
+B1%U_TANG_FS       = 0._EB
 B1%RHO_F           = RHOA
 B1%RHO_G           = RHOA
+B1%RSUM_G          = RSUM0
+B1%MU_G            = MU_AIR_0
 B1%TMP_F_OLD       = SF%TMP_FRONT
 B1%BURN_DURATION   = SF%BURN_DURATION
 
@@ -2045,13 +2053,18 @@ SUBROUTINE PACK_BOUNDARY_PROP1(NM,IC,RC,LC,OS,B1_INDEX,UNPACK_IT,COUNT_ONLY,SURF
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%TMP_F,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%TMP_G,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%TMP_B,UNPACK_IT)
+RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%TMP_FS,UNPACK_IT)
+RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%RHO_FS,UNPACK_IT)
+RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%U_TANG_FS,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%U_NORMAL,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%U_NORMAL_S,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%U_NORMAL_0,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%U_TANG,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%U_IMPACT,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%RHO_F,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%RHO_G,UNPACK_IT)
+RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%RSUM_G,UNPACK_IT)
+RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%MU_G,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%RDN,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%K_G,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%Q_DOT_G_PP,UNPACK_IT)
@@ -2723,6 +2736,81 @@ SUBROUTINE ASSIGN_PRESSURE_ZONE(NM,II,JJ,KK,I_ZONE,I_ZONE_OVERLAP)
 END SUBROUTINE ASSIGN_PRESSURE_ZONE
 
 
+!> \brief Determine distance from a boundary cell to the nearest solid object or mesh boundary
+!> \param NM Mesh number
+!> \param IOR Index of the boundary cell orientation
+!> \param IIG i-coordinate of the point (m)
+!> \param JJG j-coordinate of the point (m)
+!> \param KKG k-coordinate of the point (m)
+!> \param CLEARANCE Distance from boundary to nearest solid
+
+SUBROUTINE COMPUTE_BOUNDARY_CLEARANCE(NM,IOR,IIG,JJG,KKG,CLEARANCE)
+
+INTEGER, INTENT(IN) :: NM,IOR,IIG,JJG,KKG
+REAL(EB), INTENT(OUT) :: CLEARANCE
+INTEGER :: N,IC,II,JJ,KK
+TYPE (MESH_TYPE), POINTER :: M
+TYPE (OBSTRUCTION_TYPE), POINTER :: OB
+
+M => MESHES(NM)
+
+CLEARANCE = 0._EB
+
+II = IIG
+JJ = JJG
+KK = KKG
+
+SEARCH_LOOP: DO
+
+   IC = M%CELL_INDEX(II,JJ,KK)
+   IF (M%CELL(IC)%SOLID) EXIT SEARCH_LOOP
+
+   SELECT CASE(IOR)
+      CASE(-1)
+         CLEARANCE = CLEARANCE + M%DX(II)
+         IF (II==1) EXIT SEARCH_LOOP
+         II = II-1
+      CASE( 1)
+         CLEARANCE = CLEARANCE + M%DX(II)
+         IF (II==M%IBAR) EXIT SEARCH_LOOP
+         II = II+1
+      CASE(-2)
+         CLEARANCE = CLEARANCE + M%DY(JJ)
+         IF (JJ==1) EXIT SEARCH_LOOP
+         JJ = JJ-1
+      CASE( 2)
+         CLEARANCE = CLEARANCE + M%DY(JJ)
+         IF (JJ==M%JBAR) EXIT SEARCH_LOOP
+         JJ = JJ+1
+      CASE(-3)
+         CLEARANCE = CLEARANCE + M%DZ(KK)
+         IF (KK==1) EXIT SEARCH_LOOP
+         KK = KK-1
+      CASE( 3)
+         CLEARANCE = CLEARANCE + M%DZ(KK)
+         IF (KK==M%KBAR) EXIT SEARCH_LOOP
+         KK = KK+1
+   END SELECT
+
+   ! Look for thin obstructions bordering the current cell
+
+   DO N=1,M%N_OBST
+      OB => M%OBSTRUCTION(N)
+      SELECT CASE(IOR)
+         CASE(-1) ; IF (II  ==OB%I1.AND.II  ==OB%I2.AND.JJ>OB%J1.AND.JJ<=OB%J2.AND.KK>OB%K1.AND.KK<=OB%K2) EXIT SEARCH_LOOP
+         CASE( 1) ; IF (II-1==OB%I1.AND.II-1==OB%I2.AND.JJ>OB%J1.AND.JJ<=OB%J2.AND.KK>OB%K1.AND.KK<=OB%K2) EXIT SEARCH_LOOP
+         CASE(-2) ; IF (JJ  ==OB%J1.AND.JJ  ==OB%J2.AND.II>OB%I1.AND.II<=OB%I2.AND.KK>OB%K1.AND.KK<=OB%K2) EXIT SEARCH_LOOP
+         CASE( 2) ; IF (JJ-1==OB%J1.AND.JJ-1==OB%J2.AND.II>OB%I1.AND.II<=OB%I2.AND.KK>OB%K1.AND.KK<=OB%K2) EXIT SEARCH_LOOP
+         CASE(-3) ; IF (KK  ==OB%K1.AND.KK  ==OB%K2.AND.II>OB%I1.AND.II<=OB%I2.AND.JJ>OB%J1.AND.JJ<=OB%J2) EXIT SEARCH_LOOP
+         CASE( 3) ; IF (KK-1==OB%K1.AND.KK-1==OB%K2.AND.II>OB%I1.AND.II<=OB%I2.AND.JJ>OB%J1.AND.JJ<=OB%J2) EXIT SEARCH_LOOP
+      END SELECT
+   ENDDO
+
+ENDDO SEARCH_LOOP
+
+END SUBROUTINE COMPUTE_BOUNDARY_CLEARANCE
+
+
 !> \brief Determine the number of 1-D cells in a layer of solid material
 !> \param DIFFUSIVITY \f$ k/(\rho c) \; (\hbox{m}^2/\hbox{s}) \f$, used to determine cell size
 !> \param LAYER_THICKNESS Thickness of the material layer (m)

Source/init.f90

@@ -2878,7 +2878,7 @@ SUBROUTINE INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,TT)
 
 USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
 USE MEMORY_FUNCTIONS, ONLY: ALLOCATE_STORAGE
-USE GEOMETRY_FUNCTIONS, ONLY: SEARCH_OTHER_MESHES
+USE GEOMETRY_FUNCTIONS, ONLY: SEARCH_OTHER_MESHES, COMPUTE_BOUNDARY_CLEARANCE
 USE COMP_FUNCTIONS, ONLY: SHUTDOWN
 USE PHYSICAL_FUNCTIONS, ONLY: GET_SPECIFIC_GAS_CONSTANT
 USE CONTROL_VARIABLES, ONLY : CONTROL
@@ -2888,7 +2888,7 @@ SUBROUTINE INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,TT)
 INTEGER, INTENT(OUT) :: IERR
 REAL(EB), INTENT(IN) :: TT
 REAL(EB) :: PX,PY,PZ,T_ACTIVATE,XIN,YIN,ZIN,DIST,XW,YW,ZW,RDN,AW,TSI,&
-            ZZ_GET(1:N_TRACKED_SPECIES),RSUM_F,R1,RR,DELTA
+            ZZ_GET(1:N_TRACKED_SPECIES),RSUM_F,R1,RR,DELTA,XI,YJ,ZK,CLEARANCE
 INTEGER  :: N,SURF_INDEX_NEW,IIG,JJG,KKG,IIO,JJO,KKO,IC,ICG,ICO,NOM_CHECK(0:1),BOUNDARY_TYPE,FI,VENT_INDEX_FOUND
 LOGICAL :: ALIGNED
 TYPE (MESH_TYPE), POINTER :: M,MM
@@ -3261,6 +3261,25 @@ SUBROUTINE INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,TT)
    CASE( 3) ; BC%X1=M%X(I-1) ; BC%X2=M%X(I)   ; BC%Y1=M%Y(J-1) ; BC%Y2=M%Y(J)   ; BC%Z1=M%Z(K)   ; BC%Z2=M%Z(K)
 END SELECT
 
+! Determine "free stream" cell indices, IIS, JJS, KKS, for use in empirical boundary layer models
+
+IF (SIM_MODE==DNS_MODE) THEN
+   BC%IIS = BC%IIG
+   BC%JJS = BC%JJG
+   BC%KKS = BC%KKG
+ELSE
+   CALL COMPUTE_BOUNDARY_CLEARANCE(NM,BC%IOR,BC%IIG,BC%JJG,BC%KKG,CLEARANCE)
+   CLEARANCE = MIN(0.5_EB*CLEARANCE,SF%BOUNDARY_LAYER_THICKNESS) - 0.0001_EB
+   SELECT CASE(BC%IOR)
+      CASE(-1) ; CALL GET_IJK(BC%X-CLEARANCE,BC%Y,BC%Z,NM,XI,YJ,ZK,BC%IIS,BC%JJS,BC%KKS)
+      CASE( 1) ; CALL GET_IJK(BC%X+CLEARANCE,BC%Y,BC%Z,NM,XI,YJ,ZK,BC%IIS,BC%JJS,BC%KKS)
+      CASE(-2) ; CALL GET_IJK(BC%X,BC%Y-CLEARANCE,BC%Z,NM,XI,YJ,ZK,BC%IIS,BC%JJS,BC%KKS)
+      CASE( 2) ; CALL GET_IJK(BC%X,BC%Y+CLEARANCE,BC%Z,NM,XI,YJ,ZK,BC%IIS,BC%JJS,BC%KKS)
+      CASE(-3) ; CALL GET_IJK(BC%X,BC%Y,BC%Z-CLEARANCE,NM,XI,YJ,ZK,BC%IIS,BC%JJS,BC%KKS)
+      CASE( 3) ; CALL GET_IJK(BC%X,BC%Y,BC%Z+CLEARANCE,NM,XI,YJ,ZK,BC%IIS,BC%JJS,BC%KKS)
+   END SELECT
+ENDIF
+
 ! If the WALL cell is attached to a THIN obstruction, use the obstruction coords for the wall cell coordinates
 
 IF (OBST_INDEX>0) THEN
@@ -3626,6 +3645,9 @@ SUBROUTINE INIT_THIN_WALL_CELL(NM,I,J,K,OBST_INDEX,ITW,IOR,SURF_INDEX,IEC)
 BC%IIG = I
 BC%JJG = J
 BC%KKG = K
+BC%IIS = I
+BC%JJS = J
+BC%KKS = K
 BC%IOR = IOR
 SELECT CASE(BC%IOR)
    CASE( 1) ; BC%NVEC=(/ 1._EB, 0._EB, 0._EB/)

Source/type.f90

@@ -191,6 +191,9 @@ MODULE TYPES
    INTEGER :: IIG            !< Gas cell x index
    INTEGER :: JJG            !< Gas cell y index
    INTEGER :: KKG            !< Gas cell z index
+   INTEGER :: IIS            !< Free stream cell x index
+   INTEGER :: JJS            !< Free stream cell y index
+   INTEGER :: KKS            !< Free stream cell z index
    INTEGER :: IOR=0          !< Index of orientation of the WALL cell
 
    REAL(EB) :: X             !< x coordinate of boundary cell center
@@ -312,13 +315,18 @@ MODULE TYPES
    REAL(EB) :: TMP_G                 !< Near-surface gas temperature (K)
    REAL(EB) :: TMP_F_OLD             !< Holding value for surface temperature (K)
    REAL(EB) :: TMP_B                 !< Back surface temperature (K)
+   REAL(EB) :: TMP_FS                !< Free stream temperature (K)
+   REAL(EB) :: RHO_FS                !< Free stream density (kg/m3)
+   REAL(EB) :: U_TANG_FS             !< Free stream velocity (m/s)
    REAL(EB) :: U_NORMAL=0._EB        !< Normal component of velocity (m/s) at surface, start of time step
    REAL(EB) :: U_NORMAL_S=0._EB      !< Estimated normal component of velocity (m/s) at next time step
    REAL(EB) :: U_NORMAL_0=0._EB      !< Initial or specified normal component of velocity (m/s) at surface
    REAL(EB) :: U_TANG=0._EB          !< Tangential velocity (m/s) near surface
    REAL(EB) :: U_IMPACT=0._EB        !< Impact velocity from stagnation pressure
    REAL(EB) :: RHO_F                 !< Gas density at the wall (kg/m3)
    REAL(EB) :: RHO_G                 !< Gas density in near wall cell (kg/m3)
+   REAL(EB) :: MU_G                  !< Gas viscosity near wall cell (kg/m/s)
+   REAL(EB) :: RSUM_G                !< R/W near wall cell
    REAL(EB) :: RDN=1._EB             !< \f$ 1/ \delta n \f$ at the surface (1/m)
    REAL(EB) :: K_G=0.025_EB          !< Thermal conductivity of gas in adjacent gas phase cell near wall
    REAL(EB) :: Q_DOT_G_PP=0._EB      !< Heat release rate per unit area (W/m2)
@@ -909,6 +917,7 @@ MODULE TYPES
    REAL(EB) :: Y_S_EFF                                 !< Effective soot yield for S_pyro
    REAL(EB) :: TIME_STEP_FACTOR=10._EB                 !< Maximum amount to reduce solid phase conduction time step
    REAL(EB) :: REMESH_RATIO=0.05                       !< Fraction change in wall node DX to trigger a remesh
+   REAL(EB) :: BOUNDARY_LAYER_THICKNESS=0.051_EB       !< Nominal thickness (m) of the boundary layer
 
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: DX,RDX,RDXN,X_S,DX_WGT,MF_FRAC,PARTICLE_INSERT_CLOCK
    REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: RHO_0