Structural Mechanics Models

Bar Model

$\begin{bmatrix}F_{n} \end{bmatrix}=\begin{bmatrix} K \end{bmatrix} \begin{bmatrix}D_{n} \end{bmatrix}$

Each element (bar) has 6 degrees of freedom:

$\begin{vmatrix}F_{x1} \\\ F_{y1} \\\ F_{z1} \\\ F_{x2} \\\ F_{y2} \\\ F_{z2} \end{vmatrix}=\begin{bmatrix} K \end{bmatrix} \begin{vmatrix} X_{1} \\\ Y_{1} \\\ Z_{1} \\\ X_{2} \\\ Y_{2} \\\ Z_{2} \end{vmatrix}$

(2D Version of K matrix)

$\begin{bmatrix} K \end{bmatrix}= \frac{ E_c I}{l^3} \begin{vmatrix} 12 & 6l & -12& 6l \\\ & 4l^2 & -6l& 2l^2\\\ & & 12 & -6l\\\ & & & 4l^2 \end{vmatrix}$

Beam Model

$\begin{bmatrix}F_{n} \\\ M_{n} \end{bmatrix}=\begin{bmatrix} K \end{bmatrix} \begin{bmatrix}D_{n} \\\ \theta_{n} \end{bmatrix}$

Each element (bar) has 12 degrees of freedom:

$\begin{vmatrix}F_{x1} \\\ F_{y1} \\\ F_{z1} \\\ M_{x1} \\\ M_{y1}\\\ M_{z1} \\\ F_{x2} \\\ F_{y2} \\\ F_{z2} \\\ M_{x2} \\\ M_{y2}\\\ M_{z2} \end{vmatrix}=\begin{bmatrix} K \end{bmatrix} \begin{vmatrix} X_{1} \\\ Y_{1} \\\ Z_{1} \\\ \theta_{x1} \\\ \theta_{y1} \\\ \theta_{z1} \\\ X_{2} \\\ Y_{2} \\\ Z_{2} \\\ \theta_{x2} \\\ \theta_{y2} \\\ \theta_{z2} \end{vmatrix}$

$\begin{bmatrix} K \end{bmatrix}= \begin{vmatrix} a_1 & 0 & 0 & 0 & 0 & 0 & -a_1 & 0 & 0 & 0 & 0 & 0 \\\ & b_1 & 0 & 0 & 0 & b_2 & 0 & -b_1 & 0 & 0 & 0 & b_2 \\\ & & b_1 & 0 & -b_2 & 0 & 0 & 0 & -b_1 & 0 & -b_2 & 0 \\\ & & & a_2 & 0 & 0 & 0 & 0 & 0 & -a_2 & 0 & 0 \\\ & & & & 2b_3 & 0 & 0 & 0 & b_2 & 0 & b_3 & 0 \\\ & & & & & 2b_3 & 0 & -b_2 & 0 & 0 & 0 & b_3 \\\ & & & & & & a_1 & 0 & 0 & 0 & 0 & 0 \\\ & & & & & & & b_1 & 0 & 0 & 0 & -b_2 \\\ & & & & & & & & b_1 & 0 & b_2 & 0 \\\ & & & & & & & & & a_2 & 0 & 0 \\\ & (sym) & & & & & & & & & 2b_3 & 0 \\\ & & & & & & & & & & & 2b_3 \end{vmatrix}$

$a_1= \frac{E_c A}{l}$

$a_2= \frac{G_c J}{l}$

$b_1= \frac{12 E_c I}{l^3}$

$b_2= \frac{6 E_c I}{l^2}$

$b_3= \frac{2 E_c I}{l}$

$Bending\ Moment\ of\ Inertia \ I= \frac{bh^3}{12}$

$Polar \ Moment\ of\ Inertia\ J= \frac{bh(bb+hh)}{12}$

$Modulus \ of\ Rigidity\ G= \frac{E}{2(1+nu)}$

Dynamic Parallel Model

Based on (1) and (2), if at each time step I consider the first node of each element to located be on the primary X axis, $X_{1} , Y_{1} , Z_{1} , \theta_{x1} , \theta_{y1} , \theta_{z1}$ will all be zeros and the transformation into the X-direction stiffness matrix is a trivial calculation. We then transform $X_{2} , Y_{2} , Z_{2} , \theta_{x2} , \theta_{y2} , \theta_{z2}$ to the original rotation matrix. Therefore, the larger matrix calculation is reduced to the following equations:

$F_{x1}= -a_1 X_2$

$F_{y1}= -b_1 Y_2+b_2\theta_{z2}$

$F_{z1}= -b_1 Z_2+b_2\theta_{y2}$

$M_{x1}= -a_2 \theta_{x2}$

$M_{y1}= b_2 Z_2+b_3\theta_{y2}$

$M_{z1}= -b_2 Y_2+b_3\theta_{z2}$

$F_{x2}= -F_{x1}$

$F_{y2}= -F_{y1}$

$F_{z2}= -F_{z1}$

$M_{x2}= a_2 \theta_{x2}$

$M_{y2}= b_2 Z_2+ 2b_3\theta_{y2}$

$M_{z2}= -b_2 Y_2+2b_3\theta_{z2}$

Then the forces on each node is calculated as follows:

$F_t= \sum_{b=1}^{b=n} \vec{F_b}$

$M_t= \sum_{b=1}^{b=n} \vec{M_b}$

Damping

"damping must be included at the local interaction between voxels, not just applying a force proportional to each voxel’s global velocity which would also damp rigid body motion. The local damping between adjacent voxels ensures that modal resonances at the scale of a single voxel do not accumulate.".

For each edge, first the average position, velocity, and angular velocity are calculated. The velocity of the second voxel relative to the first:

$\vec{V}_{2 \rightarrow 1} = (\vec{V}_2-\vec{V}_a) + (\vec{D}_2-\vec{D}_a) * \vec{\omega}_a$

$\vec{V}_2$ is the velocity of the second voxel, $\vec{V}_a$ is the average velocity of two voxels, $\vec{D}_2$ is the position of the second voxel, $\vec{D}_a$is the average position, $\vec{\omega}_a$ is the average angular velocity.

For each voxel we calculate the damping force as

$F_d=2 \xi \sqrt{m k} V_r$

$F_d$ is the damping force, $m$ is the mass attached to a edge with stiffness $k$ and a relative velocity $V_r$. The damping ratio $\xi$ is normally selected to be 1, corresponding to critical damping.

The angular velocities are calculated as follows:

$M_d=2 \xi \sqrt{I k_\phi} \omega_r$

$M_d$ is the damping force with rotational inertia $I$, with edge of stiffness $k$ and relative angular velocity $\omega_r$, rotational angular ratio $\xi$ is also 1.

Timestep calculation

$dt < \frac{1}{2 \pi \omega_{0_{m}}}$

$\omega_{0_{m}}$ is the maximum natural frequency of any bond, which is calculated as follows:

$\omega_{0_{max}}=\sqrt{\frac{k_b}{m_m}}$

$k_b$ is stiffness of the bond, $m_m$ is the minimum of either mass connected to the bond.

Collision and Friction

The effective normal stiffness on the floor is limited by the maximum stiffness between any connected masses. Therefore, according to (1) I set the stiffness of each voxel contacting the floor as the the stiffness of the floor in that location. Although this allows significant floor penetration in some cases, the qualitative behavior is appropriate.

A Couloumb friction model is used, even though a standard linear model is would provide a relatively realistic simulation.

The voxel will resist motion until:

$|F_l|>\mu_s f_n$

$F_l$ is the horizontal force parallel to the ground, $\mu_s$ is the coefficient of static friction between the voxel and ground, and $F_n$ is the normal force pressing te voxel into the plane of the ground.

"A boolean flag is set indicating to the simulation that this voxel should not move laterally, but can still move in the direction normal to ground such that it can be unweighted and then moved laterally".

"Once the static friction threshold has been exceeded at any given time step, the voxel is allowed to begin motion in the appropriate lateral direction by clearing the boolean static friction flag". The voxel is allowed to move in any direction, but a friction force is applied in the opposing lateral direction

$|F_l|=\mu_d f_n$

$\mu_d$ is the dynamic coefficient of friction.

"In order to properly detect when a voxel has stopped lateral motion, a minimum motion threshold must be set".

In order to detect a stopping voxel, the voxel is halted when:

$V_l\leq \frac {f_n \mu_d dt}{m}$

"Collisions are also damped normal to the direction of contact with a user variable damping ratio ranging from zero (no damping) to 1 (critical damping)".