Problem with using blockCell command in flow/thermal engine
Hello Robert,
In order to increase timestep in advection simulation, I would like to block cells with volume less than a min volume.
I tried blockCell command [1] in your flowScenario.py example and blocked cells located on the top (z>0.04); here is what happens:
1- Temperature in the blocked cells remained as initiated; that makes sense, right? but the problem is that the heat transfer between hot flow and solid stops working. I mean, the solid does not get heat up and remain at initial temperature
2- Pressure gradient in blocked and non-blocked cells are the same. So, it seems like flow calculations still occur in all cells.
Please let me know if I am not applying blockCell in the right way.
Thanks
Zoheir
Here is what I added to the flowScenario.py example:
# Blocking cells in the upper part of the model
#++++++
flow.emulateAct
flow.updateTria
for i in range(0,
coords = flow.getCellCen
if coords[2]>.04:
flow.
#++++++
Here is the entire MWE:
# -*- encoding=utf-8 -*-
#******
# 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 ThermalEngine by permeating warm fluid
# through a cold packing. Also serves as a validation script for comparison
# with ANSYS CFD. See details in
# Caulk, R., Scholtes, L., Kraczek, M., Chareyre, B. (In Print) A
# pore-scale Thermo-
# Computer Methods in Applied Mechanics and Engineering. Accepted July 2020.
#
# note: warnings for inifiniteK and Reynolds numbers = nan for boundary
# cells in regular packings are expected. It does not interfere with the
# physics
from yade import pack, ymport
from yade import timing
import numpy as np
import shutil
timeStr = time.strftime(
num_spheres=1000# number of spheres
young=1e9
rad=0.003
mn,mx=Vector3(
thermalCond = 2. #W/(mK)
heatCap = 710. #J(kg K)
t0 = 333.15 #K
# micro properties
r = rad
k = 2.0 # 2*k*r
Cp = 710.
rho = 2600.
D = 2.*r
m = 4./3.*np.
# macro diffusivity
identifier = '-flowScenario'
if not os.path.
os.mkdir(
else:
shutil.
os.mkdir(
if not os.path.
os.mkdir(
else:
shutil.
os.mkdir(
shutil.
O.materials.
O.materials.
walls=aabbWalls
wallIds=
sp = O.bodies.
print('num bodies ', len(O.bodies))
triax=TriaxialS
maxMultiplier=
finalMaxMultip
thickness = 0,
stressMask = 7,
internalCompac
)
ThermalEngine = ThermalEngine(
newton=
intRadius = 1
O.engines=[
ForceResetter(),
InsertionSortC
InteractionLoop(
[Ig2_
[Ip2_
[Law2_
),
FlowEngine(
ThermalEngine, GlobalStiffness
triax,
VTKRecorder(
newton
]
#goal = -1e5
#triax.
for b in O.bodies:
if isinstance(b.shape, Sphere):
b.dynamic=False # mechanically static
flow.dead=0
flow.defToleran
flow.meshUpdate
flow.useSolver=4
flow.permeabili
flow.viscosity= 0.001
flow.bndCondIsP
flow.bndCondVal
flow.thermalEng
flow.debug=False
flow.fluidRho = 997
flow.fluidCp = 4181.7
flow.getCHOLMOD
flow.bndCondIsT
flow.thermalEng
flow.thermalBnd
flow.tZero=t0
flow.pZero=0
flow.maxKdivKmean=1
flow.minKdivmea
thermal.dead=0
thermal.debug=False
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.fluidK = 0.6069 #0.650
thermal.
thermal.particleT0 = t0
thermal.
thermal.particleK = thermalCond
thermal.particleCp = heatCap
thermal.
#thermal.
timing.reset()
O.dt=0.1e-3
O.dynDt=False
# Blocking cells in the upper part of the model
#++++++
flow.emulateAct
flow.updateTria
for i in range(0,
coords = flow.getCellCen
if coords[2]>.04:
flow.
#++++++
O.run(1,1)
flow.dead=0
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
print('found closest body ', cBody.id, ' at ', cBody.state.pos)
return cBody
#bodyOfInterest = bodyByPos(
bodyOfInterest = bodyByPos(
# find 10 bodies along x axis
axis = np.linspace(mn[0], mx[0], num=5)
axisBodies = [None] * len(axis)
axisTrue = np.zeros(len(axis))
for i,x in enumerate(axis):
axisBodies[i] = bodyByPos(x, mx[1]/2, mx[2]/2)
axisTrue[i] = axisBodies[
print("found body of interest at", bodyOfInterest.
from yade import plot
## a function saving variables
def history():
plot.addData(
ftemp1=
p=flow.
t=O.time,
i = O.iter,
temp1 = axisBodies[
temp2 = axisBodies[
temp3 = axisBodies[
temp4 = axisBodies[
temp5 = axisBodies[
bodyOfIntTemp = O.bodies[
)
plot.saveDataT
O.engines=
def pressureField():
flow.saveVtk(
O.engines=
def endFlux():
if O.time >= 30:
flux = 0
n=utils.
for i in flow.getBoundar
flux +=i[0]*i[3]/n # area * velocity / porosity (dividing by porosity because flow engine is computing the darcy velocity)
massFlux = flux * 997
K = abs(flow.
d=8e-3 # sphere diameter
Kc = d**2/180. * (n**3.)/(1.-n)**2
print(
print('outlet flux(with vels only):'
print('sim paused')
O.pause()
O.engines=
VTKrec.dead=0
from yade import plot
plot.plots=
plot.plot()
O.saveTmp()
print("starting thermal sim")
O.run()
[1] https:/
Question information
- Language:
- English Edit question
- Status:
- Expired
- For:
- Yade Edit question
- Assignee:
- No assignee Edit question
- Last query:
- Last reply: