# Validation of MetaVoxels Dynamics
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!
