\( \usepackage{euler} % must be loaded before fontspec for the whole doc (below); this must be kept for pngmath, however \usepackage{hyperref} \usepackage{amsmath} \usepackage{amsbsy} %\usepackage{mathabx} \usepackage{underscore} \usepackage[all]{xy} % Metadata of the pdf output \hypersetup{pdftitle={Yade Documentation}} \hypersetup{pdfauthor={Smilauer, V., et al.}} \hypersetup{pdfkeywords={Discrete element method, granular materials, contact law, hydromechanical coupling}} % symbols \renewcommand{\mat}{\boldsymbol} \renewcommand{\vec}{\boldsymbol} \renewcommand{\tens}{\boldsymbol} %normalized and locfram here \def\normalized#1{\widehat{#1}} \def\locframe#1{\widetilde{#1}} % timestep \def\Dt{\Delta t} \def\Dtcr{\Dt_{\rm cr}} % algorithm complexity \def\bigO#1{{\mathcal{O}(#1)}} % variants for greek symbols \let\epsilon\varepsilon \let\theta\vartheta \let\phi\varphi % shorthands \let\sig\sigma \let\eps\epsilon % variables at different points of time \def\prev#1{#1^-} \def\pprev#1{#1^\ominus} \def\curr#1{#1^{\circ}} \def\nnext#1{#1^\oplus} \def\next#1{#1^+} % shorthands for geometry \def\currn{\curr{\vec{n}}} \def\currC{\curr{\vec{C}}} \def\uT{\vec{u}_T} \def\curruT{\curr{\vec{u}}_T} \def\prevuT{\prev{\vec{u}}_T} \def\currn{\curr{\vec{n}}} \def\prevn{\prev{\vec{n}}} % motion \def\pprevvel{\pprev{\dot{\vec{u}}}} \def\nnextvel{\nnext{\dot{\vec{u}}}} \def\curraccel{\curr{\ddot{\vec{u}}}} \def\prevpos{\prev{\vec{u}}} \def\currpos{\curr{\vec{u}}} \def\nextpos{\next{\vec{u}}} \def\curraaccel{\curr{\dot{\vec{\omega}}}} \def\pprevangvel{\pprev{\vec{\omega}}} \def\nnextangvel{\nnext{\vec{\omega}}} \def\loccurr#1{\curr{\locframe{#1}}} \def\numCPU{n_{\rm cpu}} \DeclareMathOperator{\Align}{Align} \DeclareMathOperator{\sign}{sgn} % sorting algorithms \def\isleq#1{\currelem{#1}\ar@/^/[ll]^{\leq}} \def\isnleq#1{\currelem{#1}\ar@/^/[ll]^{\not\leq}} \def\currelem#1{\fbox{$#1$}} \def\sortSep{||} \def\sortInv{\hbox{\phantom{||}}} \def\sortlines#1{\xymatrix@=3pt{#1}} \def\crossBound{||\mkern-18mu<} \)

Part 2 of our Thermal Hands-on session will focus on the full THM coupling.

# Day 3 - Thermal Hands-on part 2

Let’s build a triaxially loaded cubic specimen:

    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)

    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=1)
    sp.toSimulation(color=(0.752, 0.752, 0.752),material='spheres')

Here we see that we are appending a sphere cloud to the simulation (we will compact them after setting the O.engines list).

Next, we need to construct our engines list as usual:

    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),
        TriaxialStressController(label='triax'),
        FlowEngine(dead=1,label="flow"),
        ThermalEngine(dead=1,label='thermal'),
        VTKRecorder(iterPeriod=500,fileName='./spheres-',recorders=['spheres','thermal','intr'],dead=1,label='VTKrec'),
        NewtonIntegrator(damping=0.5)
    ]

Now we have the full O.engines list set, which contains a TriaxialStressController() for our stress control, a FlowEngine() for the fluid fluxes and heat advection, and a ThermalEngine() for our thermal coupling.

## Compacting the specimen

Let’s setup the TriaxialStressController() for our compaction:

    triax.maxMultiplier=1.+2e4/young
    triax.finalMaxMultiplier=1.+2e3/young
    triax.thickness = 0
    triax.stressMask = 7
    triax.internalCompaction=True
    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

Here we see that we are running a loop where we run 1000 iterations of internalCompaction (particles grow in radius to achieve stress), then testing the unbalancedForce() and ultimately stopping if our stopping criteria is achieved. We can’t forget that our FlowEngine() and ThermalEngine() are both set to dead=1 in the O.engines list, so they will not activate during this compaction stage.

## Setting up the FlowEngine()

    # initial pressure condition
    flow.pZero = 10
    flow.meshUpdateInterval = 2
    # we will activate compressibility in the fluid
    flow.fluidBulkModulus = 2.2e9
    flow.useSolver = 4
    # enforcing a darcy permeability in the specimen
    flow.permeabilityFactor = -1e-5
    flow.viscosity = 0.001
    # setting the boundary conditions
    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]
    # activate the thermal engine
    flow.thermalEngine = True
    flow.thermalBndCondValue = [0,0,0,0,0,0]
    # initial temperature conditions
    flow.tZero = 25

    flow.dead=0

## Setting up the ThermalEngine()

    thermal.dead = 0
    thermal.conduction = True
    thermal.fluidConduction = True
    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
    thermal.particleT0 = 25
    thermal.particleK = 2.0
    thermal.particleCp = 710
    thermal.particleAlpha = 3.0e-5
    thermal.particleDensity = 2700
    thermal.tsSafetyFactor = 0
    thermal.uniformReynolds = 10

We won’t describe each parameter here, those descriptions can be found in the Class Reference. However, it is clear we are activating conduction, advection, the thermo-fluid mechanical coupling, the solid-fluid mechanical coupling, and fluid conduction. Each component can be deactivated in case the user does not need the full THM coupling. We also see a similar assignment of boundary conditions as we saw in the previous hands-on sessions. Some additional parameters shown here include the fluid thermal conductivity (thermal.fluidK), the coefficient of thermal expansion for the fluid (thermal.fluidBeta).

## Running the coupled simulation

The simulation is set and ready to run, first we will let FlowEngine() detect and assign the boundary conditions by running flow.emulateAction():

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

Now it is up to you to finish the script

1. collect the temperature at some interesting points in the specimen

2. plot the temperature

3. export the VTK files for viewing in paraview

## Example script

Please find a full script located in the examples folder