New question #688652 on Yade:
https://answers.launchpad.net/yade/+question/688652
I am currently trying to make the elastic calibration of the CpmMat model but
the values of the computed poisson, based on the uniax.py code, are not
changing even though I change in the CpmMat the values of 'young' and
'poisson'. Below is my simulation code on yade 2019.1a:
###################################################
#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from yade import plot,pack,timing
import time, sys, os, copy, numpy as np
savedir = './elastmat/'
# default parameters or from table
readParamsFromTable(
young = 24e9,
poisson = 0.8,
#sigmaT=18e6,
#frictionAngle=atan(0.8),
epsCrackOnset=1e-5,
relDuctility=0.1,
intRadius=1.5,
dtSafety=.4,
damping=0.4,
strainRateTension=.05,
strainRateCompression=.5,
setSpeeds=True,
# 1=tension, 2=compression (ANDed; 3=both)
doModes=1,
specimenLength=.039850,
specimenRadius=.019875,
sphereRadius=0.75e-3,
# isotropic confinement (should be negative)
isoPrestress=0,
noTableOk=True,
# Number of elements to fetch in order to compute Poisson
n=500
)
from yade.params.table import *
if 'description' in list(O.tags.keys()):
O.tags['id']=O.tags['id']+O.tags['description']
print(young, poisson)
concreteId=O.materials.append(
CpmMat(
young=young,
#frictionAngle=frictionAngle,
poisson=poisson,
density=6420,
#sigmaT=sigmaT,
relDuctility=relDuctility,
epsCrackOnset=epsCrackOnset,
#isoPrestress=isoPrestress
)
)
sps=SpherePack()
#sp=pack.randomDensePack(pack.inAlignedBox((0,0,0),(specimenLength,specimenLength,specimenLength)),
sp=pack.randomDensePack(pack.inCylinder((0,0,-.5*specimenLength),(0,0,.5*specimenLength),specimenRadius),
spheresInCell=5000,
radius=sphereRadius,
memoizeDb='/tmp/calibCylinder3.sqlite',
returnSpherePack=True)
sp.toSimulation(material=concreteId)
def getElementsCurrentPosition():
return np.array(list(map(lambda x: [x.state.pos[0], x.state.pos[1],
x.state.pos[2], x.id], O.bodies)))
ibpos = getElementsCurrentPosition()
def fetchElements(pos=(0,0,0), n=1):
pos = np.array(pos)
distances = np.sum(ibpos[:,:3]-pos, axis=1)**2
return ibpos[distances.argsort()][:n,3].tolist()
def computePosition(n=1):
elements = np.array(fetchElements(n=n)).astype(int)
current_position = getElementsCurrentPosition()
return current_position[elements]
allpos = computePosition(n)
bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
O.dt=dtSafety*PWaveTimeStep()
print('Timestep',O.dt)
mm,mx=[pt[axis] for pt in aabbExtrema()]
imm, imx = mm, mx
coord_25,coord_50,coord_75=mm+.25*(mx-mm),mm+.5*(mx-mm),mm+.75*(mx-mm)
area_25,area_50,area_75=approxSectionArea(coord_25,axis),approxSectionArea(coord_50,axis),approxSectionArea(coord_75,axis)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),],verletDist=.05*sphereRadius),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc')],
[Ip2_CpmMat_CpmMat_CpmPhys()],
[Law2_ScGeom_CpmPhys_Cpm()],
),
NewtonIntegrator(damping=damping,label='damper'),
CpmStateUpdater(realPeriod=.5),
UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
PyRunner(virtPeriod=1e-6/strainRateTension,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True),
PyRunner(realPeriod=1,command='stopIfDamaged()',label='damageChecker'),
]
# plot stresses in ¼, ½ and ¾ if desired as well; too crowded in the graph that
includes confinement, though
plot.plots={'eps':('sigma')} #,'sigma.50')},'t':('eps')}
#'sigma.25','sigma.50','sigma.75')}
O.saveTmp('initial');
O.timingEnabled=False
global mode
mode='tension' if doModes & 1 else 'compression'
def initTest():
global mode
print("init")
if O.iter>0:
O.wait();
O.loadTmp('initial')
print("Reversing plot data"); plot.reverseData()
else: plot.plot(subPlots=False)
strainer.strainRate=abs(strainRateTension) if mode=='tension' else
-abs(strainRateCompression)
try:
from yade import qt
renderer=qt.Renderer()
renderer.dispScale=(1000,1000,1000) if mode=='tension' else
(100,100,100)
except ImportError: pass
print("init done, will now run.")
O.step(); # to create initial contacts
# now reset the interaction radius and go ahead
ss2sc.interactionDetectionFactor=1.
is2aabb.aabbEnlargeFactor=1.
O.run()
def stopIfDamaged():
global mode, allpos
if O.iter<2 or 'sigma' not in plot.data: return # do nothing at the
very beginning
sigma,eps=plot.data['sigma'],plot.data['eps']
allpos = np.append(allpos, computePosition(n), 1)
print(allpos.shape)
extremum=max(sigma) if (strainer.strainRate>0) else min(sigma)
minMaxRatio=0.8 if mode=='tension' else 0.8
if extremum==0: return
import sys; sys.stdout.flush()
if abs(sigma[-1]/extremum)<minMaxRatio or abs(strainer.strain)>(5e-3 if
isoPrestress==0 else 5e-2):
if mode=='tension' and doModes & 2: # only if compression is
enabled
numpy.savetxt(savedir+str(young)+' '+str(poisson)+'
1.txt', numpy.array([sigma, eps]).T, header='sigma\teps', delimiter='\t')
numpy.savetxt(savedir+str(young)+' '+str(poisson)+'
3.txt', allpos, delimiter='\t')
mode='compression'
O.save('/tmp/uniax-tension.yade.gz')
print("Saved /tmp/uniax-tension.yade.gz (for use with
interaction-histogram.py and uniax-post.py)")
print("Damaged, switching to compression... ");
O.pause()
# important! initTest must be launched in a separate
thread;
# otherwise O.load would wait for the iteration to
finish,
# but it would wait for initTest to return and deadlock
would result
import _thread; _thread.start_new_thread(initTest,())
return
else:
numpy.savetxt(savedir+str(young)+' '+str(poisson)+'
2.txt', numpy.array([sigma, eps]).T, header='sigma\teps', delimiter='\t')
numpy.savetxt(savedir+str(young)+' '+str(poisson)+'
4.txt', allpos, delimiter='\t')
print("Damaged, stopping.")
ft,fc=max(sigma),min(sigma)
if doModes==3:
print('Strengths fc=%g, ft=%g,
|fc/ft|=%g'%(fc,ft,abs(fc/ft)))
if doModes==2:
print('Compressive strength fc=%g'%(abs(fc)))
if doModes==1:
print('Tensile strength ft=%g'%(abs(ft)))
title=O.tags['description'] if 'description' in
list(O.tags.keys()) else O.tags['params']
print('Bye.')
O.pause()
def addPlotData():
yade.plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress+isoPrestress,})
initTest()
waitIfBatch()
###################################################
If I run this code, I'll have as result two files, one for computing the Young
modulus (which is trivial) and another for computing the Poisson modulus (file
contains the 500 most central elements positions in which I compute the eps on
the three axis to then compute the poisson modulus, whose algorythm is based in
the process described on Šmilauer's thesis, page 53
(https://dspace.cvut.cz/bitstream/handle/10467/79056/F1-D-2018-Smilauer-Vaclav-thesis.pdf?sequence=-1&isAllowed=y).
The code for Poisson calculus can be found below, in which cutoff is a
parameter I manually change in order to get only the linear part of the curve:
###################################################
def computePoisson(PATH, cutoff=20):
from scipy.stats import linregress
import numpy as np
pos = np.loadtxt(PATH)
x = pos[:,::4]
y = pos[:,1::4]
z = pos[:,2::4]
dx = x-x[:,0,None]
dy = y-y[:,0,None]
dz = z-z[:,0,None]
poissons = []
for x_, dx_, y_, dy_, z_, dz_ in zip(x.T, dx.T, y.T, dy.T, z.T, dz.T):
ex = linregress(x_, dx_)[0]
ey = linregress(y_, dy_)[0]
ez = linregress(z_, dz_)[0]
v = (-0.5 * (ex + ey)) / ez
poissons.append(v)
return(np.mean(poissons[cutoff:-cutoff]))
###################################################
By doing so, I can evaluate the variation of Macro Young based on Micro Young
and Micro Poisson, but on this case my Micro Poisson will be 0.187 regardless
of the variation of these parameters.
If anyone could help on this aspect, thank you, I have no idea on what is
causing this non-variability.
I tried searching other topics like
https://answers.launchpad.net/yade/+question/685847 and
https://answers.launchpad.net/yade/+question/315620
Is there any key point I'm obviously missing in my analysis?
Thank you beforehand and best regards!
--
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 : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
