A method for modeling electromechanical coupling dynamics of an electric vertical take-off and landing aircraft based on differential-algebraic equations

Abstract

A kind of electromechanical coupling dynamics modeling method of electric vertical take-off and landing aircraft based on differential algebraic equation belongs to electric vertical take-off and landing aircraft technical field.First, the motor dynamics equation in eVTOL electric propulsion system is established according to motor parameters;Second, according to the geometric configuration and inertia parameters of eVTOL, the dynamics equation of fuselage multibody system is established, and the aerodynamic model is applied;Third, the motion driving equation is used to couple based on the general mechanical load calculation scheme of Lagrange multiplier;Finally, the electromechanical coupling dynamics model of eVTOL is established, which is discretized by trapezoidal-Newmark hybrid discrete format, and efficient solution is realized based on Newton iteration method, and it is verified by eVTOL vertical take-off.The invention avoids the defects of traditional simplified model ignoring motor transient characteristics and non-coupling modeling without real-time load feedback, improves the motion response and energy consumption calculation accuracy, and provides reliable dynamics modeling basis for eVTOL energy consumption optimization and commercial design.

Description

Technical Field

[0001] This invention belongs to the field of electric vertical takeoff and landing (EVTOL) aircraft technology, and relates to an electromechanical coupling dynamics modeling method for EVTOL aircraft based on differential algebraic equations. Background Technology

[0002] eVTOL (Electric Vertical Take-Off and Landing) aircraft, as one of the most innovative technologies in urban air mobility, has become a research hotspot in industry and academia in recent years. Chinese invention patent CN120068335A points out that against the backdrop of a continuously increasing global car ownership and increasing pressure on ground transportation systems, eVTOL offers a transformative path to solving urban traffic congestion. With its outstanding advantages such as no runway required and environmental friendliness, eVTOL is not only a cutting-edge mode of air travel but is also considered an important carrier of the low-altitude economy, currently in a crucial stage of moving from proof-of-concept to initial commercialization.

[0003] The core research directions of eVTOL cover overall size design and selection, battery system design, flight path planning, and control system development, as exemplified by Chinese invention patents CN120595615A and CN120595616A. In eVTOL research, flight energy consumption is a key factor to consider. This is because urban air traffic has stringent requirements for operating costs and sustainability, making high energy efficiency a core prerequisite for technology implementation. Energy consumption optimization allows eVTOLs to extend their cruising range and hovering time with limited battery capacity, meeting the basic needs of commercial operation; conversely, excessive energy consumption directly increases battery weight and cost, limits payload, and even threatens the safety and lifespan of the eVTOL. Therefore, energy consumption level is both a core indicator for measuring the maturity of eVTOL technology and an important benchmark for its commercial feasibility. Currently, rich research results have been accumulated in eVTOL flight energy consumption calculation, but existing models still have significant limitations in characterizing the dynamic characteristics of the system, and these limitations exhibit a hierarchical nature, progressing from superficial to profound.

[0004] First, many studies employ simplified power models that fail to incorporate the motor's inherent dynamic characteristics into their calculations. These models typically treat the motor as an ideal power source, assuming a transient linear relationship between voltage, current, and output power, while neglecting the inherent dynamic response process of the motor's physical system. In actual operating conditions, when voltage commands change, there is a significant lag in the motor's current build-up and speed increase. Even with a sudden increase in voltage, the motor's output power cannot immediately increase in sync; instead, it undergoes a dynamic phase of "gradual current increase - slow electromagnetic torque build-up - rotor speed lagging behind." This lag effect is further amplified during transient conditions such as takeoff, landing, hovering, and steering: due to drastic fluctuations in thrust demand, the frequent adjustments to voltage commands and the accumulated deviation between the actual voltage output and the actual power output lead to a significant discrepancy between the energy consumption calculations based on ideal power models and the actual energy consumption.

[0005] Furthermore, while existing research has overcome the limitations of simplified power models by introducing motor dynamics and realizing motor-rotor coupling calculations, key research gaps remain. eVTOL is essentially a strongly coupled "electromechanical-aerodynamic" system. When the motor drives the rotor to rotate, the rotor feeds back two types of loads to the motor through the transmission system: one is the rotor inertial force during the speed change phase, and the other is the aerodynamic counter-torque under the action of airflow. These two types of loads directly change the motor's electrical response, thereby affecting energy consumption. If electromechanical separation modeling is used, this bidirectional coupling effect cannot be transmitted in real time: misestimating or even ignoring the additional effects of rotor inertial force and aerodynamic counter-torque on motor current will not only cause response prediction deviations of key parameters such as current and speed, but will also further cause calculation errors in flight trajectory, forming a chain problem of "inaccurate load calculation - response deviation - misjudgment of energy consumption". However, the core objective of existing coupling modeling research is mostly focused on the impact of coupling on response characteristics, handling quality, or control logic, and has not yet systematically clarified at the positive dynamics level: the motion response prediction deviation caused by ignoring the motor-rotor dynamic coupling effect, and how this deviation will affect the accuracy of energy consumption calculation. In other words, existing work has not quantitatively demonstrated the deviation of the two modeling methods of "considering coupling" and "ignoring coupling" in motion response prediction, nor has it further linked the quantitative relationship between the deviation and the energy consumption calculation results, and has failed to reveal the specific contribution mechanism of the coupling effect to the energy loss components such as motor copper loss, iron loss, and friction loss - even if a coupling model is established, its core application scenario is not focused on energy consumption analysis. Summary of the Invention

[0006] To address the aforementioned technical problems, this invention proposes a method for electromechanical coupling dynamics modeling of electric vertical takeoff and landing (eVTOL) aircraft based on differential algebraic equations (DAE). First, this invention constructs a dynamic model of a permanent magnet synchronous motor and a fuselage multibody system dynamic model as precise sub-models. Second, relying on the bidirectional action mechanism of "motor drive-load feedback," a DAE system is integrated; and a general mechanical load calculation paradigm is constructed based on the Lagrange multiplier method, enabling load solution without decomposing the physical causes. Finally, a trapezoidal-Newmark hybrid discretization scheme is proposed to match first-order electrical variables and second-order mechanical variables, ensuring consistency in electromechanical motion coupling while achieving efficient solution of the DAE system. This method constructs a unified eVTOL electromechanical coupling framework, naturally integrating multi-physics strong coupling effects and motor energy losses into the model, improving adaptability to diverse configurations, and simultaneously solving the deficiency of uncoupled models in transmitting real-time load feedback.

[0007] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0008] A method for electromechanical coupling dynamics modeling of electric vertical takeoff and landing aircraft based on differential algebraic equations includes the following steps:

[0009] Step 1: Based on the electrical and mechanical motion parameters of the motor, establish the motor dynamics equations in the eVTOL electric propulsion system; specifically:

[0010] The surface-mounted permanent magnet synchronous motor was selected as the research object. Formula (1) is the dynamic equation of the surface-mounted permanent magnet synchronous motor considering copper loss and iron loss, which lays the foundation for subsequent electromechanical coupling dynamic analysis.

[0011] (1)

[0012] In the formula, The electromagnetic driving torque of a permanent magnet synchronous motor. This represents the number of pole pairs in the stator winding of the motor. This is the excitation flux of the rotor permanent magnet; For leakage sensation, For magnetizing inductance, and They are permanent magnet synchronous motors Stator inductance along the d-axis and q-axis in the coordinate system; and These are the voltages along the d-axis and q-axis, respectively. and These represent the d-axis and q-axis currents, respectively; t represents time. and These are the d-axis and q-axis iron loss current components, respectively; and These are the stator resistance and the core loss resistance, respectively. This refers to the rotor's mechanical angular velocity; This refers to the rotor's rotational angle position; It is the moment of inertia of the motor rotor; The viscosity coefficient of the motor bearing; This refers to the mechanical load fed back to the motor by the driven mechanism.

[0013] Further constructing a general dynamic model of the motor, the motor dynamic equations of the general dynamic model essentially include a circuit module and a mechanical motion module, which achieve energy and signal coupling and transfer through electromagnetic torque. Its general expression is written as:

[0014] (2)

[0015] The first equation describes the electrodynamic behavior of the motor, and the second equation describes the mechanical dynamic behavior of the motor. Electrical variables representing the motor Represents the mechanical motion variables of the motor. and These represent the first and second derivatives of the variable with respect to time, respectively; This represents the input voltage of the motor. This represents the mechanical load transmitted to the motor by the multibody system of the fuselage; and These are the nonlinear mapping functions for the electrodynamics and mechanical dynamics of the motor, respectively, and their specific forms are determined by the motor type and structural parameters. Step 2: Based on the geometric configuration and inertia parameters of the eVTOL, establish the dynamic equations of the eVTOL's fuselage multibody system; specifically:

[0016] This invention uses a variable-speed direct-drive quadrotor eVTOL for illustration. However, it should be noted that the modeling method employed is also applicable to other configurations. The eVTOL's multi-body system is decoupled into seven independent components, specifically including the main fuselage, left nacelle (containing the left motor), left rotor, right nacelle (containing the right motor), right rotor, upper tail rotor, and lower tail rotor (each corresponding to a built-in motor). To accurately describe the spatial motion of each rigid body, a right-handed coordinate system is established with the center of mass of each rigid body as the origin. Simultaneously, to simplify the modeling process while ensuring analytical accuracy, the following assumptions are made regarding flight characteristics:

[0017] (1) The Earth is considered an inertial frame of reference;

[0018] (2) All components of the fuselage multibody system are considered as rigid bodies, and their elastic deformation is not considered;

[0019] (3) The positions of the centers of mass of each rigid body do not change during flight;

[0020] (4) The fuselage multibody system is symmetrically distributed along the longitudinal plane of symmetry;

[0021] (5) This invention only considers the vertical takeoff phase, and therefore does not consider the aerodynamic characteristics of the wing or the aerodynamic coupling between the rotor and the wing.

[0022] (6) This invention does not consider the movement of the air itself, and the speed of the fuselage in the global coordinate system is equal to the speed of the incoming airflow.

[0023] The first in the fuselage multibody system The position of each component in space is determined by the position vector of its centroid. This indicates that the position vector represents the coordinates of the component's center of mass in the global inertial coordinate system OXYZ; the Cardan angle is chosen here. The position vector plus the rotation coordinates are used to describe the rotation of the rigid body's body-following coordinate system relative to the global inertial coordinate system. Generalized coordinates of the position and orientation of each component in space ,in Indicates the first The coordinates of the centroid position of each component. Indicates the first The rotational coordinates of each component. By combining the generalized coordinates of all components of the fuselage multibody system, a generalized coordinate vector describing the entire fuselage is obtained:

[0024] (3)

[0025] in, This represents the generalized coordinate vector of the fuselage multibody system. This represents the generalized coordinate vector of component number 1. This represents the generalized coordinate vector of component number 2. This represents the generalized coordinate vector of component number 7;

[0026] After obtaining the generalized coordinate vector of the entire fuselage multibody system, the connection relationship between the various rigid body components is explained next. The four rotors at the front and rear are all connected to their supporting structure by rotating joints, and the left and right nacelles are connected to the fuselage body by fixed support. After defining the generalized coordinates and connection methods, the dynamic equations of the fuselage multibody system are obtained by using the Cartesian modeling method, as shown in formula (4);

[0027] (4)

[0028] in, , , These represent the generalized coordinates, generalized velocity, and generalized acceleration of the fuselage multibody system, respectively. It is the total mass matrix of the fuselage multibody system, which is obtained by combining the mass matrices of each component. It is a Lagrange multiplier. It's time. These are algebraic constraint equations, representing the connection relationships between various components. It is the Jacobian matrix of the constraint equations with respect to the generalized coordinates. It refers to the generalized external forces acting on the fuselage multibody system.

[0029] Step 3: Apply an aerodynamic model that is both easy to calculate and possesses considerable accuracy to the fuselage multibody system dynamic equations established in Step 2; specifically:

[0030] Because eVTOL differs from traditional helicopters in blade aspect ratio, aerodynamic modeling can be appropriately simplified. The formulas used to calculate blade thrust and counter-torque are as follows:

[0031] (5)

[0032] (6)

[0033] in, This refers to the thrust of the propeller blades. It is the reverse torque. This represents the blade rotation speed. The diameter is the blade diameter. and These are the dimensionless tensile force coefficient and torque coefficient, respectively. It refers to the air density in the flight environment, which is related to altitude. and temperature The function;

[0034] (7)

[0035] Among them, standard atmospheric density Atmospheric pressure This can be further expressed as:

[0036] (8)

[0037] The calculation formulas for blade thrust and counter-torque of the eVTOL under vertical takeoff and landing conditions are as follows:

[0038] (9)

[0039] (10)

[0040] in:

[0041] (11)

[0042] The parameters in equations (9)-(11) are defined as follows: This is the correction factor for the effective area of ​​the blade. It refers to the number of blades. It is the lift line slope coefficient of the airfoil. It is a correction factor resulting from the downwashing effect. It is the blade pitch. It is the diameter of the blade. It is zero-lift angle of attack. It is the aspect ratio. It is the average radius proportionality coefficient, used to estimate the average incoming flow velocity on the blade. It is the airfoil drag coefficient of the blade section. It is the zero-lift drag coefficient. It is the Oswald factor. This is the aircraft's vertical velocity in the eVTOL direction, without considering the motion of the air itself. It also represents the relative speed with the airflow during eVTOL vertical takeoff.

[0043] Step 4: Couple the dynamic models from Steps 1 to 3 using the motion driving equations and a general mechanical load calculation scheme based on Lagrange multipliers; specifically:

[0044] Step 4-1: Use motion drive equations to describe the control commands of the motor established in Step 1 to the fuselage multibody system established in Step 2.

[0045] Step 4.2: Using a general mechanical load calculation scheme based on Lagrange multipliers, calculate the mechanical load on the motor established in Step 1 during the operation of the fuselage multibody system established in Step 2. Specifically:

[0046] During the process of the motor providing motion control commands to the eVTOL rotor, it is subjected to a reaction mechanical load from the drive mechanism, and the two together constitute a dynamically coupled mechanical system. The dynamic behavior of the dynamically coupled mechanical system is described by the fifth equation in formula (1). From the perspective of the load composition mechanism, the physical source of the mechanical load on the motor in this scenario is clear and can be directly disassembled, mainly consisting of the rotor's inertial torque and counter-torque. That is:

[0047] (12)

[0048] in, It is the moment of inertia of the rotor about its axis of rotation. It is the reverse torque.

[0049] The rotor inertial torque is directly obtained from the rotor's moment of inertia and angular velocity around the axis of rotation, while the aerodynamic counter-torque is related to the rotor's aerodynamic characteristics. A clear quantitative model is established based on rotor aerodynamic theory or wind tunnel test data as shown in Equation (10). Therefore, for the specific scenario of eVTOL, the mechanical load on the motor can be solved by calculating the rotor inertial torque and aerodynamic counter-torque separately and then performing vector superposition. The overall calculation process has a clear physical correspondence.

[0050] This invention proposes a general method for calculating mechanical loads based on Lagrange multipliers. This method can obtain mechanical loads under various operating conditions without dissecting the physical nature of the load, avoiding redundant derivations and handling multi-source coupling problems. It improves the versatility and adaptability of electromechanical coupling modeling. The general method for calculating mechanical loads is as follows:

[0051] Assume a rigid body Through constraints Connected to another rigid body, the specific hinge point is... The constraint equations are The Lagrange multipliers corresponding to the constraints are Acting on a rigid body The generalized constraint on is The virtual power represented by the generalized constraint force is:

[0052] (13)

[0053] In the formula, It is a variational symbol; It is a rigid body Translation coordinates; It is a rigid body The rotation coordinates; It is virtual power; It is a rigid body Generalized velocity; It is a constraint The corresponding constraint equations for rigid bodies The Jacobian matrix of the generalized coordinates; It is a rigid body The first derivative of the translational coordinates with respect to time; It is a rigid body The first derivative of the rotating coordinates with respect to time.

[0054] At the hinge point Established at the point of being fixed to a rigid body Hinge point coordinate system In the specific analysis process, the coordinate axes of the hinge point coordinate system are defined on the rotor's rotation axis vector. Above, the constraint reactions or driving forces at the hinge points will be described in the hinge point coordinate system. The hinge point coordinate system and the rigid body... Rigid body body coordinate system Directional cosine matrix Let be a constant matrix. The constraint forces are denoted in the hinge point coordinate system as:

[0055] (14)

[0056] In the formula, and Constraints The array of constraint forces and constraint couples in the hinge point coordinate system; Representing a rigid body The projection of the constraint forces onto the hinge point coordinate system.

[0057] Projecting the constraint forces onto the global coordinate system yields:

[0058] (15)

[0059] in, Representing a rigid body The projection of the constraint forces onto the global coordinate system; Representing a rigid body Rotation matrix of the body coordinate system; Indicates that it is fixed to a rigid body Hinge point coordinate system rotation matrix; Representing a rigid body The projection of the constraint forces onto the hinge point coordinate system;

[0060] The virtual power of the constraint force is expressed as:

[0061] (16)

[0062] In the formula, Representing a rigid body upper hinge point The generalized velocity matrix; definition ,and hinge point The projection of the velocity onto the global coordinate system:

[0063] (17)

[0064] in, It is the origin of the body coordinate system The first derivative of the position vector with respect to time; It is a hinge point Relative to the origin of the body coordinate system The radius vector, Indicates the symbol for an antisymmetric matrix; rigid body The angular velocity vector is written as:

[0065] (18)

[0066] in, It is a rigid body The transformation matrix for converting the derivative of the rotating coordinates to angular velocity depends on the chosen rotating coordinates; It is a rigid body The rotation coordinates.

[0067] From equations (17) and (18), we can obtain:

[0068] (19)

[0069] In the formula, Representing a rigid body Generalized velocity; Represents the generalized velocity transformation matrix; Will Transform into a rigid body The transformation matrix for generalized velocity is as follows:

[0070] (20)

[0071] in, Represents the identity matrix; This indicates a hinge point. Relative to the origin of the body coordinate system The antisymmetric matrix formed by the position vectors; Representing a rigid body The transformation matrix for converting the derivative of rotating coordinates to angular velocity is then... Representing a rigid body The inverse of the transformation matrix for converting the derivative of rotating coordinates to angular velocity;

[0072] By comparing equations (13) and (16), we can obtain the array of constraint forces in the hinge coordinate system, which is the rotor mechanical load on the motor rotor in the motor coordinate system.

[0073] Step 5: Based on the coupling scheme in Step 4, establish the eVTOL electromechanical coupling dynamic model, then discretize it using a trapezoidal-Newmark hybrid discretization scheme, and finally solve it using the Newton-Raphson iteration method; specifically:

[0074] Step 5.1: Establish the eVTOL electromechanical coupling dynamic model;

[0075] Combining formula (2) from step 1 and formula (4) from step 2, we obtain the eVTOL electromechanical coupling dynamic model, as shown below:

[0076] (twenty one)

[0077] in, This represents the first derivative of the electrical variable of the motor with respect to time. This represents the second derivative of the mechanical motion variables of the motor with respect to time. Represents the electrical variables of the motor; The nonlinear mapping function representing the electrodynamics of the motor; This indicates the mechanical load transmitted to the motor by the fuselage multibody system; Indicates time; A nonlinear mapping function representing the mechanical dynamics of an electric motor; This represents the total mass matrix of the fuselage multibody system; This represents the generalized acceleration of the fuselage multibody system; It represents the Jacobian matrix of the constraint equations with respect to the generalized coordinates; This indicates that it is a Lagrange multiplier; This refers to the generalized external forces acting on the fuselage multibody system; Represent the constraint equations; Represents the mechanical motion variables of the motor; Represents the generalized coordinates of the fuselage multibody system;

[0078] In formula (21): the first formula describes the first-order differential dynamics of the motor electrical module, the second formula characterizes the second-order differential behavior of the motor mechanical module, the third formula is the multibody dynamics equation of the fuselage, and the fourth formula is the algebraic constraint equation. The eVTOL electromechanical coupling dynamics equation shown in formula (21) contains the first-order differential of the electrical system, the second-order differential of the mechanical system, and algebraic constraints, belonging to DAE. This set of equations is passed through and To achieve the correlation between the motor and the machine body movement, through The reverse load transfer is completed, which preserves the energy loss of the motor and the multi-body dynamic characteristics of the machine body, while the DAE constraint mechanism ensures the real-time performance and accuracy of the coupling, laying the model foundation for subsequent accurate calculation of motion response and energy consumption.

[0079] Step 5.2: Discretize the eVTOL electromechanical coupling dynamics model obtained in Step 5.1 using a trapezoidal-Newmark hybrid discretization scheme;

[0080] The eVTOL electromechanical coupling dynamics model constructed in this invention is ultimately presented in the form of a DAE, which includes the differential variables of the first-order differential equations of the motor electrical module, the second-order differential variables of the motor mechanical module, the second-order dynamic differential variables of the fuselage multibody system, and Lagrange multipliers related to the algebraic constraint equations. Traditional numerical discretization schemes are difficult to apply directly.

[0081] To effectively simulate such equations numerically, this invention constructs a trapezoidal-Newmark hybrid discretization scheme, which, by specifically adapting to the dynamic characteristics of different modules, achieves efficient and accurate solutions for DAE systems. The core idea of ​​this method is: to discretize the first-order differential variables of the motor electrical module using the implicit trapezoidal method, and to discretize the second-order differential variables of the motor mechanical module, the fuselage multibody system, and the Lagrange multipliers related to the algebraic constraint equations using the Newmark method, as shown in formula (22):

[0082] (twenty two)

[0083] Among them, subscript and Represents the start and end nodes of a discrete time period; Represents the length of a discrete interval; and It is a parameter in the Newmark algorithm that controls the stability and accuracy of the algorithm. Electrical variables representing the motor at the end of the time period; Electrical variables representing the motor at the start point of the time period; This represents the first derivative of the electrical variable of the motor at the end of the time interval with respect to time. The first derivative of the electrical variable of the motor with respect to time at the starting point of the time interval; The mechanical motion variable of the motor at the end of the time period; The mechanical motion variables of the motor at the start point of the time period; The first derivative of the mechanical motion variable of the motor with respect to time at the starting point of the time interval; This represents the first derivative of the mechanical motion variable of the motor at the end of the time interval with respect to time. The second derivative of the mechanical motion variable of the motor with respect to time at the starting point of the time interval; This represents the second derivative of the mechanical motion variable of the motor with respect to time at the end of the time interval. The generalized coordinates of the fuselage multibody system representing the end node of the time period; The generalized coordinates of the fuselage multibody system representing the starting point of a time period; The generalized velocity of the fuselage multibody system at the start point of the time period; The generalized acceleration of the fuselage multibody system represents the starting point of the time period; The generalized acceleration of the fuselage multibody system at the end of the time period; This represents the generalized velocity of the fuselage multibody system at the end of the time period;

[0084] The discretization of the motor mechanical variables and the generalized coordinates of the fuselage must use the same set of parameters.

[0085] By employing the discretization strategy shown in Equation (22), the continuous-time domain eVTOL electromechanical coupling dynamics model is transformed into a system of nonlinear algebraic equations in the discrete-time domain. As the basic unknowns, equation (21) is discretized into the following form:

[0086] (twenty three)

[0087] in, and These correspond to the residuals of the electrical and mechanical modules of the motor, respectively. For the multibody dynamics equations residuals of the fuselage, For the algebraic constraint equation residuals, and for algorithm stability, Multiply by the constraint equation.

[0088] Step 5.3: The Newton-Raphson iteration method is used to solve the discrete form of the eVTOL electromechanical coupling dynamics model obtained in Step 5.2;

[0089] The Newton-Raphson iteration method is used to solve the nonlinear algebraic equations shown in formula (23) at each time step. The linear equations to be solved at each iteration step are given below.

[0090] (twenty four)

[0091] Write it in explicit format:

[0092] (25)

[0093] in, express right The Jacobian matrix; express right The Jacobian matrix; express right The Jacobian matrix; express right The Jacobian matrix; express right The Jacobian matrix; express right The Jacobian matrix; express right The Jacobian matrix; express right The transpose of the Jacobian matrix; express right The Jacobian matrix; express right The Jacobian matrix; The iterative increment of the electrical variable of the motor at the end of the time period; The iterative increment of the mechanical motion variable of the motor at the end of the time period; This represents the iterative increment of the generalized acceleration of the fuselage multibody system at the end of the time period; This represents the Lagrange multiplier iteration increment at the end of the time period;

[0094] so:

[0095] (26)

[0096] in, Indicates the previous iteration step; Indicates the next iteration step;

[0097] In summary, this step proposes a trapezoidal-Newmark hybrid discretization scheme for the constructed DAE containing first-order electrical variables, second-order mechanical variables, and algebraic constraints. The trapezoidal method is used to discretize the first-order electrical variables of the motor, while the Newmark method is used to discretize the second-order variables of the motor-mechanical and fuselage multibody systems. These are then uniformly transformed into a system of nonlinear algebraic equations and solved using Newton iterations. Simultaneously, the consistency of the discretized parameters of the motor and fuselage mechanical variables is ensured to satisfy the constraints. This scheme provides a suitable numerical solution for electromechanical coupling models, achieving coordinated calculation of variables of different orders and effective satisfaction of constraints, laying the foundation for the simulation implementation of the model.

[0098] The beneficial effects of this invention are as follows:

[0099] This invention constructs an eVTOL electromechanical coupling dynamics model based on DAE (Design for Electric Vehicles), organically integrating motor dynamics (considering energy losses such as copper and iron losses) with fuselage multibody dynamics and aerodynamic calculations. Simultaneously, it introduces the Lagrange multiplier method to achieve a unified solution for mechanical loads, ensuring efficient and high-precision dynamic input for eVTOL motion response and energy consumption analysis. By employing a trapezoidal-Newmark hybrid discretization scheme to match first-order electrical variables with second-order mechanical variables, the solution efficiency and stability of the DAE system are significantly improved. Combining the "electromechanical-aerodynamic" bidirectional coupling mechanism with real-time load feedback effectively avoids the accumulation of deviations between traditional simplified power models and uncoupled modeling, ensuring efficient generation of accurate flight trajectories and energy consumption predictions even under complex conditions such as vertical takeoff and hovering maneuvers. Attached Figure Description

[0100] Figure 1 This is a flowchart of the present invention.

[0101] Figure 2 This is a circuit model diagram of the motor of the present invention; Figure 2 (a) is the circuit diagram of the motor d-axis in this invention. Figure 2 (b) is the circuit diagram of the q-axis of the motor in this invention.

[0102] Figure 3 This is a diagram of the eVTOL model of the present invention; Figure 3 (a) is the fuselage configuration diagram of the eVTOL in this invention. Figure 3 (b) is a body coordinate system diagram of each component of the eVTOL fuselage in this invention.

[0103] Figure 4 This is a schematic diagram of motion drive in this invention.

[0104] Figure 5 This is a schematic diagram of the mechanical load in this invention.

[0105] Figure 6 This is a flowchart of the simulation calculation process for this invention.

[0106] Figure 7 This is the voltage input time history curve of the present invention.

[0107] Figure 8 This is a trajectory diagram of the fuselage altitude in this invention.

[0108] Figure 9 This is a diagram showing the energy consumption distribution of each part in this invention. Detailed Implementation

[0109] The present invention will be further described below with reference to specific embodiments.

[0110] A method for electromechanical coupling dynamics modeling of electric vertical takeoff and landing aircraft based on differential algebraic equations, such as... Figure 1 As shown, it includes the following steps:

[0111] Step 1: Based on the electrical and mechanical motion parameters of the motor, establish the motor dynamics equations in the eVTOL electric propulsion system; specifically:

[0112] eVTOL propulsion systems have stringent performance requirements. Among various traditional electric motors, permanent magnet synchronous motors (PMSMs) have become the mainstream choice for eVTOL propulsion systems due to their significant characteristics of high power density and low torque ripple, and are widely used in engineering fields. This invention selects a surface-mounted PMSM as the research object, and the following is based on... Figure 2 Based on the circuit diagrams of the motor on the d-axis and q-axis shown in (a) and (b) in Table 1, and the motor parameters in Table 1, the following motor dynamic equations are established.

[0113] (27)

[0114] In the formula, This refers to the electromagnetic driving torque of a permanent magnet synchronous motor. For leakage sensation, For magnetizing inductance, and They are permanent magnet synchronous motors Stator inductance along the d-axis and q-axis in the coordinate system. and These are the voltages along the d-axis and q-axis, respectively. and These are the d-axis and q-axis currents, respectively. and These are the iron loss current components along the d-axis and q-axis, respectively. and These are the stator resistance and the core loss resistance, respectively. This represents the number of pole pairs in the motor stator winding. ω is the rotor mechanical angular velocity. This represents the rotor's rotational angle position. It is the excitation flux of the rotor permanent magnet. It is the moment of inertia of the motor rotor. This is the viscosity coefficient of the motor bearing. This refers to the mechanical load fed back to the motor by the driven mechanism.

[0115] Table 1: Motor Parameters

[0116]

[0117] It should be emphasized that the electromechanical coupling modeling method proposed in this invention is not limited to the specific type of surface-mounted permanent magnet synchronous motor. Its core logic and modeling method are universally applicable to other motors suitable for eVTOL propulsion systems, such as asynchronous motors and brushless DC motors. To improve the universality of the model and simplify the complexity of the expression in the subsequent analysis process, this invention will further construct a general dynamic model of the motor. In essence, the motor dynamic equation of the general dynamic model includes a circuit module and a mechanical motion module. The two achieve the coupling and transfer of energy and signals through electromagnetic torque. Its general expression is shown in formula (2).

[0118] Step 2: Based on the eVTOL's geometric configuration and inertia parameters, establish the eVTOL's fuselage multibody system dynamic equations; specifically:

[0119] This embodiment selects, as follows: Figure 3 The variable-speed direct-drive quadrotor eVTOL shown in (a) is used for illustration. However, it should be noted that the modeling method used is also applicable to other configurations. The eVTOL's multi-body system is decoupled into seven independent components, specifically including the main fuselage, left nacelle (containing the left motor), left rotor, right nacelle (containing the right motor), right rotor, upper tail rotor, and lower tail rotor (each corresponding to a built-in motor). For example... Figure 3 As shown in (b), to accurately describe the spatial motion state of each rigid body, a right-handed coordinate system is established with the center of mass of each rigid body as the origin. Meanwhile, to simplify the modeling process while ensuring analytical accuracy, the following assumptions are made regarding the flight characteristics:

[0120] (1) The Earth is considered an inertial frame of reference;

[0121] (2) All components of the fuselage multibody system are considered as rigid bodies, and their elastic deformation is not considered;

[0122] (3) The positions of the centers of mass of each rigid body do not change during flight;

[0123] (4) The fuselage multibody system is symmetrically distributed along the longitudinal plane of symmetry;

[0124] (5) This invention only considers the vertical takeoff phase, and therefore does not consider the aerodynamic characteristics of the wing or the aerodynamic coupling between the rotor and the wing.

[0125] (6) This invention does not consider the movement of the air itself, and the speed of the fuselage in the global coordinate system is equal to the speed of the incoming airflow.

[0126] The first in the fuselage multibody system The position of each component in space is determined by the position vector of its centroid. This indicates that the position vector represents the coordinates of the component's center of mass in the global inertial coordinate system OXYZ. Furthermore, this invention selects Cardan angles. The rotational coordinates, used to describe the rigid body's body-following coordinate system relative to the global inertial coordinate system, are used as the generalized coordinates to describe the position and orientation of the i-th component in space. ,in Indicates the first The coordinates of the centroid position of each component. Indicates the first The rotation coordinates of each component. By combining the generalized coordinates of all components of the fuselage multibody system, a generalized coordinate vector describing the entire fuselage is obtained, as shown in formula (3).

[0127] After obtaining the generalized coordinate vector of the entire fuselage multibody system, the connection relationships between the various rigid body components are explained next. The four rotors at the front and rear are all connected to their supporting structures by rotating joints, while the left and right nacelles are connected to the main fuselage body by fixed supports. After defining the generalized coordinates and connection methods, the dynamic equations of the fuselage multibody system are obtained using the Cartesian modeling method, as shown in formula (4). The mass, moment of inertia, initial position of the center of mass, and initial values ​​of the rotational coordinates of each fuselage component are shown in Table 2.

[0128] Table 2: Parameters and Initial Simulation Values ​​of the Airframe

[0129]

[0130] Step 3: Apply an aerodynamic model that is both easy to calculate and possesses considerable accuracy to the fuselage multibody system dynamic equations established in Step 2; specifically:

[0131] Since eVTOL differs from traditional helicopters in blade aspect ratio, aerodynamic modeling can be appropriately simplified. The calculation formulas for blade thrust and counter-torque are shown in formulas (5) to (11), and the specific values ​​of each parameter are shown in Table 3.

[0132] Table 3: Aerodynamic Parameters

[0133]

[0134] Step 4: Couple the dynamic models from Steps 1 to 3 with the motion driving equations and a general mechanical load calculation scheme based on Lagrange multipliers.

[0135] Step 4-1: Using motion drive equations, describe the control commands of the motor established in Step 1 to the fuselage multibody system established in Step 2. Specifically:

[0136] Motion drive, as the interface for converting motor control commands into the motion of the multibody system, is crucial to the accuracy of the system's dynamic response. This invention employs algebraic constraint equations to achieve precise transmission of motion drive, offering strong engineering scalability. A standardized constraint library can be established based on multibody dynamics constraint classification. For different eVTOL configurations such as tiltrotor and compound airfoil, modeling can be completed simply by calling the corresponding constraint module, significantly reducing the repetitive development costs during configuration iteration.

[0137] like Figure 4 As shown, Let be three mutually perpendicular body unit vectors on the rotor. Let be three mutually perpendicular body unit vectors on the motor. and Both are located on the rotation axes of the motor rotor and the rotor. The motion commands from the motor to the rotor can be realized in the form of motion drive equations, which can be expressed by algebraic constraint equations as follows:

[0138] (28)

[0139] in It is the angle variable in equation (27).

[0140] Step 4.2: Using a general mechanical load calculation scheme based on Lagrange multipliers, calculate the mechanical load on the motor established in Step 1 during the operation of the fuselage multibody system established in Step 2. Specifically:

[0141] During the process of the motor providing motion control commands to the eVTOL rotor, it is subjected to a reaction mechanical load from the drive mechanism, and the two together constitute a... Figure 5 The dynamically coupled mechanical system shown is a dynamic system whose dynamic behavior is described by the fifth equation in formula (27). From the perspective of the load composition mechanism, the physical source of the mechanical load on the motor in this scenario is clear and can be directly disassembled, mainly consisting of the rotor inertial torque and the counter-torque. As shown in formula (12).

[0142] The rotor inertial torque is directly obtained from the rotor's moment of inertia and angular velocity around the axis of rotation, while the aerodynamic counter-torque is related to the rotor's aerodynamic characteristics. A clear quantitative model is established based on rotor aerodynamic theory or wind tunnel test data as shown in Equation (10). Therefore, for the specific scenario of eVTOL, the mechanical load on the motor can be solved by calculating the rotor inertial torque and aerodynamic counter-torque separately and then performing vector superposition. The overall calculation process has a clear physical correspondence.

[0143] However, load calculation in general electromechanical systems such as industrial robotic arms and multi-joint transmission systems faces complexity: load sources include inertial forces, joint coupling forces, and operational loads, and the model needs to be re-decomposed and re-modeled when system parameters or operating conditions change. If the traditional approach of "source analysis - individual calculation - vector superposition" is used, it not only lacks versatility but also easily introduces errors in multi-source load coupling scenarios, limiting modeling efficiency and accuracy. To address this, this invention proposes a general mechanical load calculation method based on Lagrange multipliers. This method can obtain the mechanical load under various operating conditions without decomposing the physical nature of the load, avoiding redundant derivations and handling multi-source coupling problems, thus improving the versatility and adaptability of electromechanical coupling modeling. The general mechanical load calculation method is as follows:

[0144] Assume a rigid body Through constraints Connected to another rigid body, the specific hinge point is... The constraint equations are The Lagrange multipliers corresponding to the constraints are Acting on a rigid body The generalized constraint on is The virtual power created by the generalized constraint force is shown in formula (13).

[0145] At the hinge point Established at the point of being fixed to a rigid body Hinge point coordinate system In the specific analysis process, the coordinate axes of the hinge point coordinate system are defined on the rotor's rotation axis vector. Above, the constraint reactions or driving forces at the hinge points will be described in the hinge point coordinate system. The hinge point coordinate system and the rigid body... Rigid body body coordinate system Directional cosine matrix It is a constant matrix. The constraint forces in the hinge point coordinate system are as shown in formula (14).

[0146] Project the constraint forces onto the global coordinate system, as shown in formulas (15) to (20).

[0147] By comparing equations (13) and (16), the array of constraint forces in the hinge coordinate system can be obtained, which is the rotor mechanical load on the motor rotor in the motor coordinate system.

[0148] Step 5: Based on the coupling scheme in Step 4, establish the eVTOL electromechanical coupling dynamic model, then discretize it using a trapezoidal-Newmark hybrid discretization scheme, and finally solve it using the Newton-Raphson iteration method; specifically:

[0149] Step 5.1: Establish the eVTOL electromechanical coupling dynamic model;

[0150] By combining formula (2) in step 1 and formula (4) in step 2, we obtain the eVTOL electromechanical coupling dynamic model, as shown in formula (21).

[0151] Step 5.2: Discretize the eVTOL electromechanical coupling dynamics model obtained in Step 5.1 using a trapezoidal-Newmark hybrid discretization scheme;

[0152] The eVTOL electromechanical coupling dynamics model constructed in this invention is ultimately presented in the form of a DAE, which includes the differential variables of the first-order differential equations of the motor electrical module, the second-order differential variables of the motor mechanical module, the second-order dynamic differential variables of the fuselage multibody system, and Lagrange multipliers related to the algebraic constraint equations. Traditional numerical discretization schemes are difficult to apply directly.

[0153] To effectively simulate such equations numerically, this invention constructs a trapezoidal-Newmark hybrid discretization scheme. By specifically adapting to the dynamic characteristics of different modules, it achieves efficient and accurate solution of the DAE system. The core idea of ​​this method is: to discretize the first-order differential variables of the motor electrical module using the implicit trapezoidal method, and to discretize the second-order differential variables of the motor mechanical module, the fuselage multibody system, and the Lagrange multipliers related to the algebraic constraint equations using the Newmark method, as shown in formula (22). In formula (22) of this embodiment, the simulation time is 8s, and the calculation step size is set to seconds, the calculation parameters for the Newmark algorithm are selected as follows: and .

[0154] It is particularly important to note that the discretization of the motor mechanical variables and the generalized coordinates of the fuselage must use the same set of parameters; otherwise, the calculation deviation of the mechanical motion response will lead to excessively large residuals in the algebraic constraint equations, which will disrupt the motion coupling consistency between the motor and the fuselage and ultimately make convergence difficult.

[0155] By employing the discretization strategy shown in Equation (22), the continuous-time domain eVTOL electromechanical coupling dynamics model is transformed into a system of nonlinear algebraic equations in the discrete-time domain. As the basic unknowns, equation (21) is discretized into the following form:

[0156] (29)

[0157] in, and These correspond to the residuals of the electrical and mechanical modules of the motor, respectively. For the multibody dynamics equations residuals of the fuselage, For the algebraic constraint equation residuals, and for algorithm stability, Multiply by the constraint equation.

[0158] Step 5.3: The Newton-Raphson iteration method is used to solve the discrete form of the eVTOL electromechanical coupling dynamics model obtained in Step 5.2;

[0159] The Newton-Raphson iteration method is used to solve the nonlinear algebraic equation system shown in formula (23) at each time step. The initial values ​​of all variables are 0. The linear equation system to be solved at each iteration step is given below, as shown in formulas (24) to (26).

[0160] Figure 6 A flowchart of the computational process based on the above method is presented. First, the simulation duration, step size, and algorithm parameters are set. Then, the state variables of the motor and fuselage are initialized. Next, the basic variables to be solved during the iteration process are initialized. Then, the large time step loop is entered, followed by the small iteration loop. Within the iteration loop, a discretization scheme is applied, and then the residual vector is calculated. If the magnitude of the residual vector is less than a set value, the variable values ​​are updated and the process proceeds to the next time step. If the magnitude of the residual vector is greater than or equal to the set value, the iteration continues to update the variables.

[0161] The voltage input curve of this embodiment is as follows: Figure 7 As shown, Represent Figure 3 In example (a), the input voltage of the drive motors corresponding to rotors 3, 5, 6, and 7 is adjusted as follows: a small voltage is applied during the first 0-2 seconds to keep the eVTOL hovering; then, the voltage is increased simultaneously on all four motors during the second second, causing the eVTOL to begin climbing. The altitude time-history curve of the eVTOL obtained in this example is shown below. Figure 8 As shown. In this example, Figure 3 The energy consumption distribution of the motor corresponding to rotor No. 3 in (a) is as follows: Figure 9 As shown.

[0162] The above-described embodiments are merely illustrative of the implementation methods of the present invention, but should not be construed as limiting the scope of the present invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the protection scope of the present invention.

Claims

1. A method for modeling the electromechanical coupling dynamics of an electric vertical takeoff and landing aircraft based on differential algebraic equations, characterized in that, The electromechanical coupling dynamics modeling method for electric vertical takeoff and landing aircraft includes the following steps: Step 1: Based on the electrical and mechanical motion parameters of the motor, establish the motor dynamics equations in the eVTOL electric propulsion system; Step 2: Based on the geometric configuration and inertia parameters of the eVTOL, establish the dynamic equations of the eVTOL fuselage multibody system; Step 3: Apply an aerodynamic model that is easy to calculate and has considerable accuracy to the fuselage multibody system dynamic equations established in Step 2; Step 3 specifically includes: The formulas for calculating the blade thrust and counter-torque are as follows: （5）； （6）； in, For the blade tension; For reverse torque; The blade rotation speed; The diameter of the blade; and These are the dimensionless tensile force coefficient and torque coefficient, respectively; It refers to the air density in the flight environment; The calculation formulas for blade thrust and counter-torque of the eVTOL under vertical takeoff and landing conditions are as follows: （9）； （10）； in: （11）； The parameters in equations (9)-(11) are defined as follows: This is the correction factor for the effective area of ​​the blade; It is the number of blades; It is the lift slope coefficient of the airfoil; It is a correction factor resulting from the downwashing effect; It is the blade pitch; It is the diameter of the blade; It is zero-lift angle of attack; It is the aspect ratio; It is the average radius proportionality coefficient, used to estimate the average incoming flow velocity on the blade; It is the airfoil drag coefficient of the blade section; It is the zero-lift drag coefficient; It is the Oswald factor; It is the aircraft velocity of eVTOL in the vertical direction; Step 4: Couple the dynamic models of Steps 1 to 3 with the motion driving equations and the general mechanical load calculation scheme based on Lagrange multipliers; Step 5: Based on the coupling scheme in Step 4, establish the eVTOL electromechanical coupling dynamic model, then discretize it using a trapezoidal-Newmark hybrid discretization scheme, and finally solve it using the Newton-Raphson iteration method.

2. The electromechanical coupling dynamics modeling method for an electric vertical takeoff and landing aircraft based on differential algebraic equations according to claim 1, characterized in that, Step 1 specifically includes: A surface-mounted permanent magnet synchronous motor was selected as the research object, and the dynamic equation of the surface-mounted permanent magnet synchronous motor shown in formula (1) was established. （1）； In the formula, The electromagnetic driving torque of a permanent magnet synchronous motor. This represents the number of pole pairs in the stator winding of the motor. This is the excitation flux of the rotor permanent magnet; For leakage sensation, For magnetizing inductance, and They are permanent magnet synchronous motors Stator inductance along the d-axis and q-axis in the coordinate system; and These are the voltages along the d-axis and q-axis, respectively. and These represent the d-axis and q-axis currents, respectively; t represents time. and These are the d-axis and q-axis iron loss current components, respectively; and These are the stator resistance and the core loss resistance, respectively. This refers to the rotor's mechanical angular velocity; This refers to the rotor's rotational angle position; It is the moment of inertia of the motor rotor; The viscosity coefficient of the motor bearing; The mechanical load fed back to the motor by the driven mechanism; A general dynamic model of the motor is constructed, and its general expression is written as: （2）； in, Electrical variables representing the motor Represents the mechanical motion variables of the motor. and These represent the first and second derivatives of the variable with respect to time, respectively; This represents the input voltage of the motor. This represents the mechanical load transmitted to the motor by the multibody system of the fuselage; and These are the nonlinear mapping functions of the electric motor's electrodynamics and mechanical dynamics, respectively, and their specific forms are determined by the motor type and structural parameters.

3. The electromechanical coupling dynamics modeling method for an electric vertical takeoff and landing aircraft based on differential algebraic equations according to claim 2, characterized in that, Step 2 specifically includes: The eVTOL multibody system is decoupled into seven independent components, including the main fuselage, left nacelle, left rotor, right nacelle, right rotor, upper tail rotor, and lower tail rotor; and the following assumptions are made: (1) The Earth is considered an inertial frame of reference; (2) All components of the fuselage multibody system are considered as rigid bodies, and their elastic deformation is not considered; (3) The positions of the centers of mass of each rigid body do not change during flight; (4) The fuselage multibody system is symmetrically distributed along the longitudinal plane of symmetry; (5) Only the vertical takeoff phase is considered, and therefore the aerodynamic characteristics of the wing and the aerodynamic coupling between the rotor and the wing are not considered; (6) Without considering the motion of the air itself, the speed of the fuselage in the global coordinate system is equal to the speed of the incoming airflow; The first in the fuselage multibody system The position of each component in space is determined by the position vector of its centroid. This indicates that the position vector represents the coordinates of the component's center of mass in the global inertial coordinate system OXYZ; this choice uses Cardan angles. The position vector plus the rotation coordinates are used to describe the rotation of the rigid body's body-following coordinate system relative to the global inertial coordinate system; the position vector plus the rotation coordinates are used to describe the first... Generalized coordinates of the position and orientation of each component in space ,in Indicates the first The coordinates of the centroid position of each component. Indicates the first The rotational coordinates of each component; by combining the generalized coordinates of all components of the fuselage multibody system, a generalized coordinate vector describing the entire fuselage is obtained: （3）； in, This represents the generalized coordinate vector of the fuselage multibody system. This represents the generalized coordinate vector of component number 1. This represents the generalized coordinate vector of component number 2. This represents the generalized coordinate vector of component number 7; The dynamic equations of the fuselage multibody system are obtained by using the Cartesian modeling method, as shown in equation (4); （4）； in, , , These represent the generalized coordinates, generalized velocity, and generalized acceleration of the fuselage multibody system, respectively. It is the total mass matrix of the fuselage multibody system, which is obtained by combining the mass matrices of each component. It is a Lagrange multiplier; It is time; These are algebraic constraint equations, representing the connection relationships between the various components; It is the Jacobian matrix of the constraint equations with respect to the generalized coordinates; It refers to the generalized external forces acting on the fuselage multibody system.

4. The electromechanical coupling dynamics modeling method for an electric vertical takeoff and landing aircraft based on differential algebraic equations according to claim 1, characterized in that, Step 4 specifically includes: Step 4.1: Use motion drive equations to describe the control commands of the motor established in Step 1 to the fuselage multibody system established in Step 2. Step 4.2: Using a general mechanical load calculation scheme based on Lagrange multipliers, calculate the mechanical load on the motor established in Step 1 during the operation of the fuselage multibody system established in Step 2. Specifically: During the process of the motor providing motion control commands to the eVTOL rotor, it is subjected to a reaction mechanical load from the drive mechanism, and the two together constitute a dynamically coupled mechanical system; the dynamic behavior of the dynamically coupled mechanical system is described by the fifth equation in formula (1), which consists of two parts: rotor inertial torque and counter-torque; that is: （12）； in, It is the moment of inertia of the rotor about its rotation axis. It is the reverse torque; The mechanical load on the motor is solved by calculating the rotor inertial torque and aerodynamic counter-torque separately and then performing vector superposition. A general mechanical load calculation method based on Lagrange multipliers is used to obtain the mechanical load under various working conditions.

5. The electromechanical coupling dynamics modeling method for an electric vertical takeoff and landing aircraft based on differential algebraic equations according to claim 4, characterized in that, In step 4.2, the general calculation method for mechanical load is as follows: Assume a rigid body Through constraints Connected to another rigid body, the specific hinge point is... The constraint equations are The Lagrange multipliers corresponding to the constraints are Acting on a rigid body The generalized constraint on is The virtual power represented by the generalized constraint force is: （13）； In the formula, It is a variational symbol; It is a rigid body Translation coordinates; It is a rigid body The rotation coordinates; It is virtual power; It is a rigid body Generalized velocity; It is a constraint The corresponding constraint equations for rigid bodies The Jacobian matrix of the generalized coordinates; It is a rigid body The first derivative of the translational coordinates with respect to time; It is a rigid body The first derivative of the rotating coordinates with respect to time; At the hinge point Established at the point of being fixed to a rigid body Hinge point coordinate system Define the coordinate axes of the hinge point coordinate system on the rotor's rotation axis vector. Above; hinge point coordinate system and rigid body Rigid body body coordinate system Directional cosine matrix The constraint force is a constant matrix; the constraint force is denoted as follows in the hinge point coordinate system: （14）； In the formula, and Constraints The array of constraint forces and constraint couples in the hinge point coordinate system; Representing a rigid body The projection of the constraint forces onto the hinge point coordinate system; Projecting the constraint forces onto the global coordinate system yields: （15）； in, Representing a rigid body The projection of the constraint forces onto the global coordinate system; Representing a rigid body Rotation matrix of the body coordinate system; Indicates that it is fixed to a rigid body Hinge point coordinate system rotation matrix; Representing a rigid body The projection of the constraint forces onto the hinge point coordinate system; The virtual power of the constraint force is expressed as: （16）； In the formula, Representing a rigid body upper hinge point The generalized velocity matrix; definition ,and hinge point The projection of the velocity onto the global coordinate system: （17）； in, It is the origin of the body coordinate system The first derivative of the position vector with respect to time; It is a hinge point Relative to the origin of the body coordinate system The radius vector, Indicates the symbol for an antisymmetric matrix; rigid body Given the angular velocity vector, then: （18）； in, It is a rigid body The transformation matrix for converting the derivative of the rotating coordinates to angular velocity depends on the chosen rotating coordinates; It is a rigid body The rotation coordinates; From equations (17) and (18), we obtain: （19）； In the formula, Representing a rigid body Generalized velocity; Represents the generalized velocity transformation matrix; Will Transform into a rigid body The transformation matrix for generalized velocity is as follows: （20）； in, Represents the identity matrix; This indicates a hinge point. Relative to the origin of the body coordinate system The antisymmetric matrix formed by the position vectors; Representing a rigid body The transformation matrix for converting the derivative of rotating coordinates to angular velocity is then... Representing a rigid body The inverse of the transformation matrix for converting the derivative of rotating coordinates to angular velocity; By comparing equations (13) and (16), the constraining force array in the hinge coordinate system is obtained, which is the rotor mechanical load on the motor rotor in the motor coordinate system.

6. The electromechanical coupling dynamics modeling method for an electric vertical takeoff and landing aircraft based on differential algebraic equations according to claim 5, characterized in that, Step 5 specifically involves: Step 5.1: Establish the eVTOL electromechanical coupling dynamic model; Combining formula (2) from step 1 and formula (4) from step 2, we obtain the eVTOL electromechanical coupling dynamic model, as shown below: （21）； in, This represents the first derivative of the electrical variable of the motor with respect to time. This represents the second derivative of the mechanical motion variables of the motor with respect to time. Represents the electrical variables of the motor; The nonlinear mapping function representing the electrodynamics of the motor; This indicates the mechanical load transmitted to the motor by the fuselage multibody system; Indicates time; A nonlinear mapping function representing the mechanical dynamics of an electric motor; This represents the total mass matrix of the fuselage multibody system; This represents the generalized acceleration of the fuselage multibody system; It represents the Jacobian matrix of the constraint equations with respect to the generalized coordinates; This indicates that it is a Lagrange multiplier; This refers to the generalized external forces acting on the fuselage multibody system; Represent the constraint equations; Represents the mechanical motion variables of the motor; Represents the generalized coordinates of the fuselage multibody system; Step 5.2: Discretize the eVTOL electromechanical coupling dynamics model obtained in Step 5.1; The eVTOL electromechanical coupling dynamics model includes the differential variables of the first-order differential equations of the motor electrical module, the second-order differential variables of the motor mechanical module, the second-order dynamic differential variables of the fuselage multibody system, and Lagrange multipliers related to the algebraic constraint equations. A trapezoidal-Newmark hybrid discretization scheme is used: the first-order differential variables of the motor electrical module are discretized using the implicit trapezoidal method, while the second-order differential variables of the motor mechanical module, the fuselage multibody system, and the Lagrange multipliers related to the algebraic constraint equations are discretized using the Newmark method. This transforms the continuous-time domain eVTOL electromechanical coupling dynamics model into a system of nonlinear algebraic equations in the discrete-time domain. Step 5.3: Using the Newton-Raphson iteration method, solve the system of nonlinear algebraic equations obtained in step 5.2 at each time step.

7. The electromechanical coupling dynamics modeling method for an electric vertical takeoff and landing aircraft based on differential algebraic equations according to claim 6, characterized in that, In step 5.2, the trapezoidal-Newmark hybrid discretization scheme is used for discretization as shown in formula (22): （22）； Among them, subscript and Represents the start and end nodes of a discrete time period; Represents the length of a discrete interval; and It is a parameter in the Newmark algorithm that controls the stability and accuracy of the algorithm; Electrical variables representing the motor at the end of the time period; Electrical variables representing the motor at the start point of the time period; This represents the first derivative of the electrical variable of the motor at the end of the time interval with respect to time. The first derivative of the electrical variable of the motor with respect to time at the starting point of the time period; The mechanical motion variable of the motor at the end of the time period; The mechanical motion variables of the motor at the start point of the time period; The first derivative of the mechanical motion variable of the motor with respect to time at the starting point of the time interval; This represents the first derivative of the mechanical motion variable of the motor at the end of the time interval with respect to time. The second derivative of the mechanical motion variable of the motor with respect to time at the starting point of the time interval; This represents the second derivative of the mechanical motion variable of the motor with respect to time at the end of the time interval. The generalized coordinates of the fuselage multibody system representing the end node of the time period; The generalized coordinates of the fuselage multibody system representing the starting point of a time period; The generalized velocity of the fuselage multibody system at the start point of the time period; The generalized acceleration of the fuselage multibody system represents the starting point of the time period; The generalized acceleration of the fuselage multibody system at the end of the time period; This represents the generalized velocity of the fuselage multibody system at the end of the time period; Through the discretization strategy shown in Equation (22), the eVTOL electromechanical coupling dynamics model is transformed into a system of nonlinear algebraic equations, in order to As the basic unknowns, equation (21) is discretized into the following form: （23）； in, and These correspond to the residuals of the electrical and mechanical modules of the motor, respectively. For the multibody dynamics equations residuals of the fuselage, For the algebraic constraint equation residuals, and for algorithm stability, Multiply by the constraint equation.

8. The electromechanical coupling dynamics modeling method for an electric vertical takeoff and landing aircraft based on differential algebraic equations according to claim 7, characterized in that, Step 5.3 specifically refers to: The system of linear equations that needs to be solved in each iteration step is as follows: （24）； Write it in explicit format: （25）； in, express right The Jacobian matrix; express right The Jacobian matrix; express right The Jacobian matrix; express right The Jacobian matrix; express right The Jacobian matrix; express right The Jacobian matrix; express right The Jacobian matrix; express right The transpose of the Jacobian matrix; express right The Jacobian matrix; express right The Jacobian matrix; The iterative increment of the electrical variable of the motor at the end of the time period; The iterative increment of the mechanical motion variable of the motor at the end of the time period; This represents the iterative increment of the generalized acceleration of the fuselage multibody system at the end of the time period; This represents the Lagrange multiplier iteration increment at the end of the time period; but: (26) ； in, Indicates the previous iteration step; This indicates the next iteration step.

Images