Phase change / switch variable

[XisZero]

Hi Anil,

Unfortunately I didn´t try the time dependence any further. I Fear it will take some weeks until i can finish this, if at all...

Currently i am stuck in getting the time into the UDF (With "variable temperature, time" in SIF). thats fairly easy to do, but it looks like only the actual timestep (with time and temperature) is saved, so i will need to write a file or something, to store and recall all the values.
And: In another thread i asked about the compiler from the newest ELMER GUI release. It creates an exe, which i can rename to xyz.dll and use, but it would be easier to create the dll directly. is there any option?

Any help is appreciated

Greetings,
X=0

[mzenker]

Hi,

maybe you could create a new variable with a value associated to every node in your mesh, and make it contain information about history (if that is what you intend to do).
If the compiler call in ElmerGUI does not fit your needs, you can also do it in a DOS box with the command

Code: Select all

elmerf90 myfunction.f90 -shared -o myfunction.dll

HTH,

Matthias

[XisZero]

Hi Matthias,

Yes, that was what i wanted to do. This far i couldn´t get to the nodal values of temperature due to lack of programming-knowledge of myself
But i thought, if i get the time and temperature via the "variable ..." in the material parameter, that i already have the nodal values. i just didnt manage to store the values (in form of an array, like arr(time, temp) ).
I think i will have to make an array of higher order (arr(time, temp, x,y,z)). But: Will this array be stored or is it deleted when calling the UDF every timestep?

Running the compiler in DOS is an option. It was just really easy in the previous elmer version (cant remember the Rev number, sorry), i clicked compiler, selected the UDF and voila, i got the MyUDF.dll So i could instantly change the UDF and start the test case.

Greetings,

X=0

[mzenker]

Hi X=0,

I see two options:
1. You create an array inside your udf and save it in the udf. You can use the SAVE command in FORTRAN for that.
2. You create a nodal variable in the Elmer structure and store your information there. You can use the Elmer subroutines (like VariableAdd()) for that - you will need to get information on how this is done best from the documentation and the code of existing solvers.

I agree that it would be nice if the compiler command was fixed. This has to be done by the person in the Elmer team taking care of the Windoze installer...

HTH,

Matthias

[annier]

Hi X=0,
How far have you gone in the modeling of undercooling(subcooling) and recrystallization?


Yours Sincerely
Anil Kunwar

[XisZero]

Hi Anil,

Unfortunately I haven´t managed to get any further.
I restarted the whole UDF to some very simple things today. This will take some time to get at least a few results, i could perhaps present.

By the way:
Is it allowed to use Elmer FEM and results thereof in a paper or Thesis? Of course with linking to the project page / appropriate indication.

Greetings,
X=0

[raback]

Hi

To access previous values they are stored for the current timestepping strategy. For higher order schemes you have more values but always in transient case at least the last value. This resides in Var % PrevValues. There are also some Default routines for obtaining the previous values locally...

Elmer is licensed under (L)GPL. The license does not limit the use of code in any way. As developers we are of course happy always when Elmer is used.

-Peter

[XisZero]

Hi everyone,

It took quite a while, but i finally managed to restart my project (again), with some changes:

The working UDF:

Code: Select all

FUNCTION phase(Model, n,  temp) RESULT(ent)
USE DefUtils
IMPLICIT None
TYPE(Model_t) :: Model
TYPE(Solver_t) :: Solver
TYPE(ValueList_t), POINTER :: material
TYPE(Variable_t), POINTER :: TemperatureVariable

Integer :: n,i, ph, Tm, T1, H, T2, T3, T4, Te, DofIndex
REAL(KIND=dp) :: temp, ent, NodalTemperature
Logical, Allocatable :: state(:)
Logical :: ph0 , GotIt, start = .True.

save ph0, Start, state

!Temperature at Nodes

TemperatureVariable => VariableGet(Model % Variables, 'Temperature')
DofIndex = TemperatureVariable % Perm(n)
NodalTemperature = TemperatureVariable % Values(dofIndex)



! Get Material Properties and melting peaks

material => GetMaterial()
IF (.NOT. ASSOCIATED(material)) THEN
CALL Fatal('Subcool', 'No material found')
END IF

T1 = GetConstReal( material, 'T1', GotIt)
IF(.NOT. GotIt) THEN
CALL Fatal('Subcool', 'Parameter T1 not declared')
END IF

Tm = GetConstReal( material, 'Tm', GotIt)
IF(.NOT. GotIt) THEN
CALL Fatal('Subcool', 'Parameter Tm not declared')
END IF

T2 = GetConstReal( material, 'T2', GotIt)
IF(.NOT. GotIt) THEN
CALL Fatal('Subcool', 'Parameter T2 not declared')
END IF

T3 = GetConstReal( material, 'T3', GotIt)
IF(.NOT. GotIt) THEN
CALL Fatal('Subcool', 'Parameter T3 not declared')
END IF

Te = GetConstReal( material, 'Te', GotIt)
IF(.NOT. GotIt) THEN
CALL Fatal('Subcool', 'Parameter Te not declared')
END IF

T4 = GetConstReal( material, 'T4', GotIt)
IF(.NOT. GotIt) THEN
CALL Fatal('Subcool', 'Parameter T4 not declared')
END IF

H = GetConstReal( material, 'H', GotIt)
IF(.NOT. GotIt) THEN
CALL Fatal('Subcool', 'Parameter H not declared')
END IF


!Set Starting conditions

if (temp > T2) then
ph0 = .False.
end if

if (temp < T4 ) then 
ph0= .true.
end if

if (start) then
Allocate(state(10000))
state(n) = ph0
start = .false.
end if


! Check phase of material, Melting

if (state(n)) then

if (NodalTemperature > T1 .and. NodalTemperature < Tm) then 
ent = (H*(NodalTemperature-T1))
end if

if (NodalTemperature >= Tm .and. NodalTemperature <= T2) then
ent = H*(T2-NodalTemperature)
end if

if (NodalTemperature > T2) then
ent = 0
ph0 = .false.
state(n) = ph0
end if

end if


! Freezing

if (.not. state(n)) then

if (NodalTemperature < T3 .and. NodalTemperature > Te) then
ent =  (-H*(T3-NodalTemperature))
end if

if (NodalTemperature <= Te .and. NodalTemperature >= T4) then
ent = (-H *(NodalTemperature-T4))
end if

if (NodalTemperature < T4) then
ent = 0
ph0 =.true.
state(n) = ph0
end if

end if


END FUNCTION phase

I don´t know if everything of the code is really needed, but i just didn´t want to change anything, the UDF runs and does what it is supposed to (and it should stay so ).

here is the relevant part of the SIF:

Code: Select all

...
Material 4
  Name = "test PCM"
  Enthalpy = variable temperature;real procedure "PCM" "phase"
H = real 120e6
T1 = real 348
Tm = real 349
T2 = real 350

T3 = real 347
Te = real 346
T4 = real 345
...

Now, some explanation for those who are interested or have similar problems:

The UDF:
----------------------
First partis to read the Solution of the Heat Solver, the temperature, at every node the solver is applied to.

Second part: read in the melting and freezing temperatures for the material (defined as stated in the SIF code snippet). I assume a triangular shape of the enthalpy function. The triangles base is T2 - T1 long (melting and T3-T4 long (freezing). read the Enthalpy maximum (has to be calculated, in this case: triangle area is the total enthalpy, so the enthalpy maximum point is: 2 * (total enthalpy in J/m³) / (T2-T1).

Third part: The starting conditions. I had several problems with that, but now it works. Allocation of the array to store the state of the material for every node.

Fourth part: The check, if the phase (state) of a node has changed and calculate the enthalpy.
here are the actual triangular functions for the enthalpy and the state variable. the state is set when the nodal temperature rises above T2 (so to say, the material is now all liquid). Now a temperature decrease below the point T3 results in negative enthalpy (energy is released) and the material freezes again.
----------------------

For some testing and parameter changing and evaluation, the code is functional and i can work with it.
By now, the melting and freezing should work only if temperature is higher than T2, and therefore the whole material is in "liquid" state.
What i want to add to the code is an "energy check": The total enthalpy after the "freezing interval" (that is from T4 down to T3) is too high (have to check that). I think that is due to the temperature staying the same for some time. The enthalpy is still calculated with the triangular function, and stays at a high level (I´m trying to explain with an example: Temperature is 300K, the enthalpy is-1e4, due to the enthalpy and energy released by that, the temperature doesn´t decrease. So enthalpy enthalpy, calculated via temperature, stays the same, and more energy is released).
That is physically not correct, so i will have to put in a variable, that checks the total amount of energy that has been stored/released.

i´m not sure when or if i can correct this, but i will post a new UDF if changes are made.

I attached some Pictures made with ParaView from the first simulation run. The setup is simply a heat generation in the middle block fit into a surrounding material, bottom is a cooling water flow, and the top block is the PCM material with the UDF working.
The filenames are the times of the simulation (s for seconds). The heat is generated for a certain time, and then switched off, water and convection is cooling back.
You can see the moving "enthalpy front" for melting in 60s.jpg and 80s.jpg. And the enthalpy of freezing in 240s.jpg (while cooling down).

If you want to use the UDF, feel free to do so. And please let everyone know if you made changes that could be useful for others!

So thanks to all who have contributed to this work!

Oh, and feel free to ask if something is unclear.

Greetings and thanks,

X=0

[annier]

Hi XisZero,
Thank you for your remarkable contribution.

Yours Sincerely
Anil Kunwar