yade-users team mailing list archive

[Túlio Mumić Cunha <question688652@xxxxxxxxxxxxxxxxxxxxx>]

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.