A Nonlinear CFD/Multibody Incremental-Dynamic Model for A Constrained Mechanism

: Numerical analysis of a multibody mechanism moving in the air is a complicated problem in computational ﬂuid dynamics (CFD). Analyzing the motion of a multibody mechanism in a commercial CFD software, i.e., ANSYS Fluent ® , is a challenging issue. This is because the components of a mechanism have to be constrained next to each other during the movement in the air to have a reliable numerical aerodynamics simulation. However, such constraints cannot be numerically mod-eled in a commercial CFD software, and needs to be separately incorporated into models through the programming environment, such as user-deﬁned functions (UDF). This study proposes a nonlinear-incremental dynamic CFD/multibody method to simulate constrained multibody mechanisms in the air using UDF of ANSYS Fluent ® . To testify the accuracy of the proposed method, Newton–Euler dynamic equations for a two-link mechanism are solved using Matlab ® ordinary differential equations (ODEs), and the numerical results for the constrained mechanisms are compared. The UDF results of ANSYS Fluent ® shows good agreement with Matlab ® , and can be applied to constrained multibody mechanisms moving in the air. The proposed UDF of ANSYS Fluent ® calculates the aerodynamic forces of a ﬂying multibody mechanism in the air for a low simulation cost than the constraint force equation (CFE) method. The results could have implications in designing and analyzing ﬂying robots to help human rescue teams, and nonlinear dynamic analyses of the aerodynamic forces applying on a moving object in the air, such as airplanes, birds, ﬂies, etc.


Introduction
Analyzing the equations of motion for a multibody mechanism is an intricate issue in the field of dynamic mechanics, since several dynamic equations with high-order nonlinearities are essentially required to be solved [1,2]. If a multibody mechanism enlarges with some active and passive joints moving in the air, the difficulties of analysis of the system will be amplified due to the initial and boundary conditions of the model [3]. Such mechanisms can drift away from their joint connections during their movement in the air owing to the set of complicated aerodynamic forces, such as drag, shear, and normal [4].
Constraint force equation (CFE) has been widely used for dynamic analysis of multibody mechanisms [5][6][7][8][9][10][11]. In CFE, the constraint forces and momentums are calculated using multibody dynamic equations, and then incorporated into the model as a separateexternal force on each part/component of a multibody mechanism. Although CFE is in wide using for various types of multibody mechanisms under different initial and boundary conditions, the main disadvantage of that is the external forces applying on a mechanism should be determined during the procedure as a known function and cannot be dependent on the spatial location of each part/component of a multibody mechanism. The constraint forces in the CFE method are computed based on the external forces, so the

Materials and Methods
A new form of governing equations is required to apply the algorithm of the problemsolving process. To obtain this new form, three stages need to be done, including the derivation of the governing equations based on the Newton-Euler dynamics, the timediscretization of velocities and accelerations, and rewriting the governing equations based on the generalized coordinates at the new moment of (t n+1 ). The rewritten equations are the function of generalized coordinates and external aerodynamic forces at the present moment of (t n ). By using UDF, the problem is formulated such that the values of aerodynamic forces computed by ANSYS Fluent ® at t n are called in the UDF and replaced as the inputs of the rewritten equations. Therefore, the spatial position of the bodies at the new time-step of (t n+1 ), which is introduced by the generalized coordinates, is calculated. Then, the bodies in ANSYS Fluent ® environment are transmitted to their new positions by UDF coding, and this process continues until the last time-step.

General Problem-Solving Algorithm
The dynamic equations of motion of the bodies were obtained based on the inertial, centrifugal/Coriolis, gravitational, and aerodynamic forces acting on the center of gravity (COG) of the bodies. Then, the equations of motion are discretized, as displayed in Figure 1. The method of multibody system dynamics was used to analyze the dynamics of constrained bodies, and ANSYS Fluent ® software employed for the aerodynamic analysis. aerodynamic equations were coupled in UDF, as the proposed method simulates the whole components of a multibody system straightforward and does not need to apart the components for simulation.

Materials and Methods
A new form of governing equations is required to apply the algorithm of the problem-solving process. To obtain this new form, three stages need to be done, including the derivation of the governing equations based on the Newton-Euler dynamics, the timediscretization of velocities and accelerations, and rewriting the governing equations based on the generalized coordinates at the new moment of ( 1 n t  ). The rewritten equations are the function of generalized coordinates and external aerodynamic forces at the present moment of ( n t ). By using UDF, the problem is formulated such that the values of aerodynamic forces computed by ANSYS Fluent ® at n t are called in the UDF and replaced as the inputs of the rewritten equations. Therefore, the spatial position of the bodies at the new time-step of ( 1 n t  ), which is introduced by the generalized coordinates, is calculated. Then, the bodies in ANSYS Fluent ® environment are transmitted to their new positions by UDF coding, and this process continues until the last time-step.

General Problem-Solving Algorithm
The dynamic equations of motion of the bodies were obtained based on the inertial, centrifugal/Coriolis, gravitational, and aerodynamic forces acting on the center of gravity (COG) of the bodies. Then, the equations of motion are discretized, as displayed in Figure  1. The method of multibody system dynamics was used to analyze the dynamics of constrained bodies, and ANSYS Fluent ® software employed for the aerodynamic analysis. Once the dynamic equations are derived, they can be introduced in the general form presented by Equation (1) as follows:

M(y)y + C(y, y) + G(y) = R
(1) Once the dynamic equations are derived, they can be introduced in the general form presented by Equation (1) as follows:
where y is the vector of generalized coordinates corresponding to the number of system degrees of freedom; M(y) is the symmetric, positive, and definite mass matrix; C(y, . y) is the vector of generalized centrifugal and Coriolis forces; G(y) is the vector of generalized gravitational forces; and R is the vector of generalized aerodynamic forces. Therefore, the system dynamics can be represented by the generalized coordinates (y), their first and second derivatives ( forces. The equations were discretized by the Euler method and represented by the generalized coordinates at times of t n and t n+1 , and the aerodynamic forces at time of t n . Then, the new system of equations was solved for the generalized coordinates at the new time step (y n+1 ) and represented as a function of generalized coordinates and forces at the previous time step (y n , R n ). At the initial moment (n = 0), the spatial position, velocity of the system, and aerodynamic forces are known, i.e., the initial condition is determined. Therefore, by replacing the initial condition in the new system of equations, it is possible to compute the generalized coordinates at the new time-step incrementally. The process continued until the desired time level was achieved. The coding of this algorithm was written in the UDF of ANSYS Fluent ® .
To represent Equation (1) at the time-step of n + 1 as a function of previous time-step (n), the forward finite-difference has been used, as follows: ..
where ∆t is the time interval. A more precise numerical solution requires smaller time intervals. Since the initial condition includes the initial position and velocities, at the time of t = 0, i.e., n = 0 the values of y 0 and . y 0 are known. Therefore, the value of y 1 can be computed by Equation (2). Finally, the values of generalized coordinates in the new time step (y 2 ) can be computed by replacing Equations (2) and (3) into Equation (1) and solving Equation (1) for y 2 . Similarly, the entire trajectory of the motion can be calculated.
Since the discretization was implemented for the coding process in ANSYS Fluent ® software, and the outputs of this software in the UDF environment are velocities and not positions, the discretization was carried out for the velocities of generalized coordinates, as follows: ..
Therefore, Equation (1) should be solved for . y n+1 as a function of y n , . y n , and R n ; and the same incremental procedure was applied.

Algorithm and Problem-Solving for a Four-Degree of Freedom System
A practical example of a planner system with four-degree of freedom based on the algorithm is shown in Figure 1. All the steps of dynamic and aerodynamic analyses and coding in the UDF environment were investigated. This example approves the credibility of the algorithm.
To emphasize the proposed algorithm, a simple model was chosen. This model includes two bodies with the masses of m 1 = 2 kg and m 2 = 2 kg, the lengths of l 1 = 0.5 m and l 2 = 0.5 m, the center of gravity (COG) lengths of l c 1 = 0.25 m and l c 2 = 0.25 m, and the principal inertias of I 1 = 1.4 kg·m 2 and I 2 = 1.4 kg·m 2 . For the CFD solution, the bodies were imposed to the air flow with velocity of v = −3 m/s and v = −5 m/s for constrained and unconstrained systems, respectively. These bodies were connected to each other by a revolute joint. Due to the flow of fluid, aerodynamic forces and torques were applied to the COGs as presented in Figure 2a. In this figure, the vectors of [X 1o/2o , Y 1o/2o ], [X 1I/2I , Y 1I/2I ], and [X 1/2 , Y 1/2 ] show the origin, insertion, and COG of the links 1 and 2, respectively. The vectors of R 1 and R 2 demonstrate the external aerodynamic forces acting on the links of 1 and 2, respectively. The values of d 1 and d 2 are the distances between the applied external forces and the COGs of links 1 and 2. It is desired to compute the planer configurations of the two bodies, in which the following equality is satisfied during the simulation: the planer configurations of the two bodies, in which the following equality is satisfied during the simulation: The satisfaction of Equation (5) would guarantee the connection of two bodies as a revolute joint. The dynamic equations of the motion can be derived using the Newton-Euler method for the free body diagram of the system as indicated in Figure 2b. In this figure, the external forces of where: The satisfaction of Equation (5) would guarantee the connection of two bodies as a revolute joint.
The dynamic equations of the motion can be derived using the Newton-Euler method for the free body diagram of the system as indicated in Figure 2b. In this figure, the external forces of R 1 and R 2 are transmitted to the COGs for bodies 1 and 2, respectively. This creates the torques of τ 1 and τ 2 . The external forces were rewritten based on their components (R 1x , R 1y , R 2x and R 2y ). The values of F rx and F ry are the internal forces in the directions of x and y, and the variables of θ 1 and θ 2 are the angles of bodies 1 and 2 with respect to the horizon, respectively. The system of equations of motion for Figure 2 is presented as follows: (q 1 sin(θ 1 ) − q 2 cos(θ 1 ))l 1c + τ 1 = I 1 . . where: . .
If the external forces and torques were determined, there are four unknown variables (generalized coordinates), including X 1 , Y 1 , θ 1 and θ 2 , which can be computed by solving the system of motion Equations (6)- (9).
The dynamics governing Equations (6)-(9) can be generally solved by two different numerical methods. One is to use Runge-Kutta fourth order to solve the system of ordinary differential equations (ODEs). To use this method, the external forces as known functions of time and generalized coordinates are required. Since in this study, an incremental solution has been adopted due to the dependency of spatial configurations on the aerodynamic forces and the nature of the numerical solution of CFD, another method is proposed. The alternative method is based on the discretization of the system dynamics and reformulation of the equations of motion by the algorithm presented in Figure 1. The first derivative of generalized coordinates at t n+1 can be calculated by using Equation (4) for the discretization of the second derivatives of generalized coordinates and solving Equations (6)-(9) for the time-step of n + 1. The system of equations was solved by Wolfram Mathematica ® software (Wolfram Research, Champaign, IL, USA), and the code in detail are described in Supplementary file A (Solving of the equations of the motions). The code itself is also attached to the manuscript for the readers' information.
To compare the results of the proposed method, the problem was solved by both methods for an arbitrary external force field, and the results were compared. The Matlab ® ODE was applied for the first method and the ANSYS Fluent ® environment using UDF for the second one. The external force field in detail is presented in Supplementary file B (Arbitrary external force field).

Results
The streamlines of the multibody mechanisms in the constrained and unconstrained conditions at different simulation times are illustrated in Figure 3. In the unconstrained condition, there were no constrained forces acting on the mechanisms' joint, and the system was moved only under the aerodynamic forces from the fluid, which applies to the mechanism from right to left. Both bodies 1 and 2 at the beginning of the simulation were placed at the same displacements and angles in respect to each other. The results revealed that the unconstrained mechanisms due to their body weight distributions, which act at their centers, tend to change their alignment into the vertical position. The spatial position of body 2 as the aerodynamic forces act from right to left also influences the position of body 1 in a way that at the simulation time of 0.48 s, body 1 rotates, and consequently brings about a considerable distance between two bodies of 1 and 2. However, when it comes to the constrained mechanisms, the most important point is that two bodies kept their inner distance until the end of the simulation time. Another point is that the role of the spatial location of body 2 here does not become dominant in changing the location of body 1. It implies that in a constrained system, the bodies show more realistic moves under the aerodynamic forces, which is preferable in CFD simulations.
The contours of the velocity of the multibody mechanisms in the constrained and unconstrained conditions at different simulation times are illustrated in Figure 4. In the unconstrained condition, since the static pressure is placed in front of body 2 due to its vertical location, its angle is increased by the passage of simulation time, and it led to a conversion of pressure in the front and back sides of the body. Regarding body 1, after changing its spatial condition to horizontal, the same pressure conversions in both the front and back sides of this body occurred. This issue is observable at the simulation times of 0.48 s and 1.14 s between the unconstrained and constrained bodies of 1 and 2. At the simulation time of 0.48 s, in the unconstrained condition, the highest velocities were at the end of bodies 1 and 2, while in the constrained condition, the highest velocities were at the joint of bodies 1 and 2. This implies that in unconstrained mechanisms, the distribution of velocities were dramatically affected by the aerodynamic forces acting on the system, whereas in a constrained system, the joint can well manage the aerodynamic forces following by a realistic drop mode analyses.
The aerodynamics forces and torques acted on the center of gravities of bodies 1 and 2 in the X and Y directions for both the unconstrained and constrained conditions are presented in Figure 5. Under the aerodynamic forces in the unconstrained condition, at the simulation times of 0.6 s for body 1 and 0.25 s for body 2, the motion of the mechanism was changed. In these simulation times, the bodies were at the minimum velocity values. By changing the directions of the force vectors, the bodies of 1 and 2 relatively deviate to left, which is related to the forces in Y direction. These forces are always positive with the maximum value of 20 N. Meanwhile, as the gravity force of the bodies is −20 N, the summation of these two forces led to a negative acceleration in Y direction for the unconstrained mechanism. In the constrained mechanism, in most of the simulation time the forces in X direction were negative, which is related to the direction of the aerodynamic forces acting on the mechanism. When these forces are being compared to those of unconstrained mechanism, their values are two times bigger, which related to the velocity of the flow here (−5 m/s compared to −3 m/s for the unconstraint mechanism). However, the forces in the Y direction were not altered (~20 N). The momentums here with the passage of simulation time were altered from positive to negative values, which address the rotation of the bodies. The contours of the velocity of the multibody mechanisms in the constrained and unconstrained conditions at different simulation times are illustrated in Figure 4. In the unconstrained condition, since the static pressure is placed in front of body 2 due to its vertical location, its angle is increased by the passage of simulation time, and it led to a conversion of pressure in the front and back sides of the body. Regarding body 1, after changing its spatial condition to horizontal, the same pressure conversions in both the front and back sides of this body occurred. This issue is observable at the simulation times of 0.48 s and 1.14 s between the unconstrained and constrained bodies of 1 and 2. At the simulation time of 0.48 s, in the unconstrained condition, the highest velocities were at the end of bodies 1 and 2, while in the constrained condition, the highest velocities were at the joint of bodies 1 and 2. This implies that in unconstrained mechanisms, the distribution of velocities were dramatically affected by the aerodynamic forces acting on the system, whereas in a constrained system, the joint can well manage the aerodynamic forces following by a realistic drop mode analyses. The X and Y forces acting on the joint of bodies 1 and 2 of the constrained mechanism are displayed in Figure 6. These forces were calculated by derivation of the V x and V y (velocities) at the center of gravities of bodies 1 and 2 and eventually the A x and A y (accelerations). The positive value of F x implies that the force acts from left to right, while the negative value shows that the force acts from right to left at the joint. The results revealed that the motion of body 2 is dominant and body 1 just is following that. The positive value of F y at the center of gravity of body 1 implies that the force acts from the down to up. However, since F y is almost negative during the simulation time it means that body 1 experiences a higher force from the upper side, which leads to a rotation in that and take it to an upper right position compared to body 2. The aerodynamics forces and torques acted on the center of gravities of bodies 1 and 2 in the X and Y directions for both the unconstrained and constrained conditions are presented in Figure 5. Under the aerodynamic forces in the unconstrained condition, at the simulation times of 0.6 s for body 1 and 0.25 s for body 2, the motion of the mechanism was changed. In these simulation times, the bodies were at the minimum velocity values. By changing the directions of the force vectors, the bodies of 1 and 2 relatively deviate to left, which is related to the forces in Y direction. These forces are always positive with the maximum value of 20 N. Meanwhile, as the gravity force of the bodies is −20 N, the summation of these two forces led to a negative acceleration in Y direction for the unconstrained mechanism. In the constrained mechanism, in most of the simulation time the forces in X direction were negative, which is related to the direction of the aerodynamic forces acting on the mechanism. When these forces are being compared to those of unconstrained mechanism, their values are two times bigger, which related to the velocity of the flow here ( 5 m/s  compared to 3 m/s  for the unconstraint mechanism). However, the To verify the developed UDF code, an arbitrary flow moving from right to left (Supplementary file B: Arbitrary external force field) was simulated using Matlab ® ODE, and the results were compared to that of the UDF of ANSYS Fluent ® . Since in Matlab ® ODE the problem was simulated on the basis of the fourth order Runge-Kutta, the resulted error was three orders smaller than that of the UDF. The axial and angular velocities for both bodies 1 and 2 were computed using UDF and Matlab ® ODE, and the results and their deviations are plotted in Figure 7. The mean percentage errors between the UDF and Matlab ® ODE for the V x , V y , and ω of body 1 were 1.46%, 0.84%, and 1.87%, respectively and for body 2 were 1.15%, 0.82%, and 0.43%, respectively. These trivial errors imply the appropriateness of UDF code in solving a complex multibody mechanism.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 9 of 13 forces in the Y direction were not altered (~20 N). The momentums here with the passage of simulation time were altered from positive to negative values, which address the rotation of the bodies. The X and Y forces acting on the joint of bodies 1 and 2 of the constrained mechanism are displayed in Figure 6. These forces were calculated by derivation of the F implies that the force acts from left to right, while the negative value shows that the force acts from right to left at the joint. The results revealed that the motion of body 2 is dominant and body 1 just is following that. The positive value of y F at the center of gravity of body 1 implies that the force acts from the down to up. However, since y F is almost negative during the simulation time it means that body 1 experiences a higher force from the upper side, which leads to a rotation in that and take it to an upper right position compared to body 2. To verify the developed UDF code, an arbitrary flow moving from right to left (Supplementary file B: Arbitrary external force field) was simulated using Matlab ® ODE, and the results were compared to that of the UDF of ANSYS Fluent ® . Since in Matlab ® ODE the problem was simulated on the basis of the fourth order Runge-Kutta, the resulted error was three orders smaller than that of the UDF. The axial and angular velocities for  The X and Y forces acting on the joint of bodies 1 and 2 of the constrained mechanism are displayed in Figure 6. These forces were calculated by derivation of the F implies that the force acts from left to right, while the negative value shows that the force acts from right to left at the joint. The results revealed that the motion of body 2 is dominant and body 1 just is following that. The positive value of y F at the center of gravity of body 1 implies that the force acts from the down to up. However, since y F is almost negative during the simulation time it means that body 1 experiences a higher force from the upper side, which leads to a rotation in that and take it to an upper right position compared to body 2. To verify the developed UDF code, an arbitrary flow moving from right to left (Supplementary file B: Arbitrary external force field) was simulated using Matlab ® ODE, and the results were compared to that of the UDF of ANSYS Fluent ® . Since in Matlab ® ODE the problem was simulated on the basis of the fourth order Runge-Kutta, the resulted error was three orders smaller than that of the UDF. The axial and angular velocities for Selection of the time step size is one of the most important factors in the numerical integration of differential equation systems. To determine the appropriate size of the time step, first, the problem was solved by Matlab ® ODE for an arbitrary force field with high frequencies (Supplementary file B: Arbitrary external force field) compared to a real force field. The method of ODE45 was selected. All ODE solvers provided by Matlab ® employ a variety of variable-step methods. The Matlab ® ODE solvers apply these methods by taking a step, computing the error at this step, checking if the value is greater than or less than the tolerance, and changing the step size accordingly. Utilizing a variable step guarantees that a large step size is utilized for low frequencies and a small step size is utilized for high frequencies. The ODE solvers within MATLAB are improved for a variable step, run faster with a variable step size, and clearly, the outcomes are more accurate. After solving the problem using ODE45 the smallest time step size was selected, and the differential equations of motion Equations (6)-(9) were solved by forward finite-difference Equations (2) and (3) in UDF environment to be sure that the error is less than the tolerance. The problem was solved for both unrealistic and realistic force fields. The reason a variable time step size was not selected for UDF, was that the ANSYS Fluent ® would not accept variable values for time intervals, but only fixed step size. The proposed methodology can be described in short as computing equations of motion in the presence of external and internal forces applying on a mechanism as aerodynamic and constraint forces, respectively. Then, discretizing and reformulating them based on the next time step. Therefore, the next positions can be estimated by the previous values of positions and forces. The difficult part is to compute the equations of motion. Discretizing and reformulating the equations can be simply done by the Matlab ® and Mathematica ® . Applying the new discretized and reformulated equation in UDF can be a little difficult [25]. Even though UDF programming is difficult, the physical and mathematical equations used are still based on Newton's law and can be discretized.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 10 of 13 both bodies 1 and 2 were computed using UDF and Matlab ® ODE, and the results and their deviations are plotted in Figure 7. The mean percentage errors between the UDF and Matlab ® ODE for the x V , y V , and  of body 1 were 1.46%, 0.84%, and 1.87%, respectively and for body 2 were 1.15%, 0.82%, and 0.43%, respectively. These trivial errors imply the appropriateness of UDF code in solving a complex multibody mechanism. Selection of the time step size is one of the most important factors in the numerical integration of differential equation systems. To determine the appropriate size of the time step, first, the problem was solved by Matlab ® ODE for an arbitrary force field with high frequencies (Supplementary file B: Arbitrary external force field) compared to a real force field. The method of ODE45 was selected. All ODE solvers provided by Matlab ® employ a variety of variable-step methods. The Matlab ® ODE solvers apply these methods by taking a step, computing the error at this step, checking if the value is greater than or less than the tolerance, and changing the step size accordingly. Utilizing a variable step guarantees that a large step size is utilized for low frequencies and a small step size is utilized for high frequencies. The ODE solvers within MATLAB are improved for a variable step, run faster with a variable step size, and clearly, the outcomes are more accurate. After solving the problem using ODE45 the smallest time step size was selected, and the differential equations of motion Equations (6)-(9) were solved by forward finite-difference Equations (2) and (3) in UDF environment to be sure that the error is less than the tolerance. The problem was solved for both unrealistic and realistic force fields. The reason a Figure 7. The comparison of UDF and Matlab ® ODE, including (a) velocity in X direction of body 1, (b) velocity in Y direction of body 1, (c) angular velocity of body 1, (d) velocity deviation between UDF and Matlab ® ODE in X direction of body 1, (e) velocity deviation between UDF and Matlab ® ODE in Y direction of body 1, (f) angular velocity deviation between UDF and Matlab ® ODE of body 2, (g) velocity in X direction of body 2, (h) velocity in Y direction of body 2, (i) angular velocity of body 2, (j) velocity deviation between UDF and Matlab ® ODE in X direction of body 2, (k) velocity deviation between UDF and Matlab ® ODE in Y direction of body 2, and (l) angular velocity deviation between UDF and Matlab ® ODE of body 2.

Limitations
The study is limited by the following considerations. First, a 2D model was used to validate the accuracy of the proposed algorithm, which might not be a real case of a flying object in the air. However, the proposed approach can simply be extended to a 3D model, which will be the subject of a future report. Second, the solution was fulfilled through a stepwise approach. In this approach, the forces that are being calculated at each step are being used to calculate the location of the links at one-step ahead of the time. However, this could impose a step error into the simulation results, having a very small ∆t would minimize the error in the simulation results. The experimental validation of this work was also a challenge for us due to the lack of enough resources and instruments. This is the authors' objective to extend the code for a 3D model and experimentally verify which all will be the subject of a future report. The UDF code of the proposed algorithm is freely available to use as supplementary material (Supplementary file C: UDF code) and can extend our understanding of the mechanisms of the flying objects in the air.

Conclusions
This study proposed a new UDF coding to analyze the dynamics of a multibody mechanism with different degree-of-freedom in ANSYS Fluent ® software. We could address the shortcomings of the proposed models so far, especially CFE, as it unable to model a multibody mechanism in the presence of various aerodynamics loadings related to the spatial configurations of the bodies as fulfilled in ANSYS Fluent ® software. The mean percentage errors between the UDF and Matlab ® ODE for the V x , V y , and ω of body 1 were found to be 1.4%, 0.8%, and 1.8%, respectively, and for body 2 were 1.1%, 0.8%, and 0.4%, respectively. The negligible errors found here represent the powerfulness of the proposed UDF code in addressing the complex mechanical response of a multibody mechanism flying in the air. In addition, the proposed models so far have used coupling of ANSYS Fluent ® and Matlab ® to deal with the equations; however, the simulation time in this case dramatically increases. The proposed method can directly address this issue using a single programming code, namely, UDF of ANSYS Fluent ® . In this method, the equations of motion were discretized on the basis of the velocities and accelerations of a multibody mechanism under the aerodynamic forces. The equations were solved incrementally according to the location of the rigid body and its aerodynamic forces acting on it. The equations were first solved in ANSYS Fluent ® and then updated at each increment in the UDF code. The proposed method here has implications in analyses of flying robots, airplanes, birds, flies, etc., to provide a more accurate dynamic concept of their movements and induced forces during their flies.