yade-users team mailing list archive

[Eleonora Piersanti <eleonorapiersanti@xxxxxxxx>]

Hi, 
I am trying to run a simulation using a predicate filled with spheres with a hexagonal packing plus more spheres added with a for loop around it.I added wall to the simulation as well so that I could impose boundary conditions for the FlowEngine.From what I have seen the triangulation is created only after at least one time step run. And to run the simulation it isneeded to impose boundary conditions that can be done only if walls are added (since there are no cells and it is not possible 
to impose a certain pressure or flux on some cells).   
Anyway, I still obtain the error Infinite K1!Do you have any idea why this can happen? I am using a Cohesive type material.
Also, the triangulation that I obtain for the fluid is not the one I desire. Some points are connected even if I would like them to be not.Is there any way to control this?
I attach the surface and the script.
Thanks for your help.

Attachment: 2DBrain.gts
Description: Binary data

from yade import pack
import gts, os.path, locale
from yade import qt
import csv
import numpy as np

locale.setlocale(locale.LC_ALL, 'en_US.UTF-8')   #gts is locale-dependend.  If, for example, german locale is used, gts.read()-function does not import floats normally

surf=gts.read(open('2DBrain.gts'))

Rin = 0.03
Rout = 0.1
rp = 0.0014000


idTissue=O.materials.append(CohFrictMat(young=500.0, poisson=.35, density=1000.0, fragile=False, frictionAngle=.6,label="concrete"))
wallMaterial=O.materials.append(CohFrictMat(young=50000.0, poisson=.35, density=1000.0, fragile=False, frictionAngle=.6,label="walls"))

minX = -0.10671062880416388
maxX = 0.10671062880416388
minY = -0.2363490390055643
maxY = 0.10670230537773978
minZ = 0.0
maxZ = 0.0059723808531952655
mn,mx=Vector3(minX,minY,minZ),Vector3(maxX,maxY,maxZ) # corners of the initial packing
walls=aabbWalls([mn,mx],thickness=0, material="walls")
wallIds=O.bodies.append(walls)

pred=pack.inGtsSurface(surf)
spheres=pack.regularHexa(pred, radius=rp, gap=0.0, color=(1,0,0), material=idTissue)
O.bodies.append(spheres)

radius = np.arange(Rin+rp, Rout-rp, 2.0*rp)
rSAS = 0.1 + 2.0 * rp 
rSkull = rSAS + 2.0*rp 
alphaSAS = 4.0*np.arcsin(rp/(2*rSAS)) 
alphaSkull = 4.0*np.arcsin(rp/(2*rSkull))
toll = 0.0
ntheta = 250
theta = np.linspace(0, 2.0*np.pi, ntheta, endpoint=False)
dt = theta[1] - theta[0]
tlimit = []

rb = []
for r in [rSAS, rSkull]:
    a = r*np.sin((theta[1] - theta[0])/2.0)
    b = r - r*np.cos((theta[1] - theta[0])/2.0)
    rb += [np.sqrt(a*a + b*b)]

print rb
for z in [1, 3]:
    for t in theta:
        if (t < 3.0/2.0 *np.pi - 3*dt or t > 3.0/2.0 *np.pi + 3*dt):
            s = utils.sphere(center=[rSAS*np.cos(t), rSAS*np.sin(t), rb[0]*z], radius=rb[0], material=idTissue)
            O.bodies.append(s)
            s = utils.sphere(center=[(rSAS+rb[0]+rb[1])*np.cos(t), (rSAS+rb[0]+rb[1])*np.sin(t), rb[0]*z], radius=rb[0], material=idTissue)
            O.bodies.append(s)
        if (t >= 3.0/2.0 * np.pi - 5*dt and t <= 3.0/2.0 * np.pi - 3*dt) or (t >= 3.0/2.0 * np.pi + 3*dt and t <= 3.0/2.0 * np.pi + 5*dt):
            tlimit += [t]
            for i in range(1,50):
                s = utils.sphere(center=[(rSAS+rb[0]+rb[1])*np.cos(t), (rSAS+rb[0]+rb[1])*np.sin(t)-2*i*rb[1], rb[0]*z], radius=rb[0], material=idTissue)
                O.bodies.append(s)
          

qt.View()

newton=NewtonIntegrator(damping=0.5)

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
    [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
    [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
),
FlowEngine(label="flow"),#introduced as a dead engine for the moment, see 2nd section
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
# triax,
newton,
VTKRecorder(fileName='./VTK/',recorders=['all'],iterPeriod=1)
]

#setContactFriction(radians(finalFricDegree))

#B. Activate flow engine and set boundary conditions in order to get permeability
# flow.dead=0
print "defining flow parameters"
flow.debug=True
O.dt=0.1e-6
flow.dead=0
flow.defTolerance=0.3
flow.meshUpdateInterval=200
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=10

flow.bndCondIsPressure=[1,1,1,1,1,1]
flow.bndCondValue=[0,0,0,0,0,0]

O.run(1, True)
flow.saveVtk("./VTK")