MATC enthalpy def for phase change breaks latent heat check

Numerical methods and mathematical models of Elmer
NJank
Posts: 99
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

MATC enthalpy def for phase change breaks latent heat check

Post by NJank »

so "Check Latent Heat Release = true" turns on an adaptive time stepping routine that will reduce the timestep by half and then reruns the timestep if it detects that the phase front has moved farther than a single element. If I use a table to define enthalpy, this works well. It runs 50 timesteps in ~12 seconds, and I see the latent heat release check reducing a number of the earlier timesteps. i.e.,:

Code: Select all

Enthalpy = Variable Temperature
  Real
   0.0     0.0
   $T1     E1
   $T2     E2
   1000    $E1000
  End
Phase Change Intervals(2,1) = $T1 T2$
however, if I want something more complex than a linear enthalpy model (maybe a polynomial for 'smooth' behavior), I need to either specify a ton of points to do a piecewise linear approximation, or use a MATC definition of the polynomial. to recreate the interval, however, I have to use if statements within the MATC definition. (documentation didn't indicate the availability of a step or heaviside function). this somehow breaks the adaptive routine. i.e., with the following code:

Code: Select all

 Enthalpy = Variable Temperature
    MATC "if (tx<=T1) ((E1/T1)*tx); else if (tx>=T2) (((E1000-E2)/(1000-T2))*(tx-T2)+E2); else (((E2-E1)/(T2-T1))*(tx-T1)+(E1));"
  Phase Change Intervals(2,1) = $T1 T2$
which should recreate exactly the linear profile from the table above, the adaptive routine goes into an endless loop of continually repeating and decreasing the timestep time. nominally I had it set to 1000s, and with the table it reduced some initial timesteps all the way down to maybe 50 or so. but that was only the initial few, and then as melting slows it used larger timesteps and eventually the default 1000. With the adaptive timestep, it never gets past the first timestep (let it run for hours when I first saw this behavior, timesteps in the microseconds when I looked in on it.)

I attached simplified cases grd and sif files (they're for solving the neumann 1d melting problem). am I misusing the if statements in the MATC somehow? is there a better way to do this? both enthalpy definitions are in the file, one commented out.

*note: verified that without adaptive timestep, both enthalpy forms produce very similar (but not numerically identical) results. MATC is a lot slower, but the form seems valid.

EDIT: corrected typo in attached sif file, added note.
Attachments
Model.sif
sif file, both enthalpy definitions, one commented with !
(2.5 KiB) Downloaded 555 times
Model.grd
1d bar, grd file
(539 Bytes) Downloaded 481 times
NJank
Posts: 99
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

Re: MATC enthalpy def for phase change breaks latent heat check

Post by NJank »

playing around with this more, tried making a UDF to do the same thing. here's the modified Materal section with all 3 options

Code: Select all

Material 1
   Name = "PCM"
   Density = $rhos$
   Heat Conductivity = $kths$
   
!simple linear enthalpy change
  
!Enthalpy = Variable Temperature
!  Real
!   0.0     0.0
!   $T1     E1$
!   $T2     E2$
!   1000    $E1000$
! End    

!Enthalpy = Variable Temperature
!    MATC "if (tx<=T1) ((E1/T1)*tx); else if (tx>=T2) (((E1000-E2)/(1000-T2))*(tx-T2)+E2); else (((E2-E1)/(T2-T1))*(tx-T1)+(E1));"

 ENTH1 = Real $E1$
 ENTH2 = REAL $E2$
 ENTH1000 = REAL $E1000$
 TEMP1 = REAL $T1$
 TEMP2 = REAL $T2$

Enthalpy = Variable Temperature
 Real Procedure "EnthCalc" "LinEnthCalc"

Phase Change Intervals(2,1) = $T1 T2$
End
and the UDF is below. After all the program bloat, it's basically just the same if-else structure used in the MATC code.

EnthCalc.f90:

Code: Select all

FUNCTION LinEnthCalc(model,n,temp) RESULT(enthout)
USE DefUtils
IMPLICIT None

!variables in function header
TYPE(Model_t) :: model
INTEGER :: n
REAL(KIND=dp) :: temp, enthout

!variables needed inside function
REAL(KIND=dp) :: T1, T2, enth1, enth2, enth1000
Logical :: GotIt
TYPE(ValueList_t), POINTER :: material

! get pointer on list for material
material => GetMaterial()
IF (.NOT. ASSOCIATED(material)) THEN
		  CALL Fatal('LinEnthCalc', 'No material found')
END IF

!read in reference temperatures
T1 = GetConstReal( material, 'TEMP1', GotIt)
IF(.NOT. GotIt) THEN
		  CALL Fatal('LinEnthCalc', 'TEMP1 not found')
END IF
T2 = GetConstReal( material, 'TEMP2', GotIt)
IF(.NOT. GotIt) THEN
		  CALL Fatal('LinEnthCalc', 'TEMP2 not found')
END IF

! read in reference enthalpies
enth1 = GetConstReal( material, 'ENTH1', GotIt)
IF(.NOT. GotIt) THEN
		  CALL Fatal('LinEnthCalc', 'ENTH1 not found')
END IF

enth2 = GetConstReal( material, 'ENTH2', GotIt)
IF(.NOT. GotIt) THEN
		  CALL Fatal('LinEnthCalc', 'ENTH2 not found')
END IF

enth1000 = GetConstReal( material, 'ENTH1000', GotIt)
IF(.NOT. GotIt) THEN
		  CALL Fatal('LinEnthCalc', 'ENTH1000 not found')
END IF

!caclulate actual enthalpy
if (temp .LE. T1) then
		  enthout = temp * enth1 / T1
else if (temp .GE. T2) then
		  enthout = ((enth1000-enth2)/(1000.0-T2))*(temp-T2)+enth2
else
		  enthout = ((enth2-enth1)/(T2-T1))*(temp-T1)+enth1
end if

END FUNCTION LinEnthCalc
Running with the UDF produces the same results as the MATC. It works without the adaptive Check Latent Heat Release, but with it turned on it gets stuck in a loop. Just for kicks, I timed execution for the different methods:

table no adaptivity= 9.45s
table with adaptivity = 14.13s
matc no adaptivity= 107.76sec
udf no adaptivity = 32.71sec
matc and udf with adaptivity = infinite
NJank
Posts: 99
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

Re: MATC enthalpy def for phase change breaks latent heat check

Post by NJank »

some more info.

Trying to debug, I grabbed the HeatSolve.scr and inserted a Pause at 1191 inside the IF after checking for the condition to reduce the timestep.

I had it print out the nodal temperatures and the temps from the previous time step. When using the adaptive stepping, and the enthalpy table, the last 4 'nodes of interest' produces the following after the first iteration:

Code: Select all

110.42895584970664        117.80300116068474        130.39104344623897        150.00000000000003
100.000000000000000       100.000000000000000       100.000000000000000        150.00000000000000
because the delta on the 2nd to last node stepped completely over the phase change interval (117.75-118.25), it reduces dt and reruns. that 2nd to last node gets closer together after a few iterations, and eventually it moves on.

if I use either the MATC or UDF, though, I get the following:

Code: Select all

149.84999999998834        149.89999999999222        149.94999999999609        149.99999999999997
100.000000000000000       100.000000000000000       100.000000000000000        150.00000000000000
something pulls almost all the nodes immediately up to T_hot, and they stay there no matter how many iterations. seems like a 'zero' heat capacity condition or something. Thus, it goes into an infinite loop. so, now just need to find out why it works fine without the Check Latent Heat release, but this happens with it turned on.

Also, just tried and with MATC/UDF, even with Check Latent Heat Release = False, the first time step has the temperatures up at 149.999..., i.e., heat capacity of zero. It starts correcting after the first time step. Will dig into how it's pulling enthalpy/heat capacity info and follow up.
NJank
Posts: 99
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

Re: MATC enthalpy def for phase change breaks latent heat check

Post by NJank »

ok, put a bunch of breaks in at the various points it pulls enthalpy/heat capacity. it appears the problem is within the EffectiveHeatCapacity subroutine.

the first time through the transient model, because there is no previous timestep, even with Phase Change Model = 'Temporal' it instead uses 'Spatial 1' to calcualte the initial material enthalpy. I.e., in Line 1594, TSolution is not yet allocated, so it executes the ELSE PhaseChangeModel= PHASE_SPATIAL_1 in line 1604.

This affects how it pulls the enthalpy from the sif file definitions. In Spatial 1, first at line 1635 it checks to see if an "Effective Heat Capacity" is defined. It's not, so then it checks at 1638 to see if properties are defined per mass (specific) or per volume. Mine isn't specific, so it goes on to execute 1643 to get a heat capacity:

Code: Select all

  HeatCapacity(1:n) = ListGetDerivValue( Material, &
                'Enthalpy', n,Element % NodeIndexes )      
Now, when Enthalpy = Variable Temperature is defined with a table, at this point Heat Capacity is populated with legitimate values. It finishes all elements for this first iteration. then, for the second iteration, because there's a previous timestep, it sets PhaseChangeModel =PHASE_TEMPORAL, and from then on it uses lines 1658-1676 to determine HeatCapacity, which apparently works fine.

Doing nothing more than switching to the MATC definition, however, it never gets any value of HeatCapacity during the first iteration when it goes through PHASE_SPATIAL_1. It goes through the same program path, but I guess that ListGetDerivValue function fails in some way when it's not a table. I noticed there's no logical 'GotIt" check on the function, so it carries on desipite getting nothing.

If CheckLatentHeatRelease = False, then it will go through with bad values to the second iteration, and at the second iteration it will go through the PHASE_TEMPORAL structure, and from then on it works correctly. HOWEVER, if you set "CheckLatentHeatRelease = True", it sees the large error, reduces the timestep, and reruns the first iteration. Thus it goes through PHASE_SPATIAL_1 again, HeatCapacity stays at zeros, it checks the large error, reduces the timestep, reruns the first iteration... etc., etc. ad infinitum.

SO, THERE's THE PROBLEM! I'll start looking to tweak it to see if I can find a fix.
NJank
Posts: 99
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

Re: MATC enthalpy def for phase change breaks latent heat check

Post by NJank »

I'll just keep a running commentary here:

tracing back the function call, the main problem I see:

Code: Select all

!------------------------------------------------------------------------------
!> Gets a real derivative from. This is only available for tables with dependencies.
!------------------------------------------------------------------------------
   RECURSIVE FUNCTION ListGetDerivValue(List,Name,N,NodeIndexes) RESULT(F)
so... this call only works for tables. I guess that explains the error. If that can't be changed, would be nice to have it throw an error.

EDIT:
- looking at how the other cases grab the enthalpy, it appears that the SPATIAL_2 case grabs Enthalpy, but never actually makes an assignment to HeatCapacity? Looks like that enthalpy assignment is never used either?

- looking at ListGetDerivValue, it looks like it could be expanded to other types of data. currently the actual math is done under a SELECT CASE, where the only case listed is: CASE( LIST_TYPE_VARIABLE_SCALAR )
not sure what other CASE TYPES are, but will keep digging. Again, this would be a good place to add a DEFAULT: CALL FATAL(list type not defined), assuming there's nowhere that it is used expecting something to pass through the case without execution.
NJank
Posts: 99
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

Re: MATC enthalpy def for phase change breaks latent heat check

Post by NJank »

since the Temporal part works fine, the problem is getting past the Spatial_1 part. This part needs a local derivative (dH/dT), which I guess is hard to guess at numerically for data that isn't definitely piecewise linear. Hence the ListGetDerivValue calls DerivateCurve which uses the points in a table to get the linear slope, which it gives back as the derivative. For anything nonlinear you'd have to guess at the delta-T to use. I guess you could make an iterative convergence (start with some small value epsilon, find H's, divide eps by 2, find new H's. set a convergence limit on the deltaH/deltaT. but that would be another arbitray value to set, and would probably need a keyword assigned if it was to be adjustable for arbitrary poorly behaved functions.

the quick fix, I didn't realize there was an "Effective Heat Capacity" keyword. it's only ever called once in Spatial_1, I don't see any other reference to it anywhere. So, if I defined that, it would get used for the first step of a transient sim with "Temporal" set. It would also get used any other time Spatial_1 gets called, but that should never happen after the first iteration for the transient case. Effective Heat Capacity is retrieved using a ListGetReal, which should work fine for either a tabular, MATC or UDF definition.

Actually, you probably _shouldn't_ use a tabular definition for effective heat capacity if you're using a linear enthalpy model. The table would force you to lose the jump discontinuity over phase change, which may deviate from your intended linear enthalpy model.

Will test this with MATC and UDF and post results.
NJank
Posts: 99
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

Re: MATC enthalpy def for phase change breaks latent heat check

Post by NJank »

end of the narrative:
using either a MATC or UDF defined "Effective Heat Capacity" fixes the problem. It allows adaptive time stepping to be used with arbitrary enthalpy profiles. Running the linear test as above all three ways both with and without "Check Latent Heat Release" worked fine, and the three methods produces nearly equal results. for some reason the matc version was different from the others by a fraction of a degree.

So, to do an effective capacity phase change model with Check Latent Heat Release you can:

A) use a linear enthalpy model defined in a table

B) use MATC or a UDF to define an arbitrary enthalpy function, using if-else statements to capture any discontinuities in the the solid-mushy-liquid enthalpy profile, AND provide an "Effective Heat Capacity" function defined as the derivative of the enthalpy function (can either be MATC or UDF).

also, apparently "Effective Heat Capacity" is not in the keywords list? I get a warning on that whenever the solver runs. wolud also recommend either Warnings or Errors get thrown whenever it tries to run something that only works on a table, and you haven't provided a table.

Finally, just for fun I noted the times for 50 1000s timesteps and 1000 1D elements:

table not adaptive 9.91
table adaptive 7.24
matc not adaptive 88.57
matc adaptive 79.78
udf not adaptive 11.26
udf adaptive 8.21

I recommend the UDF.

Also here's the Material section of the SIF file containing all 3 options (all but 1 commented out) and I attached the EnthCalc.f90 file which includes both the LinEnthCalc and LinEffCapCalc UDFs.

Code: Select all

Material 1
   Name = "PCM"
   Density = $rhos$
   Heat Conductivity = $kths$
   Phase Change Intervals(2,1) = $T1 T2$   

!simple linear enthalpy change  
Enthalpy = Variable Temperature
  Real
  0.0     0.0
  $T1     E1$
   $T2     E2$
   1000    $E1000$
 End    

!Enthalpy = Variable Temperature
!	Real MATC "if (tx<=T1) ((E1/T1)*tx); else if (tx>=T2) (((E1000-E2)/(1000-T2))*(tx-T2)+E2); else (((E2-E1)/(T2-T1))*(tx-T1)+(E1));"

!Effective Heat Capacity = Variable Temperature
!  Real MATC "if (tx<=T1) (E1/T1); else if (tx>=T2) (E1000-E2)/(1000-T2); else (E2-E1)/(T2-T1);"

! ENTH1 = Real $E1$
! ENTH2 = REAL $E2$
! ENTH1000 = REAL $E1000$
! TEMP1 = REAL $T1$
! TEMP2 = REAL $T2$

!Enthalpy = Variable Temperature
! Real Procedure "EnthCalc" "LinEnthCalc"
 
!Effective Heat Capacity = Variable Temperature
! Real Procedure "EnthCalc" "LinEffCapCalc"

End
Attachments
EnthCalc.f90
User defined functions for linear enthalpy model
(2.9 KiB) Downloaded 511 times
raback
Site Admin
Posts: 4812
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: MATC enthalpy def for phase change breaks latent heat check

Post by raback »

Hi Njank

And that you for your diligent work. I just added 'Effective Heat Capacity' in the keyword database and some error detection to the ListGetDerivValue. Do you any other proposals on what should be altered? Basically the derivative could be computed from the other dependency types as well using numerical derivation but maybe there is no dire need for that.

The timings are illuminating. Of course in 1D case the linear system is tridiagonal and therefore extremely cheap to solve and therefore the time used to evaluate the enthalphy with the MATC function dominates. But even generally MATC evaluations parameters in the bulk elements are better used only if speed is not crictical. The boundary conditions are usually not that critical.

-Peter
NJank
Posts: 99
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

Re: MATC enthalpy def for phase change breaks latent heat check

Post by NJank »

Was it just Lists.src that was changed? i'll pull and recompile that and any other from the trunk so I don't have to pull the whole fresh install (running off windows).

So right now EffectiveHeatCapacity (spatial 1) only works for an Enthalpy table or an "Effective Heat Capacity". Transient works, except that the first step is handled by spatial 1, so the same rules apply. And I'm not sure if Spatial 2 actually works at all. Looks like it's part way there, but stops at getting Enthalpy. Probably should have some placeholder error-messages for now until that gets filled in. Seems that stops short of handling the derivative, too, but its also a bit more complicated. other messages would include any time PHASE_TEMPORAL is true, it would check for the existance of both EHC and Enthalpy, and throw a warning if both aren't there. I think a while back I mentioned something about Check Latent Heat Release and Adaptive Time Stepping throwing an error if you don't have Simulation / BDF Order = 1.

using ListGetDerivValue works very efficiently for a table, so wouldn't want to change that. is there any easy way of telling whether it's a table? seems the only type check involves looking for LIST_TYPE_VARIABLE_SCALAR, and that doesn't differentiate source type in this way. Are there other functions that call ListGetDerivValue? Wouldn't want to break anything for them. having to have both Enthalpy and Effective HC is a bit redundant, would be nice to work with just Enthalpy. you'd know better than I what the best way would be to get a numerical derivative from an arbitrary well behaved function.

don't know if Spatial 2 could pull from same derivative function. looks like it needs gradH and gradT. gradT should be straightforward, I assume there are other models using this, but it would have to reach beyond the local element solution. same for H. Looks nontrivial, but I guess that would wind up just use the ListGetReal function, though.
NJank
Posts: 99
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

Re: MATC enthalpy def for phase change breaks latent heat check

Post by NJank »

NJank wrote:Was it just Lists.src that was changed? i'll pull and recompile that and any other from the trunk so I don't have to pull the whole fresh install (running off windows).
or not. any tips for those of us relying on the elmerf90 script on recompiling a mod file as opposed to just a solver dll? hoping its just a different gfortran switch...?
Post Reply