# Validation of MetaVoxels Dynamics
What raised the question?
I put together a tendon model that calculated reactions at each node in the tendon path including friction in the hopes of recreating the experimental data that alfonso acquired. The model has just one parameter, coefficient of friction, and captures the behavior of varying deflection along the beam. But the dynamic behavior seemed very off, the beam doesn't make it to equilibrium and there is a very slow compression wave that travels down the beam under initial loading.
<img src="json/tutorial/output.mp4" /><br></br>
Here's a running notebook of my work to perform validation of the dynamic characteristics of the MetaVoxels beam solver. The key questions are:
1. Are we calculating the mass matrix correctly?
- The issue here probably crops up from the implementation of voxcad routines where nodes are the "1st class" object and beams are the abstraction, vs. a true slender beam solver where the nodes are an abstraction and the beams are physical
- Lumped vs distributed approximations?
- Probably have to excavate from Frame3DD :-(
- Are we ignoring important couplings?
2. How are we calculating our various forms of damping?
- There is both local damping of bonds and global damping of motion to calculate and validate
3. Accuracy of the solution over time
- Current integrator is a symplectic euler scheme, should preserve total energy (in absence of damping)
- Can we reframe the code to leverage DifferentialEquations.jl library?
## The proposed model system
As a first pass for validation let's use a cantilevered beam made of a line of voxels. Our one hard constraint on beam parameters is that the length/height or width ratio be > 10, so lets use 100. Properties of the beam are as follows:
- L = 1 m
- w = h = 0.01 m
- A = 0.001 m^2
- Iyy, Izz = 0.01^4/12 m^4
- Izz = 0.01^4/6 m^4
- E = 6.9e10 N/m^2
- nu = 0.33
- rho = 2700 kg/m^3
- m = 2.7 kg/m
We'll discretize the beam into n segments from 5 to 10 to 20. To test the beam dynamics we'll apply a series of initial unit impulses to the beam tip using the externalforces! method. The impulses will be a transverse impulse to excite beam bending, an axial impulse to excite the axial modes, and a torsional impulse about the x-axis.
## How will we do the comparison?
There are analytical solutions to these beam loadings which we can compare to, but we must decide where we want to spend our effort. One option is to calculate the exact time solution to these unit impulsive loadings and compare pointwise through time. The other option, which may be quicker, is to extract the natural frequencies and mode shapes from the metavoxels response and comparing to analytical frequencies and mode shapes. Because metavoxels is matrix-free we can't solve the eigenvalue problem directly.
### Predicted natural frequencies
Bending : wn = an^2 \* sqrt(EI/(mL^4)) : an^2 = 3.516, 22.035, 61.697 ...
Axial : wn = n\*pi/2 \* sqrt(EA/(mL^2))
Torsional : wn = n\*pi/2 \* sqrt(GJ/(IzzL^2))
### Measurements
*1st Bending*
| n | Exact | Solidworks Beam | MetaVoxels (Old)| MetaVoxels (Chris) |
| :--: | :--: | :--------: | :---: | :---: |
| 5 | 8.17 Hz | 8.02 Hz | 10.6 Hz | - |
| 10 | - | 8.13 Hz | 8.52 Hz | - |
| 20 | - | 8.16 Hz | 6.37 Hz | 8.3 Hz* |
\* Resolution limited
*2nd Bending*
| n | Exact | Solidworks Beam | MetaVoxels (Old)| MetaVoxels (Chris) |
| :--: | :--: | :--------: | :---: | :---: |
| 5 | 51.2 Hz | 48.1 Hz | 68.9 Hz | :---: |
| 10 | - | 50.4 Hz | 52.4 Hz | :---: |
| 20 | - | 51.0 Hz | 38.6 Hz | 50.9 Hz |
*3rd Bending*
| n | Exact | Solidworks Beam | MetaVoxels (Old)| MetaVoxels (Chris) |
| :--: | :--: | :--------: | :---: | :---: |
| 5 | 143 Hz | 130 Hz | 194 Hz | :---: |
| 10 | - | 140 Hz | 147 Hz | :---: |
| 20 | - | 142 Hz | 109 Hz | 142 Hz |
*1st Axial*
| n | Exact | Solidworks Beam | MetaVoxels|
| :--: | :--: | :--------: | :---: |
| 5 | 1264 Hz | 1259 Hz | |
| 10 | - | 1263 Hz | |
| 20 | - | 1264 Hz | |
## What's Wrong?
The issue is how the lumped mass is being calculated. Solidworks is a good reference for us because it uses a lumped model (indicated by frequencies converging from below). Metavoxels is failing because the lumped masses are based solely on nodal cross-section properties and not the edges which comprise the actual structure. This introduces 2 errors:
1. Lumped mass doesn't change with number of elements (length of edges) in a beam
The code for calculating the lumped mass in run.jl is :
```Python
nomSize = round(sqrt(node["material"]["area"]);digits=10)
nomSize = nomSize * 2.0
mass = round(nomSize * nomSize * nomSize * rho;digits=10)
```
This explains why adding nodes decreases the natural frequencies so much, the total mass of the beam is increasing as we add nodes.
2. The nodes at beam ends are not treated appropriately
The first and last element of the beam should have half the lumped mass as the interior elements because they have only one edge attached instead of two. MetaVoxels doesn't do this currently. We can confirm that with the code and by comparing our MetaVoxels natural frequencies with the expected values. For n=5 elements the metavoxels lumped mass was 0.0216 kg/node while the correct lumped mass was 0.054 kg/interior node. The natural frequency should be proportional to mass^(1/2) and we would expect the metavoxels frequencies to be ~1.6x higher than theory but due to the end mass they were only ~1.3x higher.
## The Fix
I've written a silly first pass fix but I think we can still do better. Here's my fix outlined with the rest of the setup code
1. Setup nodes using the input file (including nodal properties like mass!)
2. Setup edges using the input file AND the nodal material properties (i.e. averaging E, rho etc. between endpoints)
3. *Now that we have edge properties recalculate nodal mass and reinitialized nodal material properties because the struct is immutable*
This could be cleaner if we can delay calculating and initializing nodal mass properties until after the edges are completed. But hey look the bending frequencies match now!
## But wait, there's more!
We still haven't adressed calculating moments of interia. For the beam above with 100:1 slenderness ratio and low frequency bending this doesn't really matter. And actually as the beams get stockier we have bigger issues with accounting for shear effects ala timoshenko formulations (see https://www.sciencedirect.com/science/article/pii/S0895717708002112). Regardless, we should try to get accurate moments of inertia. This means, similar to the masses looping over all edges attached to a node and calculating their contribution to the rotational inertia of said node. The process looks like (python pseudocode)
```python
for node in nodes:
for edge in node.edges:
inertia0 = inertia of 1/2 edge about node coord in edge coords
node.inertia += transform_edge_to_global(inertia0)
node.inertia_inverse = inverse(node.inertia)
```
Now this gets messy. The inertia0 term is aligned with the edge coordinates and will return a diagonal array of moments of inertia. When we transform to global coordinates we will (in all but a few cases) get products of inertia or off-diagonal terms in the inertia matrix. To be consistent we should keep these terms and calculate the inverse inertial matrix. Update the nodal rotational velocity then becomes a matrix-vector product of the applied moment and the inertia_inverse matrix. Additionally if we're allowing for large deformations we'll need to recalculate node.inertia and node.inertia_inverse in global coordinates (unless I misunderstand and we do our moment*inertia_inverse in nodal coords before transforming to global coords).
## Validation against voxel beam
I've begun work collecting the impulse response of a cantilevered 5x1 voxel beam. The data shown below was acquired using the microcontroller-imu setup from the voxelcopter, but the next step is setting up the laser displacement sensor for higher rate, non-contact measurement.
<img src="images/voxel_dynamic.png" width=500><br></br>
