calculation of flexural stress in 2D condition

Hi Steve,

Please note that 2D solutions are used either for plane streets or plain strain. In such a case the load is actually 'linear pressure' and the force is distributed over assumed thickness (e.g. unitary). Thanks to this, you can you use your formula, by assuming the thickness is known. In may opinion there is no good answer for your question. I would suggest to take step back and rethink whether 2D is a good idea.
If you want to stick to 2D you could also 'measure' forces in the symmetry o plane and 'smear' them over a known cross-section of the sample that you are modeling. By smearing I mean you could calculate the bending moment and apply it to the theoretical beam.

Best wishes,
Karol

Thanks Karol,
Do you mean that the force (f) I obtained from 2D is already a linear load and so that the total force over the assumed thickness (b) would be f*b, which makes the formula become sigma=f*b*L/(b*d^2)=f*L/d^2 ?
And I don't really understand the second part like what is the symmetry o plane, could you give more elaboration on this?
Best regards,
Steve

Hi,

sorry, forget the second part. It would be good for the statically indeterminate beam. Since you have a four-point bending test, you know the moment already.

Regarding the first part, I meant that in actual 2D formulation you do not have to worry about thickness.

The load that you obtain (apply?) in the simulation has a force dimension, thus I think it is related to some thickness. Whatever thickness you assume, you should treat the force as integral over this thickness. For example, if your sample has b = 40 mm, and your model has b' = 1 mm, you should apply only f'=f*b'/b = f/40.
In general, b' doesn't have to correspond to the thickness of your sample, you could also assume b'=b. This would mean that your sample is 'compressed' to the thickness of the model. In such a case, however, the stiffness considerations become harder.

In my opinion, assuming thickness to solve this problem is quite intuitive, and could be a good solution in some cases. However, it is probably oversimplified. It would be useful to know the purpose of the simulation and what is the advantage of using DEM? Is this advantage affected by the 2D assumption?

Besides all that 'fundamental questions', I think that both solutions proposed by you are quite decent (if you stick to 2D simulation ;) ). Only, remember to apply the force corresponding to the assumed thickness.

Cheers,
Karol

Thanks Karol,
I am sorry if I my statement was confusing, let me just explain it again. So I've already obtained the flexural stress & strain form the experiments with a known dimension b=100mm. And then I construct a 2D model in yade in order to simulate this experiments and see the crack pattern. However, during the calibration process, I have some problem in matching the 2D simulation results to the experiment data. Actually, I am not applying force in both experiment and simulation, but control the loading rate of the loading plate and obtain the force that acting on the plates. So in this case, is the force obtained in the simulation supposed to multiply a scale 100/b' (100 is the thickness in the experiment, b' is the thickness of the model) to restore the force to the 3D condition and then use the formula to calculate the stress? If so, then what would be the thickness of the model? is it the average diameter of the particle? or is there no thickness at all because it is 2D simulation and both force and displacement in the out-of -plane direction has been restricted?
And I test both methods in simulation which are by formula sigma=fL/d^2 and by directly get from particles:
###
sigma=bodyStressTensors()
tensile_stress=sum(sigma[i][0,0] for i in p) ## p is the particle ids located in the bottom of the beam between two loading point
###
Although the shape of curves are similar, but the peak flexural stress obtained directly from particles are over two times larger than that of formula. I am really confused where this difference come from if both methods are good enough.
Best regards,
Steve

> both force and displacement in the out-of -plane direction has been restricted?

This is not possible. In all mechanics, you know/prescribe/restrict exactly one quantity of the couple force/displacement (or some related quantity, like stress/strain). The other quantity is unknown and as such cannot be restricted. Of course it can be computed, but is result of the system.
In plane stress, you know the out-of-plane stress is zero, but you cannot say anything about out-of-plane strain.
In plane strain, you know the out-of-plane strain is zero, but you cannot say anything about out-of-plane stress.
....

> If so, then what would be the thickness of the model? is it the average diameter of the particle?

sounds reasonable

> or is there no thickness at all because it is 2D simulation

this is another approach. Then your model is "per unit thickness".
So for transfer model -> experiment, multiply your model data y real beam thickness.
Or for experiment -> model, divide real data by beam thickness

> the peak flexural stress obtained directly from particles are over two times larger than that of formula. I am really confused where this difference come from if both methods are good enough.

Have you done the "porosity correction" for bodyStressTensors [1]?

cheers
Jan

[1] https://yade-dem.org/doc/yade.utils.html#yade._utils.bodyStressTensors

Hi, Jan
Thanks for the correction. I am a bit unclear why there can be other approaches for the transfer. I mean the thickness of the model must be either 0 or some specific value(eg. the average diameter). So the unit of the force I obtained from the simulation is Newton per unit thickness or just newton?
And yes I forgot the porosity correction. But after adding the porosity which is 1-(area of the all particles/ sample area) in 2D, and multiply stress by the 1-porosity, the value is still 1.7 times larger than from formula since the porosity is only around 0.2. Please correct me if there's something wrong.
Best regards,
Steve

> I am a bit unclear why there can be other approaches for the transfer

it is a model. For the same real situation, you can have many different models...

> I mean the thickness of the model must be either 0 or some specific value(eg. the average diameter). So the unit of the force I obtained from the simulation is Newton per unit thickness or just newton?

Yes.
For real thickness, you get real forces, Newtons.
For "zero" thickness (I would not call the thickness zero, but rather something like planar or 2D case, since the thickness is not zero, it is just not important for the solution itself), you get force per unit thickness, Newton per unit thickness.

> But after adding the porosity which is 1-(area of the all particles/ sample area) in 2D
> the porosity is only around 0.2

"area" porosity, or "volume" porosity?
As I mentioned, bodyStressTensors are for 3D, are used (and assumes working with) with volumes, not areas.
If you have significant differences, in particle sizes, "area" correction is wrong, since the particle "volume" and "area" are not proportional.
Maybe you can use the same math theory, only adjusted to the 2D case.

For (almost) static case, the "stress" approach should fit well with the theory.

Anyway, I think we came to the place, where a MWE [2] is needed for reasonable continuation.
If you prepare a code, please focus on the "M=minimal" part (both from code size and running time point of view)

cheers
Jan

[2] https://www.yade-dem.org/wiki/Howtoask

Hi Jan,
Sorry that I thought the simulation was 2D so I used the 2D porosity. But I forgot that stress is also obtained in 3D.
1. After correction, the peak values do become closer to each other, but there is still some gap as well as the shape.
2. the error in volume calculation is that I didn't use aabbExtrema() to calculate the volume because if I add cement particle before the roller(moving line78,79 to line 69) , an error called 'Assertion `!swap' failed' occurred.
3. I tried the TesslationWrap() but the value is way much smaller, maybe it's due to the reason that the function also considers the rollers and supports?
line 1-59 is just making 2D packing, line 60-156 is 4-point bending model
Many thanks,
Steve

#################code attached
## compact 2D packing
from yade import pack,plot
nRead=readParamsFromTable(
 num_spheres=2000,
 compFricDegree = 30,
 key='_triax_base_',
 unknownOk=True
)
from yade.params import table
num_spheres=table.num_spheres
key=table.key
targetPorosity = 0.55
compFricDegree = table.compFricDegree
finalFricDegree = 30
rate=-0.02
damp=0.2
stabilityThreshold=0.01
young=5e6
mn,mx=Vector3(0,-.1,0),Vector3(0.35,.1,.1)
O.switchScene();O.resetThisScene()
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([(0,-.1,0),(0.35,.1,.1)],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(0.35,0,.1),-1,0.3,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = 0,
 stressMask = 7,
 internalCompaction=True,
)
newton=NewtonIntegrator(damping=damp)
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 newton
]
triax.goal1=triax.goal3=-500
triax.goal2=0
while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print ('unbalanced force:',unb,' mean stress: ',triax.meanStress)
  if unb<stabilityThreshold and abs(-1000/3-triax.meanStress)/1000/3<0.001:
    break
sp=SpherePack(); sp.fromSimulation()
O.switchScene()
## constuct 4-point bending model
import numpy as np
from yade import pack,plot
## define parameters
strainRate=-.02
sphereRadius=.0005
young=2.8e9
sigmaT=1e6
frictionAngle=atan(0.6)
relDuctility=17
omegaThreshold=0.999
epsCrackOnset=.004
damping=.6
damLaw = 0
relFuzz=0.3
factor=1.5
isoPrestress=0
poisson= 0.1
pred=pack.inAlignedBox((-0.4,-1,-1),(0.4,1,1))
idConcrete=O.materials.append(CpmMat(young=young,poisson=poisson,frictionAngle=frictionAngle,density=3120,sigmaT=sigmaT,epsCrackOnset=epsCrackOnset,isoPrestress=isoPrestress,relDuctility = relDuctility,damLaw = damLaw))
#

idRoller=O.materials.append(FrictMat(young=9e9,frictionAngle=atan(.5),poisson=.3,density=3e3,label='walls'))
## create top roller clumps
TL=O.bodies.append([sphere(center=(.125-0.175,0,.115-0.05+0.0002),radius=.015,material=idRoller)])
TR=O.bodies.append([sphere(center=(.225-0.175,0,.115-0.05+0.0002),radius=.015,material=idRoller)])
## creat bottom roller
O.bodies.append([sphere(center=(.025-0.175,0,-.015-0.05),radius=.015,material=idRoller,fixed=True)])
O.bodies.append([sphere(center=(.325-0.175,0,-.015-0.05),radius=.015,material=idRoller,fixed=True)])
#
spp=filterSpherePack(pred,sp,material=idConcrete,returnSpherePack=True)
cement=spp.toSimulation(material=idConcrete)
#### calculate the porosity for volume
volume=0
for i in cement:
  volume+=4/3*pi*pow(O.bodies[i].shape.radius,3)
porosity=1-volume/(0.3*0.0045*0.1)
print(porosity)
T=[O.bodies[i] for i in cement if (O.bodies[i].state.pos[2]<(-.05+2*O.bodies[i].shape.radius) and O.bodies[i].state.pos[0]>(-2*O.bodies[i].shape.radius) and O.bodies[i].state.pos[0]<(2*O.bodies[i].shape.radius))]
initial_dis= T[0].state.pos[0]-T[1].state.pos[0]
####grab the particle at the bottom fibre
p=[]
for i in cement:
  if (O.bodies[i].state.pos[2]<(-.05+2*O.bodies[i].shape.radius) and O.bodies[i].state.pos[0]>-.05 and O.bodies[i].state.pos[0]<.05):
    p.append(i)

O.engines=[
  ForceResetter(),
  InsertionSortCollider([
    Bo1_Sphere_Aabb(aabbEnlargeFactor=factor,label='bo1s'),
  ]),
  InteractionLoop(
    [
    Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=factor,label='ig2ss'),
    ],
    [
    Ip2_CpmMat_CpmMat_CpmPhys(),
    Ip2_FrictMat_CpmMat_FrictPhys(),
    Ip2_FrictMat_FrictMat_FrictPhys()
    ],
    [
    Law2_ScGeom_CpmPhys_Cpm(),
    Law2_ScGeom_FrictPhys_CundallStrack()
    ]
  ),
  GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.3),
  NewtonIntegrator(damping=0.5),
  CpmStateUpdater(realPeriod=.5),
  PyRunner(command='initiateStrain()',iterPeriod=1,label='move'),
  PyRunner(command='finish()',realPeriod=5),
  PyRunner(command='addPlotData()',iterPeriod=100),
]
O.bodies[TL[0]].state.blockedDOFs='xyzXYZ'
O.bodies[TR[0]].state.blockedDOFs='xyzXYZ'
O.step()
bo1s.aabbEnlargeFactor=1
ig2ss.interactionDetectionFactor=1
def initiateStrain():
   O.bodies[TL[0]].state.vel[2]=strainRate
   O.bodies[TR[0]].state.vel[2]=strainRate
def finish():
  if (O.bodies[TL[0]].state.refPos[2]-O.bodies[TL[0]].state.pos[2])*60/13>0.03:
     O.pause()
     pass
def addPlotData():
  f_up=utils.sumForces(ids=TL+TR,direction=(0,0,-1))
  ### directly obtain stress on the bottom fibre
  TW=TesselationWrapper()
  TW.setState()
  TW.computeVolumes()
  sigma=bodyStressTensors()
  tensile_stress=sum(sigma[i][0,0] for i in p)
  tensile_stress2= sum(sigma[i][0,0]*4/3*pi*O.bodies[i].shape.radius**3/TW.volume(i) for i in p)
  plot.addData(s=-f_up*0.3/(.01*.01),e2=60/13*(O.bodies[TL[0]].state.refPos[2]-O.bodies[TL[0]].state.pos[2]),s2=tensile_stress*(1-porosity)/len(p),s3=tensile_stress2/len(p)) ## e2 is the tensile strain derived from the load point displacement
plot.plots={'e2':'s','e2 ':'s2','e2 ':'s3'} ## s is from formula fl/d^2 which neglect the width, s2 is directly obtained from bottom fibre and scaled by factor, s3 is using TesslationWrap?
plot.plot()
O.run()