ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

implement numbered equations and equation references

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-09-29 12:04:43 -04:00
parent 96c0d39be2
commit c6eb6d858b
2 changed files with 53 additions and 25 deletions
Download Patch File Download Diff File Expand all files Collapse all files

70
doc/src/Howto_peri.rst
View File

@ -130,13 +130,15 @@ notation of :ref:`(Silling 2007) <Silling2007_2>`, we write the
peridynamic equation of motion as
.. math::
:label: periNewtonII
\rho(\mathbf{x}) \ddot{\mathbf{u}}(\mathbf{x},t) =
\int_{\mathcal{H}_{\mathbf{x}}} \left\{ \underline{\mathbf{T}}\left[
\mathbf{x},t \right]\left<\mathbf{x}^{\prime}-\mathbf{x} \right>
- \underline{\mathbf{T}}\left[\mathbf{x}^{\prime},t \right]\left<\mathbf{x}-\mathbf{x}^{\prime} \right> \right\}
{d}V_{\mathbf{x}^\prime} + \mathbf{b}(\mathbf{x},t),
\int_{\mathcal{H}_{\mathbf{x}}} \left\{
\underline{\mathbf{T}}\left[\mathbf{x},t
\right]\left<\mathbf{x}^{\prime}-\mathbf{x} \right> -
\underline{\mathbf{T}}\left[\mathbf{x}^{\prime},t
\right]\left<\mathbf{x}-\mathbf{x}^{\prime} \right> \right\}
{d}V_{\mathbf{x}^\prime} + \mathbf{b}(\mathbf{x},t),
where :math:`\rho` represents the mass density,
:math:`\underline{\mathbf{T}}` the force vector state, and
@ -147,7 +149,7 @@ where :math:`\rho` represents the mass density,
radius :math:`\delta>0` centered at :math:`\mathbf{x}`. :math:`\delta`
is called the *horizon*, and is analogous to the cutoff radius used in
molecular dynamics. Conditions on :math:`\underline{\mathbf{T}}` for
which \eqref{eqn:NewtonII} satisfies the balance of linear and angular
which :math:numref:`periNewtonII` satisfies the balance of linear and angular
momentum are given in :ref:`(Silling 2007) <Silling2007_2>`.
We consider only force vector states that can be written as
@ -161,6 +163,7 @@ with :math:`\underline{t}` a *scalar force state* and
state*, defined by
.. math::
:label: periM
\underline{\mathbf{M}}\left< \boldsymbol{\xi} \right> =
\left\{ \begin{array}{cl}
@ -192,6 +195,7 @@ along with the horizon :math:`\delta`.
The LPS model has a force scalar state
.. math::
:label: periLPSState
\underline{t} = \frac{3K\theta}{m}\underline{\omega}\,\underline{x} +
\alpha \underline{\omega}\,\underline{e}^{\rm d},
@ -210,6 +214,7 @@ the reference position scalar state :math:`\underline{x}` so that
defined as
.. math::
:label: periWeightedVolume
m\left[ \mathbf{x} \right] = \int_{\mathcal{H}_\mathbf{x}}
\underline{\omega} \left<\boldsymbol{\xi}\right>
@ -291,6 +296,7 @@ and the horizon :math:`\delta`.
The PMB model is expressed using the scalar force state field
.. math::
:label: periPMBState
\underline{t}\left[ \mathbf{x},t \right]\left< \boldsymbol{\xi} \right> = \frac{1}{2} f\left( \boldsymbol{\eta} ,\boldsymbol{\xi} \right),
@ -298,17 +304,21 @@ with :math:`f` a scalar-valued function. We assume that :math:`f` takes
the form
.. math::
f = c s,
where
.. math::
:label: peric
c = \frac{18K}{\pi \delta^4},
with :math:`K` the bulk modulus and :math:`\delta` the horizon, and
:math:`s` the bond stretch, defined as
.. math::
:label: peristretch
s(t,\mathbf{\eta},\mathbf{\xi}) = \frac{ \left\Vert {\mathbf{\eta}+\mathbf{\xi}} \right\Vert - \left\Vert {\mathbf{\xi}} \right\Vert }{\left\Vert {\mathbf{\xi}} \right\Vert}.
@ -320,9 +330,10 @@ appropriate for 3D models only. For more on the origins of the constant
:math:`c`, see :ref:`(Silling 2005) <Silling2005>`. For the derivation
of :math:`c` for 1D and 2D models, see :ref:`(Emmrich) <Emmrich2007>`.
Given \eqref{eqn:PMBState}, \eqref{eqn:NewtonII} reduces to
Given :math:numref:`periPMBState`, :math:numref:`periNewtonII` reduces to
.. math::
:label: periPMBEOM
\rho(\mathbf{x}) \ddot{\mathbf{u}}(\mathbf{x},t) = \int_{\mathcal{H}_\mathbf{x}} \mathbf{f} \left( \boldsymbol{\eta},\boldsymbol{\xi} \right){d}V_{\boldsymbol{\xi}} + \mathbf{b}(\mathbf{x},t),
@ -352,6 +363,7 @@ Define :math:`\mu` to be the history-dependent scalar
boolean function
.. math::
:label: perimu
\mu(t,\mathbf{\eta},\mathbf{\xi}) = \left\{
\begin{array}{cl}
@ -366,6 +378,7 @@ where :math:`\mathbf{\eta}^\prime = \textbf{u}(\textbf{x}^{\prime
critical stretch defined as
.. math::
:label: peris0
s_0(t,\mathbf{\eta},\mathbf{\xi}) = s_{00} - \alpha s_{\min}(t,\mathbf{\eta},\mathbf{\xi}), \qquad s_{\min}(t) = \min_{\mathbf{\xi}} s(t,\mathbf{\eta},\mathbf{\xi}),
@ -388,6 +401,7 @@ Following :ref:`(Silling) <Silling2005>`, we can define the damage at a
point :math:`\textbf{x}` as
.. math::
:label: peridamage
\varphi(\textbf{x}, t) = 1 - \frac{\int_{\mathcal{H}_{\textbf{x}}} \mu(t,\mathbf{\eta},\mathbf{\xi}) dV_{\textbf{x}^\prime}
}{ \int_{\mathcal{H}_{\textbf{x}}} dV_{\textbf{x}^\prime} }.
@ -395,7 +409,7 @@ point :math:`\textbf{x}` as
Discrete Peridynamic Model and LAMMPS Implementation
""""""""""""""""""""""""""""""""""""""""""""""""""""
In LAMMPS, instead of \eqref{eqn:NewtonII}, we model this equation of
In LAMMPS, instead of :math:numref:`periNewtonII`, we model this equation of
motion:
.. math::
@ -409,7 +423,7 @@ where we explicitly track and store at each timestep the positions and
not the displacements of the particles. We observe that
:math:`\ddot{\textbf{y}}(\textbf{x}, t) = \ddot{\textbf{x}} +
\ddot{\textbf{u}}(\textbf{x}, t) = \ddot{\textbf{u}}(\textbf{x}, t)`, so
that this is equivalent to \eqref{eqn:NewtonII}.
that this is equivalent to :math:numref:`periNewtonII`.
Spatial Discretization
^^^^^^^^^^^^^^^^^^^^^^
@ -422,12 +436,14 @@ denote the family of particles for which particle :math:`i` shares a
bond in the reference configuration. That is,
.. math::
:label: periBondFamily
\mathcal{F}_i = \{ p ~ | ~ \left\Vert {\textbf{x}_p - \textbf{x}_i} \right\Vert \leq \delta \}.
The discretized equation of motion replaces \eqref{eqn:NewtonII} with
The discretized equation of motion replaces :math:numref:`periNewtonII` with
.. math::
:label: peridiscreteNewtonII
\rho \ddot{\textbf{y}}_i^n =
\sum_{p \in \mathcal{F}_i}
@ -444,25 +460,29 @@ In the model discussed so far, particles interact only through their
bond forces. A particle with no bonds becomes a free non-interacting
particle. To account for contact forces, short-range forces are
introduced :ref:`(Silling 2007) <Silling2007_2>`. We add to the force in
\eqref{eqn:discreteNewtonII} the following force
:math:numref:`peridiscreteNewtonII` the following force
.. math::
:label: perifS
\textbf{f}_S(\textbf{y}_p,\textbf{y}_i) = \min \left\{ 0, \frac{c_S}{\delta}(\left\Vert {\textbf{y}_p-\textbf{y}_i} \right\Vert - d_{pi}) \right\}
\frac{\textbf{y}_p-\textbf{y}_i}{\left\Vert {\textbf{y}_p-\textbf{y}_i} \right\Vert},
where :math:`d_{pi}` is the short-range interaction distance between
particles :math:`p` and :math:`i`, and :math:`c_S` is a multiple of the
constant :math:`c` from \eqref{eqn:c}. Note that the short-range force
constant :math:`c` from :math:numref:`peric`. Note that the short-range force
is always repulsive, never attractive. In practice, we choose
.. math::
:label: pericS
c_S = 15 \frac{18K}{\pi \delta^4}.
For the short-range interaction distance, we choose :ref:`(Silling 2007)
<Silling2007_2>`
.. math::
:label: perid
d_{pi} = \min \left\{ 0.9 \left\Vert {\textbf{x}_p - \textbf{x}_i} \right\Vert, 1.35 (r_p + r_i) \right\},
@ -488,16 +508,16 @@ short-range family of particles
Modification to the Particle Volume
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The right-hand side of \eqref{eqn:discreteNewtonII} may be thought of as
a midpoint quadrature of \eqref{eqn:NewtonII}. To slightly improve the
The right-hand side of :math:numref:`peridiscreteNewtonII` may be thought of as
a midpoint quadrature of :math:numref:`periNewtonII`. To slightly improve the
accuracy of this quadrature, we discuss a modification to the particle
volume used in \eqref{eqn:discreteNewtonII}. In a situation where two
volume used in :math:numref:`peridiscreteNewtonII`. In a situation where two
particles share a bond with :math:`\left\Vert { \textbf{x}_p -
\textbf{x}_i }\right\Vert = \delta`, for example, we suppose that only
approximately half the volume of each particle is "seen" by the other
:ref:`(Silling 2007) <Silling2007>`. When computing the force of each
particle on the other we use :math:`V_p / 2` rather than :math:`V_p` in
\eqref{eqn:discreteNewtonII}. As such, we introduce a nodal volume
:math:numref:`peridiscreteNewtonII`. As such, we introduce a nodal volume
scaling function for all bonded particles where :math:`\delta - r_i \leq
\left\Vert { \textbf{x}_p - \textbf{x}_i } \right\Vert \leq \delta` (see
the Figure below).
@ -555,12 +575,12 @@ Breaking Bonds
During the course of simulation, it may be necessary to break bonds, as
described in the :ref:`Damage section <peridamage>`. Bonds are recorded
as broken in a simulation by removing them from the bond family
:math:`\mathcal{F}_i` (see \eqref{eqn:BondFamily}).
:math:`\mathcal{F}_i` (see :math:numref:`periBondFamily`).
A naive implementation would have us first loop over all bonds and
compute :math:`s_{min}` in \eqref{eqn:s0}, then loop over all bonds
compute :math:`s_{min}` in :math:numref:`peris0`, then loop over all bonds
again and break bonds with a stretch :math:`s > s0` as in
\eqref{eqn:mu}, and finally loop over all particles and compute forces
:math:numref:`perimu`, and finally loop over all particles and compute forces
for the next step of :ref:`Algorithm 1 <algvelverlet>`. For reasons of
computational efficiency, we will utilize the values of :math:`s_0` from
the *previous* timestep when deciding to break a bond.
@ -602,7 +622,7 @@ routines in :ref:`Algorithm 3 <algperilpsm>` and :ref:`Algorithm 4
| **for all** particles :math:`j \in \mathcal{F}^S_i` (the short-range family of nodes for particle :math:`i`) **do**
| :math:`r = \left\Vert {\textbf{y}_j - \textbf{y}_i} \right\Vert`
| :math:`dr = \min \{ 0, r - d \}` {*Short-range forces are only repulsive, never attractive*}
| :math:`k = \frac{c_S}{\delta} V_k dr` {:math:`c_S` *defined in \eqref{eqn:cS}*}
| :math:`k = \frac{c_S}{\delta} V_k dr` {:math:`c_S` *defined in :math:numref:`pericS`*}
| :math:`\textbf{f} = \textbf{f} + k \frac{\textbf{y}_j-\textbf{y}_i}{\left\Vert {\textbf{y}_j-\textbf{y}_i} \right\Vert}`
| **end for**
| **end for**
@ -675,7 +695,7 @@ A sketch of the PMB model implementation in the PERI package appears in
| **for all** particles :math:`j \in \mathcal{F}^S_i` (the short-range family of nodes for particle :math:`i`) **do**
| :math:`r = \left\Vert {\textbf{y}_j - \textbf{y}_i} \right\Vert`
| :math:`dr = \min \{ 0, r - d \}` {*Short-range forces are only repulsive, never attractive*}
| :math:`k = \frac{c_S}{\delta} V_k dr` {:math:`c_S` *defined in \eqref{eqn:cS}*}
| :math:`k = \frac{c_S}{\delta} V_k dr` {:math:`c_S` *defined in :math:numref:`pericS`*}
| :math:`\textbf{f} = \textbf{f} + k \frac{\textbf{y}_j-\textbf{y}_i}{\left\Vert {\textbf{y}_j-\textbf{y}_i} \right\Vert}`
| **end for**
| **end for**
@ -684,7 +704,7 @@ A sketch of the PMB model implementation in the PERI package appears in
| **for all** particles :math:`j` sharing an unbroken bond with particle :math:`i` **do**
| :math:`r = \left\Vert {\textbf{y}_j - \textbf{y}_i} \right\Vert`
| :math:`dr = r - \left\Vert {\textbf{x}_j - \textbf{x}_i} \right\Vert`
| :math:`k = \frac{c}{\left\Vert {\textbf{x}_i - \textbf{x}_j} \right\Vert} \nu(\textbf{x}_i - \textbf{x}_j) V_j dr` {:math:`c` *defined in \eqref{eqn:c}*}
| :math:`k = \frac{c}{\left\Vert {\textbf{x}_i - \textbf{x}_j} \right\Vert} \nu(\textbf{x}_i - \textbf{x}_j) V_j dr` {:math:`c` *defined in :math:numref:`peric`*}
| :math:`\textbf{f} = \textbf{f} + k \frac{\textbf{y}_j-\textbf{y}_i}{\left\Vert {\textbf{y}_j-\textbf{y}_i} \right\Vert}`
| **if** :math:`(dr / \left\Vert {\textbf{x}_j - \textbf{x}_i} \right\Vert) > \min(s_0(i), s_0(j))` **then**
| Break :math:`i`'s bond with :math:`j` {:math:`j`\ *'s bond with* :math:`i` *will be broken when this loop iterates on* :math:`j`}
@ -699,7 +719,7 @@ A sketch of the PMB model implementation in the PERI package appears in
Damage
^^^^^^
The damage associated with every particle (see \eqref{eqn:damage}) can
The damage associated with every particle (see :math:numref:`peridamage`) can
optionally be computed and output with a LAMMPS data dump. To do this,
your input script must contain the command :doc:`compute damage/atom
<compute_damage_atom>` This enables a LAMMPS per-atom compute to
@ -880,7 +900,7 @@ three times the lattice constant.) The target is a cylinder of diameter
lattice contains 103,110 particles. Each particle :math:`i` has volume
fraction :math:`V_i = 1.25 \times 10^{-10} \textrm{m}^3`.
The spring constant in the PMB material model is (see \eqref{eqn:c})
The spring constant in the PMB material model is (see :math:numref:`peric`)
.. math::
c = \frac{18k}{\pi \delta^4} = \frac{18 (14.9 \times 10^9)}{ \pi (1.5 \times 10^{-3})^4} \approx 1.6863 \times 10^{22}.
@ -927,7 +947,7 @@ the distance from the atom to the center of the indenter, and :math:`R`
is the radius of the projectile. The force is repulsive and :math:`F(r) =
0` for :math:`r > R`. For our problem, the projectile radius :math:`R =
0.05` m, and we have chosen :math:`k_s = 1.0 \times 10^{17}` (compare
with \eqref{eqn:c} above).
with :math:numref:`peric` above).
Writing the LAMMPS Input File
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

8
doc/utils/sphinx-config/_static/css/lammps.css
View File

@ -1,3 +1,11 @@
.math {
text-align: left;
}
.eqno {
float: right;
}
.wy-nav-content {
max-width: 100% !important;
}