tadd system of equations for 3d beam - 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 897289bb9cc3233e556ccd4b863cdfaf4dea29b0
(DIR) parent 01b81b008518de59f62d4471b397364dd6b78f9a
(HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Mon, 2 May 2016 09:58:42 -0700
add system of equations for 3d beam
Diffstat:
M doc/doc.pdf | 0
M doc/doc.tex | 295 ++++++++++++++++++++++++++++---
M slidergrid/slider.h | 4 ++++
3 files changed, 275 insertions(+), 24 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@@ -1,4 +1,4 @@
-\documentclass[11pt]{article}
+\documentclass[11pt,a4paper]{article}
\usepackage{a4wide}
t@@ -14,6 +14,7 @@
%\usepackage{subfig}
%\usepackage{rotating}
\usepackage{amsmath}
+\setcounter{MaxMatrixCols}{20} % allow more than 10 matrix columns
\usepackage[T1]{fontenc} % Font encoding
\usepackage{charter} % Serif body font
t@@ -37,7 +38,7 @@ maxcitenames=2, backend=bibtex8]{biblatex}
\begin{document}
\title{Lagrangian model of the elastic, viscous and plastic deformation of a
- series of bonded points moving on a frictional surface}
+ series of bonded nodes moving on a frictional surface}
\author{Anders Damsgaard}
\date{{\small Institute of Geophysics and Planetary Physics\\Scripps Institution
t@@ -46,54 +47,300 @@ of Oceanography\\University of California, San Diego}\\[3mm] Last revision:
\maketitle
-
\section{Methods}
-The Lagrangian points are connected with visco-elastic beams which are resistive
-to relative translational or rotational movement between a pair of bonded
-points. 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
-points. 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.
+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
+($\boldsymbol{a}$) are found from the sums of forces:
+\begin{equation}
+ \boldsymbol{a}_i =
+ \frac{
+ \boldsymbol{f}_i^\text{d}
+ + \boldsymbol{f}_i^\text{f}
+ + \sum^{N_c}_j \left[
+ \boldsymbol{f}_{i,j}^\text{p} +
+ \boldsymbol{f}_{i,j}^\text{s}
+ \right]
+}{m_i}
+ + \boldsymbol{g}
+\label{eq:n2-tran}
+\end{equation}
+where $\boldsymbol{f}_i^\text{d}$ is the gravitational driving stress due to
+surface slope, $\boldsymbol{g}$ is the gravitational acceleration, and
+$\boldsymbol{f}_i^\text{f}$ is the frictional force provided if the point is
+resting on the lower surface. Bonded interaction with another point $j$
+contributes to translational acceleration through bond-parallel and bond-normal
+shear forces, $\boldsymbol{f}_{i,j}^\text{p}$ and
+$\boldsymbol{f}_{i,j}^\text{s}$, respectively.
+
+The angular accelerations ($\boldsymbol{\alpha}$) are found from the sums of
+torques:
+\begin{equation}
+ \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}}
+ \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.
+
+
+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.
The deformation and reactive forces are determined relative to the orientation
of the bond. Common geometrical vectors include the inter-distance vector
-$\boldsymbol{d}$ between points $\boldsymbol{p}_i$ and $\boldsymbol{p}_j$:
+$\boldsymbol{d}$ between nodes $\boldsymbol{p}_i$ and $\boldsymbol{p}_j$:
\begin{equation}
\boldsymbol{d}_{i,j} = \boldsymbol{p}_i - \boldsymbol{p}_j
\end{equation}
-which in normalized form constitutes the bond-parallel normal vector:
+which in normalized form constitutes the bond-parallel unit vector:
\begin{equation}
\boldsymbol{n}_{i,j} = \frac{\boldsymbol{d}_{i,j}}{||\boldsymbol{d}_{i,j}||}
\end{equation}
-The points are moving by translational and rotational velocities. The combined
-relative velocity between the points is found as \citep{Hinrichsen2004,
+The nodes move by translational and rotational velocities. The combined
+relative velocity between the nodes is found as \citep[e.g.][]{Hinrichsen2004,
Luding2008}:
\begin{equation}
\boldsymbol{v}_{i,j} = \boldsymbol{v}_i - \boldsymbol{v}_j +
- \frac{d_{i,j}}{2} \times \omega_i +
- \frac{d_{i,j}}{2} \times \omega_j
+ \frac{\boldsymbol{d}_{i,j}}{2} \times \boldsymbol{\omega}_i +
+ \frac{\boldsymbol{d}_{i,j}}{2} \times \boldsymbol{\omega}_j
+\end{equation}
+The velocity can be decomposed into spatial components relative to the bond
+orientation, e.g.\ the bond-parallel and bond-shear velocity, respectively:
+\begin{equation}
+ v^\text{p}_{i,j} = \boldsymbol{v}_{i,j} \cdot \boldsymbol{n}_{i,j}
+\end{equation}
+\begin{equation}
+ \boldsymbol{v}^\text{s}_{i,j} = \boldsymbol{v}_{i,j} - \boldsymbol{n}_{i,j}
+ \left(
+ \boldsymbol{v}_{i,j}
+ \cdot
+ \boldsymbol{n}_{i,j}
+ \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{Bond tension and compression}
+\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}
-\subsection{Bond shear}
-\subsection{Bond twist}
+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{Bond bend}
+\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{Temporal integration}
+\subsection{Twisting resistance}
+\subsection{Bending resistance}
+\subsection{Temporal integration}
+Once the force and torque sum components at time $t$ have been determined, the
+kinematic degrees of freedom at time $t+\Delta t$ can be found by explicit
+temporal integration of moment balance equations~\ref{eq:n2-tran}
+and~\ref{eq:n2-ang}.
+We use an integration scheme based on the third-order Taylor expansion, which
+results in a truncation error on the order of $O(\Delta t^4)$ for positions and
+$O(\Delta t^3)$ for velocities. This scheme includes changes in acceleration as
+the highest order term, which are approximated by backwards differences. For
+the translational degrees of freedom:
+\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 +
+ \frac{1}{6} \frac{\boldsymbol{a}^i_{t}
+ - \boldsymbol{a}^i_{t - \Delta t}}{\Delta t} \Delta t^3
+\end{equation}
+\begin{equation}
+ \boldsymbol{v}^i_{t+\Delta t} =
+ \boldsymbol{v}^i_{t} +
+ \boldsymbol{a}^i_{t} \Delta t +
+ \frac{1}{2} \frac{\boldsymbol{a}^i_{t}
+ - \boldsymbol{a}^i_{t - \Delta t}}{\Delta t} \Delta t^2
+\end{equation}
+At $t=0$ the acceleration change term is defined as zero. The angular degrees
+of freedom are found correspondingly:
+\begin{equation}
+ \boldsymbol{\Omega}^i_{t+\Delta t} =
+ \boldsymbol{\Omega}^i_{t} +
+ \boldsymbol{\omega}^i_{t} \Delta t +
+ \frac{1}{2} \boldsymbol{\alpha}^i_{t} \Delta t^2 +
+ \frac{1}{6} \frac{\boldsymbol{\alpha}^i_{t}
+ - \boldsymbol{\alpha}^i_{t - \Delta t}}{\Delta t} \Delta t^3
+\end{equation}
+\begin{equation}
+ \boldsymbol{\omega}^i_{t+\Delta t} =
+ \boldsymbol{\omega}^i_{t} +
+ \boldsymbol{\alpha}^i_{t} \Delta t +
+ \frac{1}{2} \frac{\boldsymbol{\alpha}^i_{t}
+ - \boldsymbol{\alpha}^i_{t - \Delta t}}{\Delta t} \Delta t^2
+\end{equation}
+The numerical time step $\Delta t$ is found by considering the largest elastic
+stiffness in the system relative to the smallest mass:
+\begin{equation}
+ \Delta t =
+ \epsilon
+ \left[
+ \min (m_i)^{-1}
+ \max \left(
+ \max \left(
+ \frac{E A_{i,j}}{||\boldsymbol{d}_{0}^{i,j}||}
+ \right)
+ ,
+ \max \left(
+ \frac{G A_{i,j}}{||\boldsymbol{d}^{i,j}||}
+ \right)
+ \right)
+ \right]^{-1/2}
+\end{equation}
+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}:
+\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\\
+ \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]
+ \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}\\
+ \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\\
+ \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/slider.h b/slidergrid/slider.h
t@@ -31,6 +31,10 @@ typedef struct {
// moment of inertia [kg m*m]
Float moment_of_inertia;
+ // Macroscopic mechanical properties
+ Float youngs_modulus;
+ Float shear_modulus;
+
// inter-slider bond-parallel Kelvin-Voigt contact model parameters
Float bond_parallel_kv_stiffness; // Hookean elastic stiffness [N/m]
Float bond_parallel_kv_viscosity; // viscosity [N/(m*s)]