On 17 April 2010 16:39, Nasibeh Moradi <[email protected]> wrote:
> hi, > I have written two contact law(simple and cohesive) and lp functor, > I want to check them, how can I check it? > To test a contact law I would suggest you to take a very simple case, like a sphere-sphere interaction (or if you prefer a box-sphere interaction) as a case of study. What (generally) you check is that the numerical solution and the analytical solution you compare are the same. The numerical solution is that one you obtain by the code, like the normal force at contact with respect to the relative approach. The analytical solution is simply the law you implemented (still in this case the relationship between forces and displacements at contacts, both in the normal and tangential direction). According to the code, you can proceed both in displacement or force control. If you also include damping in the force law, you should check, for instance, the relative displacement over the time, so that your analytical solution is now the solution to the differential equation that represents the motion (easy to get it by hand). Read also this paper, it might be useful http://cedb.asce.org/cgi/WWWdisplay.cgi?160824 cheers, Chiara > I used UniaxialStrainer and drew the graph sigma:epsilon, and obtained > E(young module) from slip of graph, > but the result is bad. > My material have E=74e9, from graph obtained ~30e9 (from my law). > I did uniaxial test with your law, the graph is attached. the graph is > same for both of compression and tension . > I forget anything? > The crossSection from utils.uniaxialTestFeatures() is 2*Area(Area is > calculated by hand). it is possible?! > > > The code: > > import itertools, time > from yade import utils, pack,timing,plot > > planeSize=Vector3(0.08,0.08,0.006) > rSphere=0.001 > strainRate=0.001 > > #Add material > > O.materials.append(FrictMat(young=74.1e9,frictionAngle=0.179,poisson=0.2,density=2500, > label='PlaneMaterial')) > > #generate a FCC Pack > a=2*rSphere; > Spheres=[] > h=a*sqrt(2)/2 > ii,jj,kk=[range(0,10),range(0,20),range(0,30)] > > for i,j,k in itertools.product(ii,jj,kk): > x,y,z=0+rSphere+i*sqrt(2)*a,0+rSphere+j*h,0+rSphere+k*h > if j%2==0: x+= sqrt(2)*a/2. if k%2==0 else -sqrt(2)*a/2. > if k%2!=0: x+=sqrt(2)*a/2. > > Spheres+=[utils.sphere((x,y,z),radius=rSphere,color=(1,0,1),material='PlaneMaterial')] > > PlaneIds=O.bodies.append(Spheres) > > bb=utils.uniaxialTestFeatures() > > negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area'] > O.dt=0.8*utils.PWaveTimeStep() > > #O.initializers=[ > # BoundDispatcher([Bo1_Sphere_Aabb()]) > #] > O.engines=[ > ForceResetter(), > > BoundDispatcher([ > Bo1_Sphere_Aabb(aabbEnlargeFactor=1.07) > ]), > > InsertionSortCollider(), > > InteractionDispatchers( > [Ig2_Sphere_Sphere_Dem3DofGeom()], > [Ip2_FrictMat_FrictMat_FrictPhys()], > [Law2_Dem3DofGeom_FrictPhys_Basic()] > ), > > #GravityEngine(gravity=[0,0,-9.81]), > NewtonIntegrator(), > UniaxialStrainer(strainRate=-strainRate,axis=axis,asymmetry=0, > > posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea, > blockDisplacements=False,blockRotations=False, > setSpeeds=True,label='strainer'), > PeriodicPythonRunner(virtPeriod=3e-7/strainRate,realPeriod=1, > > command='addPlotData()',label='plotDataCollector',initRun=True) > ] > > plot.plots={'eps':('sigma',)} > > O.saveTmp('initial'); > def addPlotData(): > > yade.plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress}) > > > _______________________________________________ > Mailing list: > https://launchpad.net/~yade-users<https://launchpad.net/%7Eyade-users> > Post to : [email protected] > Unsubscribe : > https://launchpad.net/~yade-users<https://launchpad.net/%7Eyade-users> > More help : https://help.launchpad.net/ListHelp > >
_______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
