structureModel.md
-
Amira Abdel-Rahman authoredAmira Abdel-Rahman authored
Structural Mechanics Models
Bar Model
\begin{bmatrix}F_{n} \end{bmatrix}=\begin{bmatrix} K \end{bmatrix} \begin{bmatrix}D_{n} \end{bmatrix} \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_rF_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_rM_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)".