Question #703169 on Yade changed:
https://answers.launchpad.net/yade/+question/703169

Xu gave more information on the question:
Hi Karol,

Sorry, the format of the code changes when i copy and paste it. I have
already corrected it shown as below.  There are two scripts, first is
the consolidation and the second is the undrained test. After you run
through the first script, a file of consolidated sample will be get,
which can be directly loaded by the second script.

> - minimal = short. It is difficult to find mistakes in long pieces of code 
> and it is usually much more annoying.
Sorry for that,  I have done my best to shorten the code. 

This is the corrected script
######Consolidation#######
from yade import pack,plot,qt,export,ymport,utils
import numpy as np
IsoSigma = -1.e5
O.periodic=True
O.trackEnergy=True
idSand=O.materials.append(FrictMat(young=1e8,poisson=.8,frictionAngle=0.2,density=26500000,label='spheres'))
sp1=pack.SpherePack()
sp1.makeCloud(maxCorner=(0.1, 0.1, 0.1), psdSizes=[0.0017, 0.00191, 0.002285, 
0.0026, 0.00292, 0.00325, 0.0035], psdCumm=[0.0, 0.1, 0.3, 0.5, 0.6, 0.9, 1], 
periodic=True,num=2000,seed=1)
sp1.toSimulation(color=(0,0,1),material=idSand)
v=utils.getSpheresVolume()
print utils.getSpheresVolume()
fout = file('Packing_volume.txt','w')
fout.write(str(v))
fout.close()

relRadList2 = [1,1]
relPosList2 = [[0.4,0,0],[-0.4,0,0]]

templates= []
templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))
O.bodies.replaceByClumps(templates,[1.],discretization=10)

law=Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True)
O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb()]),
        InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [law]),
        GlobalStiffnessTimeStepper(),
        NewtonIntegrator(damping=0.20),
        PeriTriaxController(    goal=(IsoSigma,IsoSigma,IsoSigma), # Vector6 of 
prescribed final values
                                stressMask=7,
                                dynCell=True,
                                maxStrainRate=(5.e+0,5.e+0,5.e+0),
                                maxUnbalanced=0.001,
                                relStressTol=1.e-3,
                                doneHook='Finished()',
                                label='p3d'
        ),
]

def Finished():
        O.save('post_iso.xml')
        print 'Finished'
        print (O.cell.size[0]*O.cell.size[1]*O.cell.size[2]-v)/v
        print getStress()
        print unbalancedForce()
        O.pause()

O.run()

########Undrained triaxial test#########
from yade import pack,plot,qt,export
from math import fabs
import numpy as np
O.load('post_iso.xml')
O.trackEnergy=True
setContactFriction(0.5)
filename1='Packing_volume.txt'
a=np.loadtxt(filename1)
Ep_ini=O.energy['elastPotential']
Edamp_ini=O.energy['nonviscDamp']
Ekin_ini=utils.kineticEnergy()
work=0

law=Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True)
O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb()]),
        InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [law]),
        NewtonIntegrator(damping=0.2),
        GlobalStiffnessTimeStepper(),
        PyRunner(command='inputwork()',iterPeriod=1),
        PyRunner(command='plotAddData()',iterPeriod=10),
]

def inputwork():
        global work
        
dwork=(0.5*(O.cell.prevVelGrad+O.cell.velGrad)*utils.getStress()).trace()*O.dt*(O.cell.hSize).determinant()
        work=work+dwork

def plotAddData():
        global Ep_ini
        global Edamp_ini
        global Ekin_ini
        global work
        plot.addData(
                iter=O.iter,iter_=O.iter,
                sxx=utils.getStress()[0,0],
                syy=utils.getStress()[1,1],
                szz=utils.getStress()[2,2],
                exx=O.cell.getSmallStrain()[0,0],
                eyy=O.cell.getSmallStrain()[1,1],
                ezz=O.cell.getSmallStrain()[2,2],
                Z=avgNumInteractions(),
                Zm=avgNumInteractions(skipFree=True),
                voidratio=(O.cell.size[0]*O.cell.size[1]*O.cell.size[2]-a)/a,
                unbalanced=utils.unbalancedForce(),
                t=O.time,
                inwork=work,
                gWork=O.energy['gravWork'],
#               Ep=O.energy['elastPotential'],
                Ep=law.elasticEnergy()-Ep_ini,
                Edamp=O.energy['nonviscDamp']-Edamp_ini,
                Ediss=law.plasticDissipation(),
#               Ediss=O.energy['plastDissip'],
                Ekin=utils.kineticEnergy(),
                
Etot=law.elasticEnergy()-Ep_ini+O.energy['nonviscDamp']-Edamp_ini+law.plasticDissipation()+utils.kineticEnergy()-Ekin_ini,
#               Etot=O.energy.total(),**O.energy

        )
        
plot.saveDataTxt('macroFile',vars=('t','exx','eyy','ezz','sxx','syy','szz','Z','Zm','voidratio'))
        
plot.saveDataTxt('energyFile',vars=('t','unbalanced','gWork','Edamp','Ekin','Ep','Ediss','Etot','inwork'))

O.cell.trsf=Matrix3.Identity
O.cell.velGrad=Matrix3(-0.1,0,0,0,0.05,0,0,0,0.05)
O.run()


###########################################
My problom is the difference of "inwork" and "Etot" in the second script when 
using Clump particles.
I really appreciate if you can confirm the issue or help me to solve the 
problem.
Thanks a lot,

Xu

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.

_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to     : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Reply via email to