FEAP User Forum
FEAP => Programming => Topic started by: markus_h on April 18, 2021, 11:21:18 AM
-
Hello Feap-Community,
I am doing a thermal Simulation and would like to introduce an upper limit to the possible Temperature values at each node (the evaporation temperature of my material, say 1000K)
I want to check the solution-array after solving a timestep and set all temperatures above 1000K to 1000K. (the excess thermal Energy is handeled separately)
So i want so set:
U_mod=min(U_sol,1000K)
after each timestep.
I would like the time integration, so the update of velocity and acceleration, to happen based on the new values of U_mod and not U_sol.
To check if this is the case, I set the Solution-vector to 1.0 after each step in a umacr.
I expected the temperature Rate to be 0.0 because the temperature stays constant at 1.0.
However this was not the case.
My Solution commands are:
PARA
dt=1.2*10^(-9) !timestep
t=dt*8 !tmax time of simulation
BATCh
TRANsient,NEWMark
DT,,dt !set timestep
LOOP INFInite
TIME,,t
DISP all !prints U_mod (all values are 1.0 at every timestep)
VELO all !prints velocities not equal to 0.0, which i would have expected
TANGent,,1
tchk,,100 !calls umacr
NEXT
END
in my umacr i have the code:
integer :: index_U, xlen, xpre,index_d
logical flag
include 'comblk.h' ! hr array
call pgetd ("U", index_U,xlen,xpre,flag)
!set all solutions to 1 to check how it affects time integration
hr(index_U:index_U+xlen-1)=1.0d0
How can i make the time-integration use the modified values of my U_mod?
What command triggers the time-integration? is it at the TIME-Command or does it happen at TANGent?
Kind Regards
-
What happens if you try TRANsient,BACK instead of TRANsient,NEWMark? So using a backward Euler method on your first order ODE system?
Note the implementation of backward Euler is found in program/dyan02.f. program/dyna01.f implements newmark.
Some details on how time integration is handled in FEAP can be found in this report by Karl Steegar: http://faculty.ce.berkeley.edu/sanjay/FEAP/ucb_semm_2010_07.pdf (http://faculty.ce.berkeley.edu/sanjay/FEAP/ucb_semm_2010_07.pdf)
-
Dear Prof Govindjee
Thanks for the report, it really helped make sense of the arrays used in the dyna0X.f routines.
If I use the Euler-Backwarts integration method, the temperature-rate stays constant at 0.0 as expected.
For the Newmark-formula it makes sense that the rate does not go to zero, since the rate of the last timestep is part of the predictor of the rate of the next step.
When I used
TRANsient,NEWMark 0.5 0.5 1
so that gamma and beta have the same values i could also achieve a rate of zero.
I also found out how to use the modified solution-array for the time-integration
I am attaching my code for anyone else who may want to do the same thing.
If i am not mistaken, the call
call dynlib(hr(np(40)),hr(np(42)),2)
updates the time integration of the current step based on values in the U-array, np(40).
Kind regards
-
One thing you need to note is that all integers accessing either the hr(*) or mr(*) array are "integer (kind=8)" or "integer*8", otherwise you will get erratic behavior.
Check the include "p_point.h" where one is defined to use if your wish.