timprove documentation - slidergrid - grid of elastic sliders on a frictional surface
(HTM) git clone git://src.adamsgaard.dk/slidergrid
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
(DIR) commit 01f9042427f6b9e5933c803f55b9fdafe5b4a33c
(DIR) parent 897289bb9cc3233e556ccd4b863cdfaf4dea29b0
(HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Wed, 4 May 2016 13:10:02 -0700
improve documentation
Diffstat:
M doc/doc.pdf | 0
M doc/doc.tex | 295 +++++++++++++++----------------
M slidergrid/grid.c | 29 +++++++++++++++++++++++++++++
M slidergrid/slider.h | 1 +
4 files changed, 175 insertions(+), 150 deletions(-)
---
(DIR) diff --git a/doc/doc.pdf b/doc/doc.pdf
Binary files differ.
(DIR) diff --git a/doc/doc.tex b/doc/doc.tex
t@@ -19,7 +19,7 @@
\usepackage[T1]{fontenc} % Font encoding
\usepackage{charter} % Serif body font
\usepackage[charter]{mathdesign} % Math font
-\usepackage[scale=0.9]{sourcecodepro} % Monospaced fontenc
+\usepackage[scale=0.8]{sourcecodepro} % Monospaced fontenc
\usepackage[lf]{FiraSans} % Sans-serif font
\usepackage{listings}
t@@ -48,15 +48,24 @@ of Oceanography\\University of California, San Diego}\\[3mm] Last revision:
\maketitle
\section{Methods}
-The method is derived from \citet{Schlangen1996}, \citet{Radjai2011} and
-\citet{Potyondy2004} but is, relative to the cited works, adapted for three
-spatial dimensions and non-linear properties.
-
-The Lagrangian nodes are connected with visco-elastic beam elements which are
-resistive to relative translational or rotational movement. The kinematic
-degrees of freedom are determined by explicit integration of Newton's second law
-of motion for translation and rotation. For a point $i$ with bonded
-interactions to nodes $j\in N_c$, the translational accelerations
+Our approach treats the short-temporal scale behavior of stick-slip as a
+rigid-body dynamics problem. The material is represented as a discrete number
+of Lagrangian points (\emph{nodes}) which are mechanically interacting with each
+other and the boundary conditions.
+
+The Lagrangian nodes are connected with visco-elastic beam elements.
+The bonds are resistive to tension and compression, shearing, twisting, and
+bending, which ensures elastic uniformity regardless of geometric node
+arrangement \citep{Bolander1998, Radjai2011}. Alternatively, the node
+interaction could be parameterized as simple springs which exclusively provide
+resistance to tension and compression, and resistance to shearing, bending and
+twisting would be introduced by discretizing the elastic material into an
+irregular network of many more nodes \citep[e.g.][]{Topin2007, Topin2009}. Here
+we chose the former approach which allows us to keep the number of nodes low.
+
+The kinematic degrees of freedom are determined by explicit integration of
+Newton's second law of motion for translation and rotation. For a point $i$
+with bonded interactions to nodes $j\in N_c$, the translational accelerations
($\boldsymbol{a}$) are found from the sums of forces:
\begin{equation}
\boldsymbol{a}_i =
t@@ -85,24 +94,85 @@ torques:
\boldsymbol{\alpha}_i =
\sum^{N_c}_j
\left[
- \frac{\boldsymbol{t}^\text{s}_{i,j}}{I_i} +
- \frac{\boldsymbol{t}^\text{t}_{i,j}}{J_{i,i}}
+ \frac{\boldsymbol{t}^{i,j}_{\bar{x}}}{I^\text{p}_{i,j}} +
+ \frac{\boldsymbol{t}^{i,j}_{\bar{y}}}{I^\text{n}_{i,j}} +
+ \frac{\boldsymbol{t}^{i,j}_{\bar{z}}}{I^\text{n}_{i,j}}
\right]
\label{eq:n2-ang}
\end{equation}
here, $\boldsymbol{t}^\text{s}$ is the torque resulting from shearing motion of
the bond, while the torque $\boldsymbol{t}^{t}$ results from relative twisting.
-$I_i$ is the local moment of inertia at the point, and $J_{i,j}$ is polar moment
-of inertia of the bond.
+$I^\text{n}_{i,j}$ is the bond-normal mass moment of inertia at the point, and
+$I^\text{p}$ is polar mass moment of inertia of the bond. The above equation
+implies the simplifying assumption that the nodes are bonded in a configuration
+with geometric symmetry, which is a good approximation inside the grid but
+slightly worse at the grid edges.
-At the beginning of each time step the accumulated strain on each inter-point
-bond is determined by considering the relative motion of the bonded nodes. The
-bond deformation is decomposed per kinematic degree of freedom, andis determined
-by an incremental method derived from \citet{Potyondy2004}. The strain can be
-decomposed into bond tension and compression, bond shearing, bond twisting, and
-bond bending. The accumulated strains are used to determine the magnitude of
-the forces and torques resistive to the deformation.
+\subsection{Visco-elastic interaction between nodes}
+The total force ($\boldsymbol{f}$) and torque ($\boldsymbol{t}$) on two nodes
+($i$ and $j$) is found by determining the stress response of a three-dimensional
+elastic Timoshenko beam to strain \citet{Schlangen1996, Austrell2004,
+ Aastroem2013}.
+In the following the forces and torques are described for node $i$ but are
+of equal magnitude and with opposite sign for node $j$. The interaction
+accounts for resistance to tension and compression, shear, torsion, and bending.
+The equations below are derived from the stiffness matrix in
+\citet{Austrell2004}. The components for the three-dimensional force vector on
+node $i$ are:
+\begin{equation}
+ \begin{split}
+ f_{\bar{x}}^i & = \frac{EA}{L}
+ \left( p_{\bar{x}}^{*,i} - p_{\bar{x}}^{*,j} \right)\\
+%
+ f_{\bar{y}}^i & = \frac{12EI_A}{L^3}
+ \left( p_{\bar{y}}^{*,i} - p_{\bar{y}}^{*,j} \right)
+ + \frac{6EI_A}{L^2}
+ \left( \Omega_{\bar{z}}^{*,i} + \Omega_{\bar{z}}^{*,j} \right)\\
+%
+ f_{\bar{z}}^i & = \frac{12EI_A}{L^3}
+ \left( p_{\bar{z}}^{*,i} - p_{\bar{z}}^{*,j} \right)
+ - \frac{6EI_A}{L^2}
+ \left( \Omega_{\bar{y}}^{*,i} + \Omega_{\bar{y}}^{*,j} \right)
+ \end{split}
+\end{equation}
+where $\bar{x}, \bar{y}, \bar{z}$ are the bond-relative coordinates.
+The linear and angular relative displacement of the nodes is described by $p^*$
+and $\Omega^*$. $E$ is Young's modulus, $G$ is the shear modulus, $A$ is the
+beam cross-sectional area, and $L$ is the original beam length. $I_A$ is the
+area moment of inertia of the beam ($I_A = a^4 12^{-1}$ where $a$ is the beam
+side length). The torque on node $i$ is:
+\begin{equation}
+ \begin{split}
+ t_{\bar{x}}^i & = \frac{GJ}{L}
+ \left( \Omega_{\bar{x}}^{*,i} + \Omega_{\bar{x}}^{*,j} \right)\\
+%
+ t_{\bar{y}}^i & = \frac{6EI_A}{L^2}
+ \left(p_{\bar{z}}^{*,j} - p_{\bar{z}}^{*,i} \right)
+ + \frac{4EI_A}{L}
+ \left( \Omega_{\bar{y}}^{*,i} + \frac{\Omega_{\bar{y}}^{*,j}}{2}
+ \right)\\
+%
+ t_{\bar{z}}^i & = \frac{6EI_A}{L^2}
+ \left(p_{\bar{y}}^{*,i} - p_{\bar{y}}^{*,j} \right)
+ + \frac{4EI_A}{L}
+ \left( \Omega_{\bar{z}}^{*,i} + \frac{\Omega_{\bar{z}}^{*,j}}{2}
+ \right)\\
+ \end{split}
+\end{equation}
+$GJ$ is the Saint-Venant torsional stiffness, where $J$ called the torsional
+rigidity multiplier or torsion constant. For a beam with a solid square
+cross-sectional shape it can be approximated as $J \approx 0.140577 a^4$ where
+$a$ is the side length \citep{Timoshenko1951, Roark1954, Weisstein2016}.
+
+% Torsional constant:
+% https://en.wikipedia.org/wiki/Torsion_constant
+% http://mathworld.wolfram.com/TorsionalRigidity.html (K: K_v)
+% http://physics.stackexchange.com/questions/83148/where-i-can-find-a-torsional-stiffness-table-for-different-types-of-stainless-st
+% St Venant torsion: K_v = 1/G (Austrell et al. 2004, table 3) Does it make sense?
+% chapter4.pdf: K_v = J
+% https://skyciv.com/free-moment-of-inertia-calculator/
+
The deformation and reactive forces are determined relative to the orientation
of the bond. Common geometrical vectors include the inter-distance vector
t@@ -137,82 +207,7 @@ orientation, e.g.\ the bond-parallel and bond-shear velocity, respectively:
\right)
\end{equation}
-The axial strain is the bond-parallel deformation and is determined as the
-change in inter-point length relative to the initial distance:
-\begin{equation}
- \epsilon_a = \frac{
- (\boldsymbol{d}_{i,j} - \boldsymbol{d}^0_{i,j}) \cdot n_{i,j}}
- {||\boldsymbol{d}^0_{i,j}||}
-\end{equation}
-The cross-sectional area of a bond ($A_{i,j}$) varies with axial strain
-($\epsilon_a$) scaled by Poissons ratio $\nu$:
-\begin{equation}
- A_{i,j} = A^0_{i,j}
- - A^0_{i,j}
- \left(
- 1 -
- \left(
- 1 + \epsilon_a
- \right)^{-\nu}
- \right)
-\end{equation}
-The mass of point $i$ is defined as the half of the mass of each of its bonds:
-\begin{equation}
- m_i = \frac{\rho}{2} \sum^{N_c}_j A^0_{i,j} ||\boldsymbol{d}^0_{i,j}||
-\end{equation}
-The density ($\rho$) is adjusted so that the total mass of all nodes matches the
-desired value.
-
-\subsection{Resistance to tension and compression}
-Bond tension and compression takes place when the relative translational
-distance between a pair of bonded nodes changes, and is the most important
-deformational mode in this model. The current axial strain is determined with a
-second-order central difference scheme. It is determined from the previous
-point positions and projected future positions:
-\begin{equation}
- \Delta d^t_{i,j} = \frac{d_{i,j}^{*,t+\Delta t} - d_{i,j}^{t-\Delta t}}{2}
-\end{equation}
-The future point distance in the above ($d_{i,j}^{*,t+\Delta}$) is found by
-applying a second-order Taylor expansion:
-\begin{equation}
- \boldsymbol{p}_i^{*,t+\Delta t} =
- \boldsymbol{p}_i^{t} +
- \boldsymbol{v}_i^{t} \Delta t +
- \frac{1}{2}\boldsymbol{a}_i^{t} \Delta t^2
-\end{equation}
-
-
-
-The bond-parallel force is determined from Young's modulus ($E$) and the
-cross-sectional area ($A_{i,j}$) of the bond:
-\begin{equation}
- \boldsymbol{f}^{i,j}_\text{p} =
- \frac{E A_{i,j}}{|| \boldsymbol{d}^0_{i,j} ||}
- \left(
- \boldsymbol{d}_{i,j} -
- \boldsymbol{d}^0_{i,j}
- \right)
-\end{equation}
-
-\subsection{Shear resistance}
-The bond-shear force is determined incrementally for the duration of the
-interaction:
-\begin{equation}
- \boldsymbol{f}^{i,j}_\text{s} = \int^t \Delta \boldsymbol{f}^{i,j}_\text{s}
- %\, dt
-\end{equation}
-where the increment in shear force is determined from the shear modulus ($G$),
-the cross-sectional area ($A_{i,j}$) of the bond, and the
-\begin{equation}
- \Delta \boldsymbol{f}^{i,j}_\text{s} =
- \frac{G A_{i,j}}{||\boldsymbol{d}^{i,j}||}
- \Delta \boldsymbol{d}^{i,j}_\text{s}
-\end{equation}
-
-\subsection{Twisting resistance}
-
-\subsection{Bending resistance}
\subsection{Temporal integration}
Once the force and torque sum components at time $t$ have been determined, the
t@@ -278,71 +273,71 @@ stiffness in the system relative to the smallest mass:
where $\epsilon$ is a safety factor related to the geometric structure of the
bonded network. We use $\epsilon = 0.07$.
-The total force ($\boldsymbol{f}$) and torque ($\boldsymbol{t}$) on two nodes
-($i$ and $j$) with translational ($\boldsymbol{p}$) and angular
-($\boldsymbol{\Omega}$) positions interconnected with a three-dimensional
-elastic beam can be expressed as the following set of equations. The
-interaction accounts for resistance to tension and compression, shear, torsion,
-and bending. The symmetrical matrix on the right hand side constitutes the
-\emph{stiffness matrix} \citep{Schlangen1996, Austrell2004}:
+
+\appendix
+\section{Stiffness matrix}
\begin{equation}
\begin{bmatrix}
- f_\text{x}^i\\[0.6em]
- f_\text{y}^i\\[0.6em]
- f_\text{z}^i\\[0.6em]
- t_\text{x}^i\\[0.6em]
- t_\text{y}^i\\[0.6em]
- t_\text{z}^i\\[0.6em]
- f_\text{x}^j\\[0.6em]
- f_\text{y}^j\\[0.6em]
- f_\text{z}^j\\[0.6em]
- t_\text{x}^j\\[0.6em]
- t_\text{y}^j\\[0.6em]
- t_\text{z}^j\\
+ f_{\bar{x}}^i\\[0.6em]
+ f_{\bar{y}}^i\\[0.6em]
+ f_{\bar{z}}^i\\[0.6em]
+ t_{\bar{x}}^i\\[0.6em]
+ t_{\bar{y}}^i\\[0.6em]
+ t_{\bar{z}}^i\\[0.6em]
+ f_{\bar{x}}^j\\[0.6em]
+ f_{\bar{y}}^j\\[0.6em]
+ f_{\bar{z}}^j\\[0.6em]
+ t_{\bar{x}}^j\\[0.6em]
+ t_{\bar{y}}^j\\[0.6em]
+ t_{\bar{z}}^j\\
\end{bmatrix}
=
\begin{bmatrix}
\frac{EA}{L} & 0 & 0 & 0 & 0 & 0 & \frac{-EA}{L} & 0 & 0 & 0 & 0 & 0\\[0.5em]
- 0 & \frac{12EI_\text{z}}{L^3} & 0 & 0 & 0 & \frac{6EI_\text{z}}{L^2} & 0 & \frac{-12EI_\text{z}}{L^3} & 0 & 0 & 0 & \frac{6EI_\text{z}}{L^2}\\[0.5em]
- 0 & 0 & \frac{12EI_\text{y}}{L^3} & 0 & \frac{-6EI_\text{y}}{L^2} & 0 & 0 & 0 & \frac{-12EI_\text{y}}{L^3} & 0 & \frac{-6EI_\text{y}}{L^2} & 0\\[0.5em]
- 0 & 0 & 0 & \frac{GK_\text{v}}{L} & 0 & 0 & 0 & 0 & 0 & \frac{-GK_\text{v}}{L} & 0 & 0\\[0.5em]
- 0 & 0 & \frac{-6EI_\text{y}}{L^2} & 0 & \frac{4EI_\text{y}}{L} & 0 & 0 & 0 & \frac{6EI_\text{y}}{L^2} & 0 & \frac{2EI_\text{y}}{L} & 0\\[0.5em]
- 0 & \frac{6EI_\text{z}}{L^2} & 0 & 0 & 0 & \frac{4EI_\text{z}}{L} & 0 & \frac{-6EI_\text{z}}{L^2} & 0 & 0 & 0 & \frac{2EI_\text{z}}{L}\\[0.5em]
+ 0 & \frac{12EI_{\bar{z}}}{L^3} & 0 & 0 & 0 & \frac{6EI_{\bar{z}}}{L^2} &
+ 0 & \frac{-12EI_{\bar{z}}}{L^3} & 0 & 0 & 0 &
+ \frac{6EI_{\bar{z}}}{L^2}\\[0.5em]
+ 0 & 0 & \frac{12EI_{\bar{y}}}{L^3} & 0 & \frac{-6EI_{\bar{y}}}{L^2} & 0
+ & 0 & 0 & \frac{-12EI_{\bar{y}}}{L^3} & 0 & \frac{-6EI_{\bar{y}}}{L^2} &
+ 0\\[0.5em]
+ 0 & 0 & 0 & \frac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & \frac{-GJ}{L} & 0 &
+ 0\\[0.5em]
+ 0 & 0 & \frac{-6EI_{\bar{y}}}{L^2} & 0 & \frac{4EI_{\bar{y}}}{L} & 0 & 0
+ & 0 & \frac{6EI_{\bar{y}}}{L^2} & 0 & \frac{2EI_{\bar{y}}}{L} &
+ 0\\[0.5em]
+ 0 & \frac{6EI_{\bar{z}}}{L^2} & 0 & 0 & 0 & \frac{4EI_{\bar{z}}}{L} & 0
+ & \frac{-6EI_{\bar{z}}}{L^2} & 0 & 0 & 0 &
+ \frac{2EI_{\bar{z}}}{L}\\[0.5em]
\frac{-EA}{L} & 0 & 0 & 0 & 0 & 0 & \frac{EA}{L} & 0 & 0 & 0 & 0 & 0\\[0.5em]
- 0 & \frac{-12EI_\text{z}}{L^3} & 0 & 0 & 0 & \frac{-6EI_\text{z}}{L^2} & 0 & \frac{12EI_\text{z}}{L^3} & 0 & 0 & 0 & \frac{-6EI_\text{z}}{L^2}\\[0.5em]
- 0 & 0 & \frac{-12EI_\text{y}}{L^3} & 0 & \frac{-6EI_\text{y}}{L^2} & 0 & 0 & 0 & \frac{12EI_\text{y}}{L^3} & 0 & \frac{6EI_\text{y}}{L^2} & 0\\[0.5em]
- 0 & 0 & 0 & \frac{-GK_\text{v}}{L} & 0 & 0 & 0 & 0 & 0 & \frac{GK_\text{v}}{L} & 0 & 0\\[0.5em]
- 0 & 0 & \frac{-6EI_\text{y}}{L^2} & 0 & \frac{2EI_\text{y}}{L} & 0 & 0 & 0 & \frac{6EI_\text{y}}{L^2} & 0 & \frac{4EI_\text{y}}{L} & 0\\[0.5em]
- 0 & \frac{6EI_\text{z}}{L^2} & 0 & 0 & 0 & \frac{2EI_\text{z}}{L} & 0 & \frac{-6EI_\text{z}}{L^2} & 0 & 0 & 0 & \frac{4EI_\text{z}}{L}\\
+ 0 & \frac{-12EI_{\bar{z}}}{L^3} & 0 & 0 & 0 & \frac{-6EI_{\bar{z}}}{L^2}
+ & 0 & \frac{12EI_{\bar{z}}}{L^3} & 0 & 0 & 0 &
+ \frac{-6EI_{\bar{z}}}{L^2}\\[0.5em]
+ 0 & 0 & \frac{-12EI_{\bar{y}}}{L^3} & 0 & \frac{-6EI_{\bar{y}}}{L^2} & 0
+ & 0 & 0 & \frac{12EI_{\bar{y}}}{L^3} & 0 & \frac{6EI_{\bar{y}}}{L^2} &
+ 0\\[0.5em]
+ 0 & 0 & 0 & \frac{-GJ}{L} & 0 & 0 & 0 & 0 & 0 & \frac{GJ}{L} & 0 &
+ 0\\[0.5em]
+ 0 & 0 & \frac{-6EI_{\bar{y}}}{L^2} & 0 & \frac{2EI_{\bar{y}}}{L} & 0 & 0
+ & 0 & \frac{6EI_{\bar{y}}}{L^2} & 0 & \frac{4EI_{\bar{y}}}{L} &
+ 0\\[0.5em]
+ 0 & \frac{6EI_{\bar{z}}}{L^2} & 0 & 0 & 0 & \frac{2EI_{\bar{z}}}{L} & 0
+ & \frac{-6EI_{\bar{z}}}{L^2} & 0 & 0 & 0 & \frac{4EI_{\bar{z}}}{L}\\
\end{bmatrix}
\begin{bmatrix}
- p_\text{x}^i\\[0.6em]
- p_\text{y}^i\\[0.6em]
- p_\text{z}^i\\[0.6em]
- \Omega_\text{x}^i\\[0.6em]
- \Omega_\text{y}^i\\[0.6em]
- \Omega_\text{z}^i\\[0.6em]
- p_\text{x}^j\\[0.6em]
- p_\text{y}^j\\[0.6em]
- p_\text{z}^j\\[0.6em]
- \Omega_\text{x}^j\\[0.6em]
- \Omega_\text{y}^j\\[0.6em]
- \Omega_\text{z}^j\\
+ p_{\bar{x}}^i\\[0.6em]
+ p_{\bar{y}}^i\\[0.6em]
+ p_{\bar{z}}^i\\[0.6em]
+ \Omega_{\bar{x}}^i\\[0.6em]
+ \Omega_{\bar{y}}^i\\[0.6em]
+ \Omega_{\bar{z}}^i\\[0.6em]
+ p_{\bar{x}}^j\\[0.6em]
+ p_{\bar{y}}^j\\[0.6em]
+ p_{\bar{z}}^j\\[0.6em]
+ \Omega_{\bar{x}}^j\\[0.6em]
+ \Omega_{\bar{y}}^j\\[0.6em]
+ \Omega_{\bar{z}}^j\\
\end{bmatrix}
\end{equation}
-$E$ is Young's modulus, $G$ is the shear stiffnes, $A$ is the beam
-cross-sectional area, and $L$ is the original beam length. $I_\text{y}$ is the
-moment of inertia normal to the beam in the $\bar{y}$-direction, and
-$I_\text{z}$ is the moment of inertia normal to the beam in the
-$\bar{z}$-direction. $K_\text{v}$ is the Saint-Venant torsional stiffness.
-
-% Torsional constant:
-% https://en.wikipedia.org/wiki/Torsion_constant
-% http://mathworld.wolfram.com/TorsionalRigidity.html
-% http://physics.stackexchange.com/questions/83148/where-i-can-find-a-torsional-stiffness-table-for-different-types-of-stainless-st
-% St Venant torsion: K_v = 1/G (Austrell et al. 2004, table 3) Does it make sense?
-
-
(DIR) diff --git a/slidergrid/grid.c b/slidergrid/grid.c
t@@ -34,6 +34,35 @@ slider* create_regular_slider_grid(
return sliders;
}
+slider* create_regular_slider_grid_with_randomness(
+ const int nx,
+ const int ny,
+ const int nz,
+ const Float dx,
+ const Float dy,
+ const Float dz)
+{
+ slider* sliders;
+ sliders = malloc(sizeof(slider)*nx*ny*nz);
+
+ int i = 0; int ix, iy, iz, j;
+ for (iz = 0; iz < nz; iz++) {
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ sliders[i].pos.x = dx * ix;
+ sliders[i].pos.y = dy * iy;
+ sliders[i].pos.z = dz * iz;
+ i++;
+
+ for (j=0; j<MAX_NEIGHBORS; j++)
+ sliders[i].neighbors[j] = -1;
+ }
+ }
+ }
+
+ return sliders;
+}
+
/* Find neighboring sliders within a defined cutoff distance */
void find_and_bond_to_neighbors_n2(
slider* sliders,
(DIR) diff --git a/slidergrid/slider.h b/slidergrid/slider.h
t@@ -34,6 +34,7 @@ typedef struct {
// Macroscopic mechanical properties
Float youngs_modulus;
Float shear_modulus;
+ Float torsion_stiffness;
// inter-slider bond-parallel Kelvin-Voigt contact model parameters
Float bond_parallel_kv_stiffness; // Hookean elastic stiffness [N/m]