## MATC enthalpy def for phase change breaks latent heat check

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

### MATC enthalpy def for phase change breaks latent heat check

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 !
Model.grd
1d bar, grd file
NJank
Posts: 98
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

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

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

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

T1 = GetConstReal( material, 'TEMP1', GotIt)
IF(.NOT. GotIt) THEN
END IF
T2 = GetConstReal( material, 'TEMP2', GotIt)
IF(.NOT. GotIt) THEN
END IF

enth1 = GetConstReal( material, 'ENTH1', GotIt)
IF(.NOT. GotIt) THEN
END IF

enth2 = GetConstReal( material, 'ENTH2', GotIt)
IF(.NOT. GotIt) THEN
END IF

enth1000 = GetConstReal( material, 'ENTH1000', GotIt)
IF(.NOT. GotIt) THEN
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:

matc and udf with adaptivity = infinite
NJank
Posts: 98
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

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

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: 98
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

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

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: 98
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

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

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: 98
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

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

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: 98
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

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

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:

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
raback
Posts: 3886
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

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

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: 98
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

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

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: 98
Joined: 05 Dec 2009, 00:05
Location: Baltimore, MD, USA

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

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...?