Multibody Dynamics with Finite Element Simulation of Simple Pendulum.

1.1. Problem description.

A simple pendulum under the action of uniform gravity is simulated using Multibody Dynamics and Finite Element Method.
The objective is to determine the stresses induced in the material, resulting from the motion dynamics (force and torque) generated at the hinge joint.

The next two videos show the final results of the simulation. First, the stress on the pendulum, second, the stress on the pin.

1.2. Software.

Pre-post processor: FreeCAD + MotionSim and MotionView Workbenches.

MbD solver: MBDyn.

FEM solver: Calculix.

Mesher: Gmsh.

Video editor: FFmpeg.

1.3. Pre-processing: CAD.

Two FreeCAD App::Parts were created:

Pin
Pendulum

Each part contains one PartDesign::Body.
Each body contains one Part::TopoShape.

Isometric view of the CAD geometries and the global coordinate system.

1.4. Pre-processing: MbD.

Material properties of steel.
Name:	Density [kg/m^3]:	Poisson's ratio []:	Young's modulus [Pa]:
Steel	7850	0.25	200e9

Pin physical properties.
Part name:	Mass [kg]:	Center of mass matrix of inertia [kg*m^2]:
Pin	0.148	4.809e-05 0.0 0.0 0.0 4.809e-5 0.0 0.0 0.0 7.398e-6

Pendulum physical properties.
Part name:	Mass [kg]:	Center of mass matrix of inertia [kg*m^2]:
Pendulum	0.359	3.138e-05 0.0 0.0 0.0 7.548e-4 0.0 0.0 0.0 7.801e-4

The pin is grounded (all 6 degrees of freedom are fixed to the global inertial coordinate system).
The pin and pendulum each have one joint marker attached, as shown below.

One joint marker is attached to each part.

A hinge joint is defined between the two markers.

The hinge joint ensures that:

The marker origins remain coincident.
The marker Z axes remain parallel.
The initial angle (assembly angle) between the marker X axes is zero.

Left = exploded assembly. Right = solved assembly.

Earth gravity (9.807 m/s^2) is pointing in global negative Y direction.

The pre-processor generates an input file for the MbD solver, in *.mbd format.

1.5. Simulation: MbD.

The MbD solver is launched.
The results of the MbD simulation are stored as properties of the FreeCAD parts and joints, and saved in the *.FCStd file.
For each part, the following data are stored:

Exploded and assembled position and orientation.
Position and orientation time series.
Velocity and angular velocity time series.
Acceleration and angular acceleration time series.

For each joint, the following are stored:

Joint reaction force time series.
Joint reaction torque time series.

1.6. Post-processing: MbD.

The MotionView post-processor allows user to:

Inspect the simulation results data via plotting, using Matplotlib.
Export the simulation results in *.csv format.
Inspect the simulation results using the FreeCAD Spreadsheet Workbench.
Animate the motion of the mechanism within the 3D view.
Export each frame of the animation as a *.png image, and create a smooth video from the images.

The next video shows the pendulum motion.

In the next video:

The free body diagram of the pendulum is animated.
The pendulum weight vector is shown at its center of mass.
The D'Alembert force and torque vectors are shown at the center of mass of the pendulum.
The joint force vector is shown at center of the joint hole.
Force vectors are depicted by red lines. Torque by blue lines.

1.7. Pre-processing: FEM.

The FEM simulation is applied to the Pendulum.
FreeCAD meshes the pendulum using Gmsh.

Two FEM "bearing loads" are applied to the mesh at the hole of the hinge joint.
The first is from the hinge joint reaction force.
The second is from the hinge joint reaction torque. In this case, joint torque is zero.
These loads are automatically distributed among adequate nodes on the joint hole surface.
The D'Alembert Principle ensures the pendulum FEM model is in static equilibrium. For this, the total acceleration vector at the center of mass of each tetrahedron mesh element is computed, and the negative of this acceleration is applied to each element in the FEM model.
To prevent rigid-body motion, an algorithm chooses three suitable nodes and applies fixed boundary conditions using the "3-2-1" approach.

In the next video:

The joint load distribution among the FEM nodes is shown.
Gravity and D'Alembert acceleration vectors are shown acting at the center of mass of an evenly distributed sample of tetrahedron elements.
Load vectors are depicted by red lines. Acceleration vectors in orange.

1.8. Simulation: FEM of the pendulum.

The pre-processor generates a set of input files for the FEM solver, one for each time step of the MbD simulation.
The FEM input files are in *.inp format.
The FEM solver is launched.
The results of each FEM simulation are stored as properties of the FreeCAD mesh object, and are saved within the *.FCStd file. These results include:

Nodal stress matrices time series.
321-constraint reaction forces time series.

1.9. Post-preocessing: FEM.

The MotionView post-processor allows the user to:

Compute von-Mises stresses from the stress matrices.
Compute Max Principal stresses from the stress matrices.
Animate the FEM results using contour plots.
Closely examine stresses by pointing to the FEM mesh nodes.
Export the frames as images, and create fluid animations.

The next video shows the results of the MbDFEM simulation. The contour plots represent von-Mises stresses.

1.10. Simulation validation.

One way to validate the simulation is to compare the reaction forces that take place at the nodes used for 321-constraint; with the joint loads applied to the FEM nodes. If the computation of the D'Alembert accelerations is correct, then the pendulum must remain in static equilibrium through all the time-steps. Then, the 321-constraint nodes must see no reaction force. In reality, due to error introduced during numeric approximations, the reaction forces at the 321-constraint nodes must be a small fraction of the joint loads. The next plots show both the joint force (the load applied to the FEM nodes), and the forces that take place at the 321-constraint nodes.

One can see that the 321-constraint forces are three orders of magnitude less than the applied joint loads.

A second way to verify the simulation is to inspect the stress around the nodes used for 321-constraint. If the 321-constraint forces were a significant fraction of the applied loads, then the FEM solver would return a significant stress at the 321-constraint nodes. However, the stress at 321-constraint nodes is negligible compared to the stress produced by the joint load, as shown in the contour plots.

A third way to verify the simulation is to perform an equivalent simulation using another software. The pendulum was simulated using Ansys Mechanical. The images below show that both programs produce equivalent stress distribution, and that the maximum stresses occur at the same points and present negligible difference.

Max von-Mises stress from Ansys Mechanical.

Max von-Mises stress obtained from FreeCAD MotionSim Workbench.

1.11. Simulation: FEM of the pin.

In this case:

Two sets of the pin FEM nodes are grounded.
No 3-2-1 constraint is required.
The pin does not undergo linear or angular acceleration.
D'Alembert force and torque are zero.
The joint force is distributed among the adequate FEM nodes.

In order to have a face where the joint force is applied, the cylindrical face of the pin CAD geometry was split as shown below:

Pin cylindrical face split in three faces. Joint force is applied to the nodes belonging to the middle face.

All the FEM nodes belonging to the two circular faces of the pin were grounded, as shown below:

All FEM nodes belonging to the two circular faces are grounded.

The next video depicts the joint force acting on the pin:

The next video shows the joint force load distribution among the FEM nodes:

The next video shows the von-Mises stress on the pin:

1.12. FEM load from joint torque.

In the previous example, the direction of the hinge joint Z axis (the rotation axis), was set perpendicular to the gravity vector. Thus, the hinge joint torque was zero, and there was no FEM load derived from joint torque. In this example, the pin is rotated 45 degrees about global X, as shown below:

Oblique pin exploded and solved assembly.

This change in the assembly means that:

The joint reaction force becomes oblique.
The joint reaction torque is no longer zero.
The are two FEM loads derived from joint force.
The FEM loads derived from joint torque are no longer zero.

1.12.1 Joint force.

The joint-force vector is resolved into its side-force and axial-force components.
The axial force is parallel to the hinge Z axis.
The next relation is satisfied:

joint-force = side-force + axial-force.

The next video shows the joint-force vector (continous red line), as well as its side-force and axial-force components (dashed red lines). The body weight vector is also shown, acting at the center of mass of the pendulum.

The next video shows the distribution of the joint side-force component among the FEM nodes, along with gravity and d'Alembert acceleration vectors acting at the center of mass of a representative sample of tetrahedron volume elements.

The next video shows the distribution of the joint axial-force component among the FEM nodes, along with gravity and d'Alembert acceleration vectors acting at the center of mass of a representative sample of tetrahedron volume elements.