Oscillations in transient thermal results THER_NON_LIN

rodrigofarias · 2021-08-09 20:12:08

Hi.
I've been translating one of my Ansys models to Aster for some time now, and I already have interesting results.

However, I have differences between the thermal cycles (timeXtemperature) obtained in Aster.

Here's what's equal between models:
- same values for properties (LAMBDA - conductivity and BETA - enthalpy )
- very similar mesh
- 3D thermal nonlinear transient model
- Convection and radiation on external surfaces, with same coeficients
- ramp load
- Moving heat load using the same equation, that depends only on the positions of nodes/elements.
- same timestep (~0.9sec)

Differences:
- Ansys has automatic timesteping, Aster does not for thermal problems.

I attach a figure comparing the Ansys thermal cycle (considered a reference) and the Aster thermal cycle.
And a video for better understanding (youtu.be/OF-SSt0JV90)

Why those 'breaks' on the thermal cycles (circled region), instead of a more smooth cycle like Ansys?
Why am I experiencing temperature fluctuations at the points shown in the figure? Is it a convergence problem, which causes the result to fluctuate? The model is converging with no warnings.

For me, it looks like the temperature decreases too rapidly at the beginning of each loadstep, before the next ramp load "kicks in".

In the figure, for both software's, the temperature points are for each substep of each loadstep.

I can provide part of the code if necessary. Thanks for any insight!

Attachments:
fig_ansys_aster.png, Size: 26.93 KiB, Downloads: 153

RichardS · 2021-08-10 07:19:48

Hi Rodrigo,
while we don't have your input data, I can just say that oscillations are a typical behavior for thermal transient analyses with rapid change of temperatures when using standard elements.
You can use diagonalized thermal mass elements (3D_DIAG) to reduce these oscillations.
You can see a good example here on Slide 8: code-aster.org/V2/UPLOAD/DOC/Formations/06-thermal_analysis.pdf

Best,
Richard

rodrigofarias · 2021-08-10 20:39:00

RichardS wrote:

Hi Rodrigo,
while we don't have your input data, I can just say that oscillations are a typical behavior for thermal transient analyses with rapid change of temperatures when using standard elements.
You can use diagonalized thermal mass elements (3D_DIAG) to reduce these oscillations.
You can see a good example here on Slide 8: code-aster.org/V2/UPLOAD/DOC/Formations/06-thermal_analysis.pdf
Best,
Richard

Thanks for the reply Richard!

I tested with 3D_diag element (attached fig) but the behavior is almost the same, different from what happens in the pdf you pointed out. Just changed

model = AFFE_MODELE(AFFE=_F(MODELISATION=('3D', ),

to

model = AFFE_MODELE(AFFE=_F(MODELISATION=('3D_DIAG', ),

Should I have changed something in Salome's mesh generation, or something else in Aster's code?

I never had this type of oscillations in my Ansys simulations (attached fig). The element used in Ansys is the SOLID70
(www .mm.bme.hu/~gyebro/files/ans_help_v182/ans_elem/Hlp_E_SOLID70)
Its a standard 3D linear element for conduction heat transfer, with very similar setup as necessary for the Aster simulation.
My moving heat load is inside a loop, as follows (omitted some parts):

mater = DEFI_MATERIAU(THER_NL=_F(BETA=ENTHAL, LAMBDA=CONDUC))

mesh = LIRE_MAILLAGE(FORMAT='MED', UNITE=3)

model = AFFE_MODELE(AFFE=_F(MODELISATION=('3D', ),
	                         PHENOMENE='THERMIQUE',
	                         TOUT='OUI'),
	                         MAILLAGE=mesh)

for i in range(...):

    X0=hs_centers[i]

    eqf[i] = FORMULE( NOM_PARA=('X','Y','Z', ),
                  VALE=...)

    eqr[i] = FORMULE( NOM_PARA=('X','Y','Z', ),
                  VALE=...)

    times[i] = DEFI_LIST_REEL(DEBUT=tatual,
                              INTERVALLE=_F(JUSQU_A=tatual+dt,
                                            PAS=timestep,),)

    lt[i] = DEFI_LIST_INST(DEFI_LIST=_F(LIST_INST=times[i],
                                        PAS_MINI=0.001,),
                                        METHODE='AUTO')

    loadf[i] = AFFE_CHAR_THER_F(MODELE=model,
                                SOURCE=_F(GROUP_MA=(hsf_groups[i], ),
                                SOUR=eqf[i]),)

    loadr[i] = AFFE_CHAR_THER_F(MODELE=model,
                                SOURCE=_F(GROUP_MA=(hsr_groups[i], ),
                                SOUR=eqr[i]),)

    rter_0 = THER_NON_LINE(reuse=rter_0,
                            CHAM_MATER=fmat,
                            CONVERGENCE=_F(ITER_GLOB_MAXI=50,
                                           RESI_GLOB_RELA=1e-03),
                            ETAT_INIT=(_F(EVOL_THER=rter_0),),
                            EXCIT=(_F(CHARGE=loadf[i],
                                 	  FONC_MULT=RAMP,),
                                   _F(CHARGE=loadr[i],
                                 	  FONC_MULT=RAMP,),
                                   _F(CHARGE=bcs),
                                   _F(CHARGE=linkmesh)),
                            METHODE='NEWTON',
                            INCREMENT=_F(LIST_INST=lt[i]),
                            MODELE=model,
                            NEWTON=_F(REAC_ITER=0, 	
                                   RESI_LINE_RELA=1.E-3,
                                   ITER_LINE_MAXI=my_iter_line,),
                       		COMPORTEMENT=_F(RELATION='THER_NL',),
                       		SOLVEUR=_F(RENUM='AUTO',
                                    NPREC=my_nprec,
                                    ELIM_LAGR='LAGR2',
                                    STOP_SINGULIER='OUI',
                                    TYPE_RESOL='AUTO',
                                    ACCELERATION='AUTO',
                                    LOW_RANK_SEUIL=0.0,
                                    PRETRAITEMENTS='AUTO',
                                    POSTTRAITEMENTS='AUTO',
                                    PCENT_PIVOT=20,
                                    RESI_RELA=-1.0,
                                    GESTION_MEMOIRE='AUTO',
                                    FILTRAGE_MATRICE=-1.0,
                                    MIXER_PRECISION='NON',
                                    MATR_DISTRIBUEE='NON',
                                    METHODE='MUMPS',),
                       		ARCHIVAGE=_F(PRECISION=1.E-06,
                                      CRITERE='RELATIF',),
                       		)
    
    tatual = tatual+dt

IMPR_RESU(RESU=_F(RESULTAT=rter_0),
	          UNITE=80)

Besides the element type, could there be another setup in the code that is generating this type of oscillation?

Again, thanks for the insights!

Attachments:
3D_x_3D_diag.png, Size: 263.97 KiB, Downloads: 156

miib · 2021-08-10 22:56:27

sorry in advance for my replay

for curiosity

can you try ,rebuilding the functions, using this kind of times partition ???

tstep = 0.001

listr = DEFI_LIST_REEL(DEBUT=**start_time**,
INTERVALLE=_F(JUSQU_A=**end_time**,
PAS=tstep))

times = DEFI_LIST_INST(DEFI_LIST=_F(LIST_INST=listr,
PAS_MINI=1e-12),
METHODE='AUTO')

in your case is it possible ???

:-)

mf · 2021-08-11 06:19:39

Hello Rodrigo,

the timestep is also crucial, because the calculation itself is unstable. If Ansys does this automatically, then your two calculations are clearly not the same. Like miib suggests, a very fine timestep and some kind of automation with DEFI_LIST_INST will perhaps do this (you can also use DEFI_LIST_INST to 'observe' a result and cut the timestep although, at the moment, I cannot imagine how to do this in this case....).

Take a look at slide 9 of the before mentioned document (attached). I imagine Ansys does calculate the necessary timestep like in the 2nd formula, try it by hand for your example and see if you are close with your chosen timestep or not. If not, decrease the timestep. I once programmed a diffusion solver for a handheld device (forward difference solver), even in 1D it is similar, you have to choose your timestep wisely (=automatically :-) before the calculation...),

Mario.

EDIT: or try to find out the timestep in Ansys and use it in CA....if all fails..
EDIT2: in principle, you could also do an automated timestep with a little Python formula! :-)

Last edited by mf (2021-08-11 08:16:39)

Attachments:
Bildschirmfoto vom 2021-08-11 07-07-55.png, Size: 164.14 KiB, Downloads: 145

Volker · 2021-08-11 12:23:12

I have often used this construct for thermal calculations of power plant equipment

FORMULE(VALE='Python_function_what_ever_I_wanted(X,Y,Z,inst)',
NOM_PARA=('X','Y','Z','INST'),)

without a for in range loop. I find that the units in Code_Aster are a bit tricky at times.

Kind regards Volker

rodrigofarias · 2021-08-11 17:39:50

miib wrote:

sorry in advance for my replay
for curiosity
can you try ,rebuilding the functions, using this kind of times partition ???

tstep = 0.001
listr = DEFI_LIST_REEL(DEBUT=**start_time**,
INTERVALLE=_F(JUSQU_A=**end_time**,
PAS=tstep))
times = DEFI_LIST_INST(DEFI_LIST=_F(LIST_INST=listr,
PAS_MINI=1e-12),
METHODE='AUTO')

in your case is it possible ???
:-)

Hi! I think its not possible. In the original code I have dt=0.9s and timestep=0.9/50=0.018s. 'tatual' is my current time.

So first loadstep (in the loop) goes from 0 to 0.9 with 0.018 increment, the second from 0.9 to 1.8s, and so on...
If I understood your suggestion, if I set the start_time=0 and end_time= ~85s (end of my simulation), I will end up applying only my first heat load during 85s. So I need to do it in my way to move the heat along the plate.

As for the 'PAS_MINI=1e-12', I could just delete from the code, since Aster does not do automatic timestepping for THER_NON_LINE problems. The command is issued to the software, but doesn't change anything.

If I misunderstood something, please clarified a little bit. Thanks for the suggestion!

rodrigofarias · 2021-08-11 17:56:36

mf wrote:

Hello Rodrigo,
the timestep is also crucial, because the calculation itself is unstable. If Ansys does this automatically, then your two calculations are clearly not the same. Like miib suggests, a very fine timestep and some kind of automation with DEFI_LIST_INST will perhaps do this (you can also use DEFI_LIST_INST to 'observe' a result and cut the timestep although, at the moment, I cannot imagine how to do this in this case....).
Take a look at slide 9 of the before mentioned document (attached). I imagine Ansys does calculate the necessary timestep like in the 2nd formula, try it by hand for your example and see if you are close with your chosen timestep or not. If not, decrease the timestep. I once programmed a diffusion solver for a handheld device (forward difference solver), even in 1D it is similar, you have to choose your timestep wisely (=automatically :-) before the calculation...),
Mario.
EDIT: or try to find out the timestep in Ansys and use it in CA....if all fails..
EDIT2: in principle, you could also do an automated timestep with a little Python formula! :-)

Excellent comment my friend! Since the model is converging nicely, I would not expect my 'timestep' to be a problem.

I thought about using the formula you mentioned, but some variables are not clearly identified in the documentation:
rho is density, C specific heat, theta is PARM_THETA from Aster, lambda is conductivity.
But k_mesh? Is it average element size in the refined region? And delta_x? Is my plate length? Ive seem similar equations for timestepping.

Regarding your EDIT2, I was thinking exactly the same! But if is converging and I set a big timestep that leads to non convergence, how I would restart the simulation from previous point? Nonetheless, is a great idea.

I will test with smaller timestep and see what happens.
Thanks for the infos!

rodrigofarias · 2021-08-11 18:06:21

Volker wrote:

I have often used this construct for thermal calculations of power plant equipment
FORMULE(VALE='Python_function_what_ever_I_wanted(X,Y,Z,inst)',
NOM_PARA=('X','Y','Z','INST'),)
without a for in range loop. I find that the units in Code_Aster are a bit tricky at times.
Kind regards Volker

Hi Volker!

As I remember, Aster didn't allow use of position (X,Y,Z) and time(INST) in the FORMULE command at the same time. Do you use it without problems?

Regarding the function, it would have to be something in the form of a list, with an 'IF' where for that particular loadstep, one of the loads equals to one and the others equals to zero.

I'm going to have to think of some way to try and implement this type of function.
Thanks for the suggestion.

miib · 2021-08-12 14:49:55

rodrigofarias wrote:

But k_mesh? Is it average element size in the refined region? And delta_x? Is my plate length? Ive seem similar equations for timestepping.

boh! :-)

Δx is surely a length

by dimensional analysis k hasn't unit of measure.

may be statistical constant ????

σ±kμ (see attached pic)

substituting this actually hypothetical values , can you obtain a value of Δt with sense ???

:-)

Last edited by miib (2021-08-12 14:55:44)

Attachments:
sigma +- k mu.png, Size: 161.49 KiB, Downloads: 129

miib · 2021-08-12 15:02:16

may be statistical constant ????

σ±kμ (see attached pic)

:-)

Attachments:
by mesh information.png, Size: 161.25 KiB, Downloads: 115

miib · 2021-08-12 15:23:40

dimensional analysis

maybe???

:-)

miib · 2021-08-12 15:24:26

dimensional analysis

maybe???

:-)

Attachments:
IMG_3959.JPG, Size: 138.23 KiB, Downloads: 125

rodrigofarias · 2021-08-12 15:45:52

miib wrote:

boh! :-)
Δx is surely a length
by dimensional analysis k hasn't unit of measure.
may be statistical constant ????
σ±kμ (see attached pic)
substituting this actually hypothetical values , can you obtain a value of Δt with sense ???
:-)

Hi. I checked the Ansys values and the smallest timestep used is the value I'm using equal to 0.018s.
Attached a comparison using 3D and 3D_DIAG with 0.018, and 3D_DIAG with 0.009s. The heat increased a little but the oscillations continue.

Regarding the formula, looking to other formulas its usually (element_size)^2 / diffusivity . So I have no clue if k_mesh is some statistical value. Using the classical formula it gives ~0.018, exactly what I using

You think that decreasing (<0.009) even more will make any difference?

Again, thanks!

Attachments:
fig_forum_half_time.png, Size: 14.73 KiB, Downloads: 136

miib · 2021-08-12 16:05:32

rodrigofarias wrote:

You think that decreasing (<0.009) even more will make any difference?

i don't know, i started only last week with therm_linear , i'm thinking with you about this "problem" in my holiday time, for understand a know how for the future.

In my opinion, if you want , for curiosity, try with 0.001like time step, if the computation's time don't become monstrously long.

:-)

Last edited by miib (2021-08-12 16:07:08)

#1 2021-08-09 20:12:08

Oscillations in transient thermal results THER_NON_LIN

#2 2021-08-10 07:19:48

Re: Oscillations in transient thermal results THER_NON_LIN

#3 2021-08-10 20:39:00

Re: Oscillations in transient thermal results THER_NON_LIN

#4 2021-08-10 22:56:27

Re: Oscillations in transient thermal results THER_NON_LIN

#5 2021-08-11 06:19:39

Re: Oscillations in transient thermal results THER_NON_LIN

#6 2021-08-11 12:23:12

Re: Oscillations in transient thermal results THER_NON_LIN

#7 2021-08-11 17:39:50

Re: Oscillations in transient thermal results THER_NON_LIN

#8 2021-08-11 17:56:36

Re: Oscillations in transient thermal results THER_NON_LIN

#9 2021-08-11 18:06:21

Re: Oscillations in transient thermal results THER_NON_LIN

#10 2021-08-12 14:49:55

Re: Oscillations in transient thermal results THER_NON_LIN

#11 2021-08-12 15:02:16

Re: Oscillations in transient thermal results THER_NON_LIN

#12 2021-08-12 15:23:40

Re: Oscillations in transient thermal results THER_NON_LIN

#13 2021-08-12 15:24:26

Re: Oscillations in transient thermal results THER_NON_LIN

#14 2021-08-12 15:45:52

Re: Oscillations in transient thermal results THER_NON_LIN

#15 2021-08-12 16:05:32

Re: Oscillations in transient thermal results THER_NON_LIN

Board footer