Modifying Solution-array before time-integration

[markus_h]

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

[Prof. S. Govindjee]

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

[markus_h]

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

[Prof. R.L. Taylor]

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.