SWELL_RATIO computation

[drjfloyd]

I don’t think this block of code in read.f90 (~9410) is correct.

      DENSITY_MAX = 0._EB
      DENSITY_MIN = 10000000._EB
      DO N = 1,SF%N_MATL
         ML => MATERIAL(SF%MATL_INDEX(N))
         DO NN = 1,SF%N_LAYER_MATL(NL)
            IF ((ML%PYROLYSIS_MODEL==PYROLYSIS_SOLID.OR.ML%PYROLYSIS_MODEL==PYROLYSIS_SURFACE_OXIDATION) .AND. &
                 SF%LAYER_MATL_INDEX(NL,NN)==SF%MATL_INDEX(N)) THEN
               DENSITY_MAX = MAX(DENSITY_MAX,SF%MATL_MASS_FRACTION(NL,NN)*SF%LAYER_DENSITY(NL))
            ENDIF
         ENDDO
         DENSITY_MIN = MIN(DENSITY_MIN,ML%RHO_S)
      ENDDO

      IF (SF%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED .AND. DENSITY_MIN>TWO_EPSILON_EB) &
         SF%SWELL_RATIO = MAX(1._EB,DENSITY_MAX/DENSITY_MIN)

We are in a loop over material layers and for each layer we compute and store SF%SWELL_RATIO. This gets used later in init.f90. If there are multiple layers, then we could be using a different value in init.f90 for a layer than used in read.f90 since SF%SWELL_RATION in init.f90 would be that of the last layer. I think this variable should be a vector over layers that is initialized to 1 (no swelling).

I also think we don’t get the correct values for the maximum and minimum density going into it. Swelling can only happen if it is allowed (ML%ALLOW_SWELLING=T), SF%PYROLYSIS_MODEL=PYROLYSIS_PREDICTED, and there is some residue MATL possible in some reaction associated with the surface. If any is false, then the surface could not swell in wall.f90. ML%PYROLYSIS_MODEL isn’t used in wall.f90 for remeshing decisions. Additionally, we could have some complicated set of MATL reactions in a layer where the initial composition is multiple MATL with multiple reactions that become other MATL that have multiple reactions. Trying to untangle the exact worst-case swelling might not be possible; however, the absolute worst case maximum swelling would be if the initial layer composition was 100 % converted to the least dense material in the chain of reactions. I think this should be:

      IF (SF%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) THEN
         DENSITY_MIN = 10000000._EB
         DO N = 1,SF%N_MATL
            ML => MATERIAL(SF%MATL_INDEX(N))
            IF (ANY(ML%NU_RESIDUE>0._EB)) CAN_SWELL = .TRUE.
           IF (ML%ALLOW_SWELLING) DENSITY_MIN = MIN(DENSITY_MIN,ML%RHO_S)
         ENDDO
         IF (DENSITY_MIN>TWO_EPSILON_EB .AND. CAN_SWELL) SF%SWELL_RATIO(NL) = MAX(1._EB,SF%LAYER_DENSITY(NL)/DENSITY_MIN)
      ENDIF

[drjfloyd]

@mcgratta

For HT3D or HT1D with VARIABLE_THICKNESS the SURF may not have a MATL, only the OBST. If no MATL is assigned to the SURF, the SF%SWELL_RATIO will be 1 for the call to GET_N_LAYER_CELLS in init.f90 at 1667. Do we allow for an OBST to be assigned a MATL that has some reaction scheme that leads eventually to swelling? If so, then I we either need some kind of OBST SWELL_RATIO to deal with this or only allow that if the SURF has a MATL otherwise we might not allocate enough storage space.

[mcgratta]

I'll take a look.

[mcgratta]

I need to make SF%SWELL_RATIO an array and then create a ONE_D%SWELL_RATIO as well. Will that cause a conflict for you?

[drjfloyd]

I have already made SWELL_RATIO a layer array in my local repo.

[drjfloyd]

I made a wall_stuff branch and pushed up my current source. The changes in SUBROUTINE SOLID_HEAT_TRANSFER in wall.f90 are a work in progress (Jonathan's issue on runtime). The changes in PYROLYSIS seem OK so far.

[mcgratta]

Here's what you need to do:

1. Create ONE_D%SWELL_RATIO. Use ONE_D%LAYER_THICKNESS as a template for setting it up.
2. In init.f90, add the machinery to calculate ONE_D%SWELL_RATIO to the loop around 4060 starting with DO NN=1,ONE_D%N_MATL
3. Replace SF%SWELL_RATIO with ONE_D%SWELL_RATIO in the call to GET_N_LAYER_CELLS near init.f90 1667.

[drjfloyd]

Thanks for the confirmation that ONE_D%SWELL_RATIO is the way to go.

[drjfloyd]

I changed things so that in READ_MATL we create an integer MATERIAL array that tracks which materials are present due to reactions -- any residues along the reaction chain and the original material. The arrray is 1 if the material is present and zero other wise. The for any surface layer or HT3D obstruction we can simply sum the arrays for the initial materials and do a WHERE > 1 set to 1. That quickly gives the number of materials (sum of the new array) and which materials are present. Then we can repeat over all layers for a master list of all materials in surface. Now SWELL_RATIO is set as the initial density divided by the minimum material density. In cases with multiple residues or multiple reactions this might overestimate the SWELL_RATIO, if there were a material with multiple reactions leading to other materials with multple reactions, there wasn't an obvious clean way to determine how to pick a better ratio. This passed firebot but not sure if we have enough HT3D and variable thickness pyrolysis cases in Verifcation for this.

[shostikk]

logic seems good. Did you run the swell/shrink verification tests?

[drjfloyd]

Ran firebot and there were no verifcation cases out of tolerance. Made some tweaks for a suggestion by Randy so running firebot again now.