FDS Source: add analytical solutions for certain cases of CONE_MESH_INTERSECTION_VOLUME. Issue #14139

Author: ericvmueller

Source/func.f90

@@ -9436,7 +9436,8 @@ REAL(EB) FUNCTION CONE_MESH_INTERSECTION_VOLUME(NM,X0,Y0,Z0,RR0,RRI,HH0,G_FACTOR
 INTEGER, INTENT(IN) :: NM,G_FACTOR
 REAL(EB), INTENT(IN) :: X0,Y0,Z0,RR0,RRI,HH0
 INTEGER :: I,J,K,NX,NY,NZ
-REAL(EB) :: XX,YY,ZZ,R2,X_MIN,X_MAX,Y_MIN,Y_MAX,Z_MIN,Z_MAX,DX,DY,DZ,SHORT_DIMENSION
+REAL(EB) :: XX,YY,ZZ,R2,X_MIN,X_MAX,Y_MIN,Y_MAX,Z_MIN,Z_MAX,DX,DY,DZ,SHORT_DIMENSION,&
+            OUTER_BASE_AREA,INNER_BASE_AREA,HH_IN,RR0_BASE,RRI_BASE,RR0_TOP,RRI_TOP
 TYPE (MESH_TYPE), POINTER :: M
 
 CONE_MESH_INTERSECTION_VOLUME = 0._EB
@@ -9448,18 +9449,38 @@ REAL(EB) FUNCTION CONE_MESH_INTERSECTION_VOLUME(NM,X0,Y0,Z0,RR0,RRI,HH0,G_FACTOR
 Y_MAX = MIN(M%YF,Y0+RR0)
 Z_MIN = MAX(M%ZS,Z0)
 Z_MAX = MIN(M%ZF,Z0+HH0)
+HH_IN = Z_MAX-Z_MIN
+
+! Case 1: volume fully outside the mesh
 IF (X_MAX<=X_MIN .OR. Y_MAX<=Y_MIN .OR. Z_MAX<=Z_MIN) RETURN
+
+! Case 2: base area fully inside the mesh
 IF ((X_MIN-M%XS)>TWO_EPSILON_EB .AND. (M%XF-X_MAX)>TWO_EPSILON_EB &
-   .AND. (Y_MIN-M%YS)>TWO_EPSILON_EB .AND. (M%YF-Y_MAX)>TWO_EPSILON_EB &
-   .AND. (Z_MIN-M%ZS)>TWO_EPSILON_EB .AND. (M%ZF-Z_MAX)>TWO_EPSILON_EB) THEN
+   .AND. (Y_MIN-M%YS)>TWO_EPSILON_EB .AND. (M%YF-Y_MAX)>TWO_EPSILON_EB) THEN
    IF (G_FACTOR==0) THEN
-      CONE_MESH_INTERSECTION_VOLUME = PI*(RR0**2._EB-RRI**2._EB)*HH0
+      ! Subtract inner from outer cylinder volume
+      CONE_MESH_INTERSECTION_VOLUME = PI*(RR0**2._EB-RRI**2._EB)*HH_IN
    ELSEIF (G_FACTOR==1) THEN
-      CONE_MESH_INTERSECTION_VOLUME = PI*(RR0**2._EB-RRI**2._EB)*HH0/3._EB
+      RR0_BASE = RR0*(1-(Z_MIN-Z0)/HH0)
+      RRI_BASE = RRI*(1-(Z_MIN-Z0)/HH0)
+      RR0_TOP = RR0*(1-(Z_MAX-Z0)/HH0)
+      RRI_TOP = RRI*(1-(Z_MAX-Z0)/HH0)
+      ! Subtract inner from outer truncated cone volume
+      CONE_MESH_INTERSECTION_VOLUME = PI*HH_IN/3._EB*&
+         (RR0_BASE**2._EB+RR0_BASE*RR0_TOP+RR0_TOP**2._EB-RRI_BASE**2._EB-RRI_BASE*RRI_TOP-RRI_TOP**2._EB)
    ENDIF
    RETURN
 ENDIF
 
+! Case 3: partial hollow cylinder intersection volume
+IF (G_FACTOR==0) THEN
+   OUTER_BASE_AREA = CIRCLE_CELL_INTERSECTION_AREA(X0,Y0,RR0,M%XS,M%XF,M%YS,M%YF)
+   INNER_BASE_AREA = CIRCLE_CELL_INTERSECTION_AREA(X0,Y0,RRI,M%XS,M%XF,M%YS,M%YF)
+   CONE_MESH_INTERSECTION_VOLUME = (OUTER_BASE_AREA-INNER_BASE_AREA)*HH_IN
+   RETURN
+ENDIF
+
+! Case 4: partial hollow cone intersection volume (numerical integration)
 SHORT_DIMENSION = MIN(X_MAX-X_MIN,Y_MAX-Y_MIN,Z_MAX-Z_MIN)
 ! Limit resolution length scale 
 SHORT_DIMENSION = MAX(SHORT_DIMENSION, 0.1_EB*MIN(M%DXMIN,M%DYMIN,M%DZMIN))