Body temperature rise higher than thermalBndCondValue
Hello,
I'm trying to explore the thermal expansion with Yade, similar to the paper [1]. But while I heated up one side of the boundary: thermal.
###################
from yade import pack, ymport, plot, utils, export, timing
import numpy as np
young=5e6
mn,mx=Vector3(
O.materials.
O.materials.
walls=aabbWalls
wallIds=
sp=pack.
sp.makeCloud(
sp.toSimulation
triax=TriaxialS
maxMultiplier=
finalMaxMultip
thickness = 0,
stressMask = 7,
internalCompac
)
O.engines=[
ForceResetter(),
InsertionSortC
InteractionLoop(
[Ig2_
[Ip2_
[Law2_
),
GlobalStiffnes
triax,
FlowEngine(
ThermalEngine(
NewtonIntegrat
]
O.step()
ss2sc.interacti
is2aabb.
O.dt=0.1e-5
O.dynDt=False
tri_pressure = 10000
triax.goal1=
while 1:
O.run(1000, True)
unb=unbalance
print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
if unb<0.01 and abs(-tri_
break
triax.internalC
for b in O.bodies:
if isinstance(b.shape, Sphere):
b.dynamic=False
flow.dead=0
flow.defToleran
flow.meshUpdate
flow.boundaryUs
flow.useSolver=4
flow.permeabili
flow.viscosity= 0.001
flow.bndCondIsP
flow.bndCondVal
flow.thermalEng
flow.debug=False
flow.fluidBulkM
flow.fluidRho = 997
flow.fluidCp = 4181.7
flow.bndCondIsT
flow.thermalEng
flow.thermalBnd
flow.tZero=333.15
flow.pZero=0
thermal.dead=0
thermal.debug=False
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.fluidK = 0.6069
thermal.
thermal.
thermal.particleT0 = 333.15
thermal.
thermal.particleK = 2.
thermal.particleCp = 710.
thermal.
thermal.
thermal.
timing.reset()
O.dt=1e-5
O.run(1,1)
def bodyByPos(x,y,z):
cBody = O.bodies[1]
cDist = Vector3(
for b in O.bodies:
if isinstance(b.shape, Sphere):
dist = b.state.pos - Vector3(x,y,z)
if np.linalg.
cDist = dist
cBody = b
return cBody
bodyOfInterest = bodyByPos(
from yade import plot
def history():
plot.addData(
ftemp1=
p=flow.
t=O.time,
i = O.iter,
bodyOfIntTemp = O.bodies[
)
plot.saveDataT
O.engines=
plot.plots=
plot.plot(
O.run()
################
Best regards
Jiannan
[1]: https:/
Question information
- Language:
- English Edit question
- Status:
- Solved
- For:
- Yade Edit question
- Assignee:
- No assignee Edit question
- Solved by:
- Jiannan Wang
- Solved:
- 2020-09-09
- Last query:
- 2020-09-09
- Last reply:
- 2020-09-04
This question was reopened
- 2020-08-31 by Jiannan Wang
Robert Caulk (rcaulk) said : | #1 |
Sounds like a stabiltiy problem. Decrease the time step.
Jiannan Wang (jnwang) said : | #2 |
Hello Robert,
Thank you again on these series questions.
I tried to reduced the time step by 1e-6 but still no luck: the temperature grown asymptomatically toward the target temperature and stayed there flat for a while. Then jumped cross the target temperature and grown at faster speed. The pore pressure didn't decrease as expected when the temperature reached the target temperature: only decrease a little bit then went up as the temperature crossed the target temperature.
If I decrease the time step more, the computing time would be astronomically increase. Do you have any tip? Or did I missed something when setting up the thermal engine? Thanks.
Best regards
Jiannan
Robert Caulk (rcaulk) said : | #3 |
Hey Jiannan,
>If I decrease the time step more, the computing time would be astronomically increase.
Thanks for testing the stability. The most sensitive scheme you are currently using is the fluidConduction scheme, especially if you have fluid cells whose "centers" are on top of one another (not uncommon in typical triangulations for distributed particle sizes). So consider what happens to the thermal resistivity between two fluid cells when the distance between them goes to 0. That is an important point to understand if you wish to understand the stability of the system.
The current way to "stabilize" these edge cases without decreasing the time step is to limit the allowed distance between cells, thus preventing the infinite thermal resistance between cells. Practically it is called thermal.
Cheers,
Robert
https:/
Jiannan Wang (jnwang) said : | #4 |
Hello Robert,
Thank yo for your suggestions and sorry for my delayed response. I tried to use thermal.
Best regards
Jiannan
Robert Caulk (rcaulk) said : | #5 |
Hello Jiannan,
I think one of your main issues is that you are not updating the triangulation often enough (meshUpdateInte
Another problem is that you are setting all the bodies to non-dynamic. That means you do not care about the mechanics. But I don't think that is the case, since you told me you are trying to simulate THM.
Another problem is that your packing is being stochastically generated - hard to really debug anything when the packing changes each run. You can use the seed argument of makeCloud to avoid this.
I guess there might be many other issues with your script. But you are requesting to run a fully coupled THM sim and so it is probably easier I just post another example script showing that, so I went ahead and did that [2].
Cheers,
Robert
[1]https:/
[2]https:/
Jiannan Wang (jnwang) said : | #6 |
Hello Robert,
Thank you so much for all the explanations and sharing the new example script. They are incredibly helpful. I'll work through them to have better understanding. This indeed solve the irritating puzzle I have been having.
Best regards
Jiannan Wang
Jiannan Wang (jnwang) said : | #7 |
Hello Robert,
Just a follow-up update, I ran the code as in [1], but the temperature still went beyond thermalBndCondV
Also, only changing the side of the boundaries (both pressure and thermal) also seems affect the result: such as changing from z-,z+ to x-, x+.
I'll try to run it with smaller time step and seem what may happen (tried to increase the density scale to 1e20 and 1e30 but no luck yet).
Best regards
Jiannan Wang
Robert Caulk (rcaulk) said : | #8 |
>I ran the code as in [1], but the temperature still went beyond thermalBndCondV
The code works for me.
Which version of yade are you using? Which installation method? At what time/timestep during the simulation does the temperature "go beyond the boundary value"? How far beyond the boundary value does it go? Does it happen in with multi and single core runs?
Details :-).
Robert Caulk (rcaulk) said : | #9 |
>also seems affect the result
Please, add details.
Jiannan Wang (jnwang) said : | #10 |
Hello Robert,
Thank you for being so patient. I used the yadedaily (20200829-
```
Yade version : 20200829-
Yade features : BoostLog PrecisionDouble Odeint VTK OpenMP GTS GUI-Qt5 CGAL PFVFLOW PFVFLOW LINSOLV MPI TWOPHASEFLOW FEMLIKE GL2PS LBMFLOW THERMAL PotentialParticles PotentialBlocks
Yade config dir: ~/.yadedaily
Yade precision : 53 bits, 15 decimal places, without mpmath
```
Libraries used :
| library | cmake | C++ |
| ------------- | -------
| boost | 106501 | 1.65.1 |
| cgal | | 4.11 |
| clp | 1.16.11 | 1.16.11 |
| cmake | 3.10.2 | |
| coinutils | 2.10.14 | 2.10.14 |
| compiler | /usr/bin/c++ 7.5.0 | gcc 7.5.0 |
| eigen | 3.3.4 | 3.3.4 |
| freeglut | 2.8.1 | |
| gl | | 20190911 |
| ipython | 5.5.0 | |
| metis | | 5.1.0 |
| mpi | 3.1 | ompi:2.1.1 |
| mpi4py | 2.0.0 | |
| openblas | | OpenBLAS 0.2.20 |
| python | 3.6.9 | 3.6.9 |
| qglviewer | | 2.6.3 |
| qt | | 5.9.5 |
| sphinx | 1.6.7-final-0 | |
| sqlite | | 3.22.0 |
| suitesparse | 5.1.2 | 5.1.2 |
| vtk | 6.3.0 | 6.3.0 |
Linux version: Ubuntu 18.04.5 LTS
Best regards
Jiannan Wang
Jiannan Wang (jnwang) said : | #11 |
Hello Robert,
I ran the simulation as it is, the temperature went beyond the boundary value (45) at 155 s and up to 47.8035. The body temperature and the pore pressure curve looks a little bit funny. Body and water temperature in the middle looks like [1] and the pore pressure in the middle looks like [2]. Hope it is OK to attach google link here. Thank you again.
Best regards
Jiannan
[1] https:/
[2 ] https:/
Robert Caulk (rcaulk) said : | #12 |
Hello,
Thanks for the information. Will you please try reducing the meshUpdateInterval to 2, and then reporting back?
thank you,
Robert
Jiannan Wang (jnwang) said : | #13 |
Hello Robert,
Absolutely. I'll get back to you as soon as I get the result.
Best regards
Jiannan Wang
Jiannan Wang (jnwang) said : | #14 |
Hello Robert,
I tried with meshUpdateInter
Best regards
Jiannan
Robert Caulk (rcaulk) said : | #15 |
Strange, for me it works with newest yadedaily. Not sure what to tell you! here is the exact script I am using:
#******
# Copyright (C) 2019 by Robert Caulk *
# <email address hidden> *
# *
# This program is free software; it is licensed under the terms of the *
# GNU General Public License v2 or later. See file LICENSE for details. *
#******
#
# Script demonstrating the use of full Thermo-
# in ThermalEngine. Methods published in:
#Robert Caulk, Luc Sholtès, Marek Krzaczek, Bruno Chareyre,
#A pore-scale thermo–
#Computer Methods in Applied Mechanics and Engineering,
#Volume 372,
#2020,
#113292,
#ISSN 0045-7825,
#https:/
#http://
from yade import pack, ymport, plot, utils, export, timing
import numpy as np
young=5e6
mn,mx=Vector3(
identifier = '-thm_coupling'
if not os.path.
os.mkdir(
if not os.path.
os.mkdir(
O.materials.
O.materials.
walls=aabbWalls
wallIds=
sp=pack.
sp.makeCloud(
sp.toSimulation
triax=TriaxialS
maxMultiplier=
finalMaxMultip
thickness = 0,
stressMask = 7,
internalCompac
)
O.engines=[
ForceResetter(),
InsertionSortC
InteractionLoop(
[Ig2_
[Ip2_
[Law2_
),
GlobalStiffnes
triax,
FlowEngine(
ThermalEngine(
VTKRecorder(
NewtonIntegrat
]
O.step()
ss2sc.interacti
is2aabb.
tri_pressure = 1000
triax.goal1=
triax.stressMask=7
while 1:
O.run(1000, True)
unb=unbalance
print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
if unb<0.1 and abs(-tri_
break
triax.internalC
flow.debug=False
# add flow
flow.permeabili
flow.pZero = 10
flow.meshUpdate
flow.fluidBulkM
flow.useSolver=4
flow.permeabili
flow.viscosity= 0.001
flow.decoupleForces = False
flow.bndCondIsP
flow.bndCondVal
## Thermal Stuff
flow.bndCondIsT
flow.thermalEng
flow.thermalBnd
flow.tZero=25
flow.dead=0
thermal.dead=1
thermal.
thermal.
thermal.debug=0
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.fluidK = 0.650
thermal.fluidBeta = 2e-5 # 0.0002
thermal.particleT0 = 25
thermal.particleK = 2.0
thermal.particleCp = 710
thermal.
thermal.
thermal.
thermal.
thermal.
timing.reset()
O.dt=1e-4
O.dynDt=False
thermal.dead=0
flow.emulateAct
def bodyByPos(x,y,z):
cBody = O.bodies[1]
cDist = Vector3(
for b in O.bodies:
if isinstance(b.shape, Sphere):
dist = b.state.pos - Vector3(x,y,z)
if np.linalg.
cDist = dist
cBody = b
return cBody
bodyOfInterest = bodyByPos(
from yade import plot
def history():
plot.addData(
ftemp1=
p=flow.
t=O.time,
i = O.iter,
bodyOfIntTemp = O.bodies[
)
#plot.
O.engines=
def pressureField():
flow.saveVtk(
O.engines=
plot.plots=
plot.plot(
O.run()
Jiannan Wang (jnwang) said : | #16 |
Hello Robert,
Thank you for all the helps and sharing. They are very helpful. I keep play around the script and see if I can make it work.
Best regards
Jiannan
Jiannan Wang (jnwang) said : | #17 |
Hello Robert,
I think I find out where I did wrong: I set ignoreFictiousC
Best regards
Jiannan
Robert Caulk (rcaulk) said : | #18 |
Happy to hear it :-)
ignoreFictiousC