## 2 Material and Methods
### Analysis from a Mechanical Engineering Standpoint
#### 2.2.5 Discretization  
To perform the analysis, it is essential to discretize the model, which involves meshing. This procedure entails dividing the model into multiple elements defined by nodes. The analysis variables will be calculated at these nodes, and values will be interpolated between them to determine the value within each element. The accuracy of this interpolation depends on the number of nodes within an element; a higher number of nodes (i.e., smaller element size) results in a more precise solution. Following this brief explanation of FEM and its importance in addressing the problem, the automatic discretization provided by the software can be observed (Figure 7).  
Beginning with the mathematical foundation underlying discretization, the equilibrium equations for a continuous system can be formulated as partial differential equations. This equation describes the behavior of the model, specifically for a linear static analysis:  
\[\
abla\sigma+f=0 \tag{1}\]  
where \(\sigma\) is the stress tensor and \(f\) is the applied force vector [36].  
In the finite element method, the principle of virtual work is used:  
\[\int_{\Omega}\delta\epsilon:\sigma d\Omega=\int_{\Omega}\delta u\cdot fd\Omega+ \int_{\mathbb{T}_{t}}\delta u\cdot td\Gamma \tag{2}\]  
where \(\delta\epsilon\) is the virtual variation of strain, \(\delta\)u is the virtual variation of displacement, \(\Omega\) is the body domain, and \(\Gamma_{t}\) is the boundary with applied traction \(t\)[36].  
Now, continuing with the discretization process itself, the continuous domain \(\Omega\) is divided into a finite number of elements e, such that:  
\[\Omega\approx\bigcup_{\epsilon=1}^{n_{e}}\Omega_{\epsilon} \tag{3}\]  
Each \(\Omega_{\epsilon}\) is a subregion with simple geometry; in this study, it will be tetrahedral [37].  
Figure 6: Brick base fixing (**left**) and flywheel fixing (**right**).  
Within each element, displacements \(u\), are approximated using shape functions \(N_{i}\):  
\[u^{e}(x)=\sum_{i=1}^{n_{en}}N_{i}(x)u_{i} \tag{4}\]  
where \(u_{i}\) are the nodal displacements of the element and \(n_{en}\) is the number of nodes per element [36].  
The relationship between nodal forces and nodal displacements is expressed through the elemental stiffness matrix \(K^{e}\):  
\[K^{e}=\int_{\Omega_{e}}(B^{e})^{T}DB^{e}d\Omega \tag{5}\]  
where \(B^{e}\) is the strain displacement matrix and \(D\) is the constitutive matrix depending on the material [36]. This elemental stiffness matrix \(K^{e}\) is assembled into a global stiffness matrix \(K\):  
\[K=\sum_{e=1}^{n_{e}}K^{e} \tag{6}\]  
Similarly, the nodal force vectors are assembled to obtain the global force vector \(F\). The discrete problem is reduced to solving a system of linear equations:  
\[F=Ku \tag{7}\]  
where \(u\) is the vector of nodal displacement.  
Finally, the system can be solved using numerical methods such as Gaussian elimination, iterative methods (e.g., conjugate gradient), among others [36].  
As mentioned previously, in this study, quadratic interpolation elements, specifically parabolic (tetrahedral) elements, were used as they better conform to the model, particularly in curved regions, to achieve results that are as accurate as possible. The software automatically generates a coarser mesh in larger areas and a finer mesh in smaller areas,  
Figure 7: Automatic mesh of the model.  
allowing the grid to better fit the geometry of the part and improve result accuracy. For this investigation, a configuration was defined with an element size of 0.005 m, a maximum element growth ratio of 1.1, a maximum turn angle of 65\({}^{\circ}\), and a minimum element size of 0.001 m, providing greater control over the mesh. Additionally, the mesh was refined in specific areas where problems are anticipated to conduct a more precise study. This study utilized a mesh comprising 6,530,981 elements and 9,960,107 nodes.