Hello,
In the peridynamics PMB pair style source (pair_peri_pmb.cpp), when the scalar force state is calculated and added to the total force on the “node” it is first made negative:
pair_peri_pmb.cpp, line 242
“”"
stretch = dr / r0[i][jj];
rk = (kspring[itype][jtype] * vfrac[j]) * vfrac_scale * stretch;
if (r > 0.0) fbond = -(rk/r);
“”"
Here stretch is the strain (amount stretched from equilibrium), kspring is the spring constant, fbond is the bond force, and r is the magnitude of the total stretch
My question: why is fbond the negative of rk? Is this simply Hooke’s law: F=-kx?
Thanks
Hello,
In the peridynamics PMB pair style source (pair_peri_pmb.cpp), when the
scalar force state is calculated and added to the total force on the "node"
it is first made negative:
pair_peri_pmb.cpp, line 242
"""
stretch = dr / r0[i][jj];
rk = (kspring[itype][jtype] * vfrac[j]) * vfrac_scale * stretch;
if (r > 0.0) fbond = -(rk/r);
"""
Here stretch is the strain (amount stretched from equilibrium), kspring is
the spring constant, fbond is the bond force, and r is the magnitude of the
total stretch
My question: why is fbond the negative of rk? Is this simply Hooke's law:
F=-kx?
this question is rather pointless without also stating, whether this
is f_AB or f_BA. remember newton's third law?
axel.
Axel,
I understand you point. I should have been more specific but I’m trying to learn a pretty large code base so cut me a little slack 
So, to be more specific, here’s the full code:
“”"
for (i = 0; i < nlocal; i++) {
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jnum = npartner[i];
s0_new[i] = DBL_MAX;
first = true;
for (jj = 0; jj < jnum; jj++) {
if (partner[i][jj] == 0) continue;
j = atom->map(partner[i][jj]);
// check if lost a partner without first breaking bond
if (j < 0) {
partner[i][jj] = 0;
continue;
}
// compute force density, add to PD equation of motion
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
if (periodic) domain->minimum_image(delx,dely,delz);
rsq = delxdelx + delydely + delz*delz;
jtype = type[j];
delta = cut[itype][jtype];
r = sqrt(rsq);
dr = r - r0[i][jj];
// avoid roundoff errors
if (fabs(dr) < 2.2204e-016) dr = 0.0;
// scale vfrac[j] if particle j near the horizon
if ((fabs(r0[i][jj] - delta)) <= half_lc)
vfrac_scale = (-1.0/(2half_lc))(r0[i][jj]) +
(1.0 + ((delta - half_lc)/(2*half_lc) ) );
else vfrac_scale = 1.0;
stretch = dr / r0[i][jj];
rk = (kspring[itype][jtype] * vfrac[j]) * vfrac_scale * stretch;
if (r > 0.0) fbond = -(rk/r);
else fbond = 0.0;
f[i][0] += delxfbond;
f[i][1] += delyfbond;
f[i][2] += delz*fbond;
“”"
So, F=-k*displacement?
Axel,
I understand you point. I should have been more specific but I'm trying to
learn a pretty large code base so cut me a little slack
if you would be looking at a complex piece of code (and LAMMPS has
those, too), yes, i would be willing to cut you some more slack, but
this one is not a complex one.
So, to be more specific, here's the full code:
[...]
So, F=-k*displacement?
this isn't really more specific. which displacement of what are you
talking about here?
as my (final) advice, here is what i can read from the code, purely
based on looking at the code and some educated guesses. please note,
that my understanding of the peridynamics method itself is extremely
limited. however, this code seems rather straightforward:
for each atom i, you loop over its "partners", i.e. atoms that were
neighbors of i at the beginning of the calculation. if a partner
doesn't exist anymore, you continue with the next.
then you compute the (current) distance r between atoms i and j and
compute the difference dr between the current distance and the
distance at the beginning (r0). from this you compute the relative
extension or compression, stretch, and multiply it with the force
constant. finally you project on the x, y and z components via delx/r
dely/r and delz/r. the /r term is applied ahead of time for obvious
computational efficiency. as for the sign of fbond, you must keep in
mind in mind that dr can be both positive or negative depending on
whether the distance between atoms i and j has been reduced or grown.
from here on, i strongly suggest you take a piece of paper and make a
sketch of the situation and see what the sign of the (restoring) force
should be.
reading code correctly requires more than just matching individual
lines with a formula, that you recognize from somewhere.
axel.