Body temperature rise higher than thermalBndCondValue

Sounds like a stabiltiy problem. Decrease the time step.

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

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.minimumFluidCondDist [1].

Cheers,

Robert

https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.minimumFluidCondDist

Hello Robert,

Thank yo for your suggestions and sorry for my delayed response. I tried to use thermal.minimumFluidCondDist = minimum particle diameter but still no luck. The only difference is the temperature curve now went smoothly as "normal" across the pre-set boundary temperature. Maybe I did something wrong?

Best regards
Jiannan

Hello Jiannan,

I think one of your main issues is that you are not updating the triangulation often enough (meshUpdateInterval). But another is that your processes are on two incredibly different scales. Heat is slow. Elastic waves are fast. It is important to understand what regime you are in, which regime you want to be in, and how to get there. If you want to be in a quasistatic regime, then you can safely scale the density. This will let you take larger (stable) time steps to resolve the heat transfer. Please read the paper associated with ThermalEngine to become better acquainted with how these processes work and the methods driving them [1].

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://www.sciencedirect.com/science/article/pii/S0045782520304771?via%3Dihub
[2]https://gitlab.com/yade-dev/trunk/-/merge_requests/519

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

Hello Robert,

Just a follow-up update, I ran the code as in [1], but the temperature still went beyond thermalBndCondValue.

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

[1]https://gitlab.com/yade-dev/trunk/-/merge_requests/519

>I ran the code as in [1], but the temperature still went beyond thermalBndCondValue.

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 :-).

>also seems affect the result

Please, add details.

Hello Robert,

Thank you for being so patient. I used the yadedaily (20200829-4169~e627a44~bionic1) in the VritualBox on a Mac. I give 3 cores to the VirtualBox, but I didn't run it parallel. Here is the printAllVersions() output. I will re-run the exact setting again and keep you posted with more detail of how simulation ended.
```
Yade version : 20200829-4169~e627a44~bionic1
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

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://drive.google.com/file/d/1FiV8S1L-2J4mnV9l9-NWtw6zI2gvZQhM/view?usp=sharing
[2 ] https://drive.google.com/file/d/1YEaPS-lSdf5XuWVcfAAacmkF7Xoz6dMp/view?usp=sharing

Hello,

Thanks for the information. Will you please try reducing the meshUpdateInterval to 2, and then reporting back?

thank you,

Robert

Hello Robert,

Absolutely. I'll get back to you as soon as I get the result.

Best regards
Jiannan Wang

Hello Robert,

I tried with meshUpdateInterval=2, but the result is pretty much identical. Let me know.

Best regards
Jiannan

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-Hydro-Mechanical coupling
# in ThermalEngine. Methods published in:
#Robert Caulk, Luc Sholtès, Marek Krzaczek, Bruno Chareyre,
#A pore-scale thermo–hydro-mechanical model for particulate systems,
#Computer Methods in Applied Mechanics and Engineering,
#Volume 372,
#2020,
#113292,
#ISSN 0045-7825,
#https://doi.org/10.1016/j.cma.2020.113292.
#http://www.sciencedirect.com/science/article/pii/S0045782520304771

from yade import pack, ymport, plot, utils, export, timing
import numpy as np

young=5e6

mn,mx=Vector3(0,0,0),Vector3(0.05,0.05,0.05)

identifier = '-thm_coupling'

if not os.path.exists('VTK'+identifier):
 os.mkdir('VTK'+identifier)

if not os.path.exists('txt'+identifier):
 os.mkdir('txt'+identifier)

O.materials.append(FrictMat(young=young*100,poisson=0.5,frictionAngle=0,density=2600e10,label='walls'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(30),density=2600e10,label='spheres'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=0.0015,rRelFuzz=0.333,num=100,seed=11)
sp.toSimulation(color=(0.752, 0.752, 0.752),material='spheres')

triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = 0,
 stressMask = 7,
 internalCompaction=True,
)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1,label='is2aabb'),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1,label='ss2sc'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.5),
 triax,
 FlowEngine(dead=1,label="flow",multithread=False),
 ThermalEngine(dead=1,label='thermal'),
 VTKRecorder(iterPeriod=500,fileName='VTK'+identifier+'/spheres-',recorders=['spheres','thermal','intr'],dead=0,label='VTKrec'),
 NewtonIntegrator(damping=0.5)
]

O.step()
ss2sc.interactionDetectionFactor=-1
is2aabb.aabbEnlargeFactor=-1

tri_pressure = 1000
triax.goal1=triax.goal2=triax.goal3=-tri_pressure
triax.stressMask=7
while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
  if unb<0.1 and abs(-tri_pressure-triax.meanStress)/tri_pressure<0.001:
    break

triax.internalCompaction=False

flow.debug=False
# add flow
flow.permeabilityMap = False
flow.pZero = 10
flow.meshUpdateInterval=2
flow.fluidBulkModulus=2.2e9
flow.useSolver=4
flow.permeabilityFactor=-1e-5
flow.viscosity= 0.001
flow.decoupleForces = False
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondValue=[0,0,0,0,10,10]

## Thermal Stuff
flow.bndCondIsTemperature=[0,0,0,0,0,0]
flow.thermalEngine=True
flow.thermalBndCondValue=[0,0,0,0,0,0]
flow.tZero=25

flow.dead=0
thermal.dead=1

thermal.conduction=True
thermal.fluidConduction=True
thermal.debug=0
thermal.thermoMech=True
thermal.solidThermoMech = True
thermal.fluidThermoMech = True
thermal.advection=True
thermal.useKernMethod=False
thermal.bndCondIsTemperature=[0,0,0,0,0,1]
thermal.thermalBndCondValue=[0,0,0,0,0,45]
thermal.fluidK = 0.650
thermal.fluidBeta = 2e-5 # 0.0002
thermal.particleT0 = 25
thermal.particleK = 2.0
thermal.particleCp = 710
thermal.particleAlpha = 3.0e-5
thermal.particleDensity = 2700
thermal.tsSafetyFactor = 0 #0.01
thermal.uniformReynolds =10
thermal.minimumThermalCondDist=0

timing.reset()
O.dt=1e-4
O.dynDt=False
thermal.dead=0
flow.emulateAction()

def bodyByPos(x,y,z):
 cBody = O.bodies[1]
 cDist = Vector3(100,100,100)
 for b in O.bodies:
  if isinstance(b.shape, Sphere):
   dist = b.state.pos - Vector3(x,y,z)
   if np.linalg.norm(dist) < np.linalg.norm(cDist):
    cDist = dist
    cBody = b
 return cBody

bodyOfInterest = bodyByPos(0.025,0.025,0.025)

from yade import plot

def history():
 plot.addData(
  ftemp1=flow.getPoreTemperature((0.025,0.025,0.025)),
  p=flow.getPorePressure((0.025,0.025,0.025)),
  t=O.time,
  i = O.iter,
  bodyOfIntTemp = O.bodies[bodyOfInterest.id].state.temp,
)
 #plot.saveDataTxt('Record.txt',vars=('t','i','p','ftemp1','bodyOfIntTemp'))

O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]

def pressureField():
 flow.saveVtk('VTK'+identifier+'/',withBoundaries=False)

O.engines=O.engines+[PyRunner(iterPeriod=2000,command='pressureField()')]

plot.plots={'t':(('ftemp1','k-'),('bodyOfIntTemp','r-'))}

plot.plot(subPlots=False)
O.run()

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

Hello Robert,

I think I find out where I did wrong: I set ignoreFictiousConduction=True. Though I'm not clear why exactly ... But now the script seems working fine. Thank yo again for all the helps.

Best regards
Jiannan

Happy to hear it :-)

ignoreFictiousConduction is basically just saying that the spheres that touch the "fictious" cells (which are the cells that form the boundary) will not interact conductively with the boundary fluid. It was added for debugging purposes - but in some cases it could be useful such as regular packings and analytical comparisons where fictious cells may have weird geometrical representations that we want to ignore. Best to leave it at its default value for full THM models.