Decoupling multi-degree-of-freedom dynamics simulation method for automobile chassis system

Abstract

The present application belongs to the technical field of automobile chassis control system development and simulation verification, and particularly relates to a decoupling multi-degree-of-freedom dynamics simulation method of an automobile chassis system, which can realize the migration of robot whole body control to an automobile chassis. The present application splits the vehicle degrees of freedom into five independent integral channels and adopts a two-way transmission mechanism of suspension coupling force based on the least square method, while automatically calculating the multi-body dynamics library, fundamentally eliminating the non-physical cross coupling introduced by the unified integration of the traditional articulated body algorithm, and realizing the real-time simulation of vehicle dynamics with both physical correctness and calculation efficiency. The statistical comparison verification framework automatically runs the WBC control, the traditional control and the passive mode on the standard scene set covering multiple working conditions, calculates the standardized performance indicators, and introduces an out-of-control automatic detection mechanism including six indicators such as side slip and roll over, and sets the trigger and recovery threshold hysteresis.

Description

Technical Field

[0001] This invention belongs to the field of automotive chassis control system development and simulation verification technology, specifically involving a decoupled multi-degree-of-freedom dynamic simulation method for automotive chassis systems. Background Technology

[0002] Currently, vehicle dynamics simulation technology is the cornerstone of chassis controller development and performance verification. The mainstream vehicle multibody dynamics simulation methods currently fall into three main categories: Commercial multibody dynamics software, such as VI-CarRealTime, CarMaker, and CarSim, typically builds complete vehicle multibody dynamics models based on Lagrange equations or the Newton-Euler method. Their advantage is high model accuracy. However, such software has the following drawbacks: the source code is closed-source, preventing users from deeply modifying the physical model to meet specific needs; the simulation step size is usually limited to 0.1~1ms, resulting in high single-step calculation time and making it difficult to embed real-time control loops for hardware-in-the-loop (HIL) testing; furthermore, commercial licensing costs are also high.

[0003] Rigid body multibody simulation based on Pinocchio (ABA algorithm) utilizes the articulated body algorithm (ABA) in the open-source multibody dynamics library Pinocchio for forward dynamics simulation. This method abstracts the vehicle into a tree-like kinematic chain: from the floating base (6 degrees of freedom) to the suspension joints (4 prism pairs) to the steering joints (2 revolute pairs) to the wheels (4 consecutive revolute pairs). The generalized acceleration is automatically calculated by Pinocchio using a unified mass matrix and Coriolis force matrix.

[0004] However, the ABA algorithm has inherent limitations in automotive chassis simulation scenarios. ABA obtains acceleration by solving for the inverse of the generalized mass matrix of the entire system. When the rotational degrees of freedom of the wheels are included in the solution of a unified generalized mass matrix, non-physical cross-coupling effects are introduced. In actual physics, the rotational motion of the wheels is decoupled from the vertical motion of the vehicle body / suspension through the tire-ground contact force; there is no direct inertial coupling between the two. However, the unified solution by ABA mathematically results in: the vertical vibration of the suspension being coupled to the rotational acceleration of the wheels through the mass matrix, and the driving force / braking torque of the wheels being coupled to the vertical acceleration of the suspension through the mass matrix. These artifacts are caused by the ABA algorithm inverting all degrees of freedom under the same mass matrix, which reduces the simulation credibility under certain conditions.

[0005] The 14-DOF analytical model, a classic 14-DOF vehicle dynamics model, decomposes the vehicle into independent integral subsystems, including: body center of mass translation (3DOF), body rotation (pitch, roll, yaw, 3DOF), vertical compression of the four suspensions (4DOF), and rotation of the four wheels (4DOF). The coupling between these subsystems is transmitted through explicit mechanical relationships (such as suspension forces and tire forces), rather than through the inversion of a unified mass matrix. The advantages of this approach are clear physical relationships and high computational efficiency. However, the disadvantages are that all coupling relationships must be manually derived and modeled, and the ability to automatically generate Jacobian and mass matrices using existing multibody dynamics libraries cannot be utilized, resulting in higher system complexity and a higher modeling error rate.

[0006] In summary, existing vehicle dynamics simulation methods suffer from contradictions in terms of physical accuracy, real-time performance, customizability, and modeling efficiency. In particular, there is a lack of a hybrid simulation scheme that can both utilize the automated computation capabilities of multibody dynamics libraries and physically avoid the non-physical coupling introduced by the ABA unified integral. Summary of the Invention

[0007] In view of this, the present invention provides a decoupled multi-degree-of-freedom dynamics simulation method for automobile chassis systems, which can fundamentally eliminate the non-physical cross-coupling introduced by the unified integration of traditional articulated body algorithms, and realize real-time vehicle dynamics simulation that combines physical correctness and computational efficiency.

[0008] To achieve the objectives of this invention, the following technical solutions are provided.

[0009] A decoupled multi-degree-of-freedom dynamics simulation method for an automotive chassis system includes: splitting the 16 degrees of freedom of the vehicle into five independent integration channels, namely: a 3-degree-of-freedom attitude channel, a 4-degree-of-freedom suspension channel, a 4-degree-of-freedom wheel channel, a 2-degree-of-freedom steering channel, and a 3-degree-of-freedom center-of-mass channel; the attitude channel uses the Euler rigid body equation with inertia product Ixz and the quaternion integration method to independently integrate the pitch angle, roll angle, and yaw angle of the vehicle body; the suspension channel uses a spring-damping model to independently integrate the compression of the four suspension components; the wheel channel uses a tire slip ratio model to independently integrate the angular velocities of the four wheels; calculating the suspension coupling force generated by the pitch moment and roll moment, and simultaneously transferring the suspension coupling force to the suspension channel and the attitude channel to form a two-way physical feedback closed loop; and performing parallel integration on the five independent integration channels within the same simulation step, wherein each channel uses an integration method that matches its physical characteristics.

[0010] In the attitude channel integration process, the synthesized total external torques τx, τy, and τz are calculated using the following formulas, directly incorporating the roll torque Mx and pitch torque My from the suspension system: τx = Mx + (Iyy - Izz) * ωy * ωz, τy = My + (Izz - Ixx) * ωz * ωx, τz = Mz + (Ixx - Iyy) * ωx * ωy, where ωx, ωy, and ωz are the roll, pitch, and yaw angular velocities, respectively; Ixx, Iyy, and Izz are the principal inertia; and Ixz is the inertia product. The above equations are expressed in matrix form as the inertia matrix I multiplied by the angular acceleration vector equals the vector on the right. The angular accelerations of each axis are obtained by solving the 2×2 coupled subsystem (dwx and dwz are cross-coupled through Ixz, and dwy is solved independently) using Cramer's rule.

[0011] The dynamic equation of the spring-damped model of the suspension channel is: mw * dvsusp / dt = Fz_tire + (-ks * (z - z0) - cs * vsusp) + Fcoupling - mw * g, where mw is the unsprung mass (wheel-end mass), vsusp is the suspension compression velocity, ks is the suspension stiffness, cs is the suspension damping coefficient, z is the current compression, z0 is the static equilibrium position, and Fcoupling is the force from the coupling force channel.

[0012] The dynamic equation of the wheel channel is: Iw * dω / dt = Tdrive - Fx * R, where Iw is the wheel moment of inertia, ω is the wheel angular velocity, Tdrive is the driving torque from the controller, Fx is the tire-ground longitudinal force, and R is the effective radius of the wheel. The tire slip ratio model is: slip ratio s = (ω * R - vx) / max(|vx|,v_ref), where v_ref is the low-speed regularized reference speed, longitudinal force Fx = ftire(s, Fz), where ftire is the tire force model, supporting linear and Pacejka Magic Formula models, where vx is the wheel center longitudinal velocity, and Fz is the tire vertical force.

[0013] In step (e), under braking conditions, the pitching moment My is generated by the braking deceleration ax. The suspension coupling force distribution formula is as follows: The coupling force is calculated from the external moment using the least squares method to ensure zero net force constraint: For the pitching moment My, the coupling force of the i-th suspension is ΔFz_i = (My / Sxx) * (x_i - ), where Sxx is the sum of squared deviations of the longitudinal coordinates of the suspension, and x_i is the i-th longitudinal coordinate of the suspension. The mean is used; this method guarantees ΣΔFz_i=0 (pure couple) and ΣΔFz_i*x_i=My; for the roll moment Mx, the coupling force of the i-th suspension is ΔFz_i += (Mx / Syy) * (y_i - ), where Syy is the sum of squared deviations of the suspension lateral coordinates, and y_i is the i-th suspension lateral coordinate. The value is the average; the external torque is generated by acceleration: during braking, My = m*ax*h_cg, and during tilting, Mx = m*ay*h_cg, where h_cg is the height of the center of mass.

[0014] The method also integrates the multibody dynamics library Pinocchio, but the Pinocchio library is not used to perform forward simulation of the articulated body algorithm ABA, but only to perform the following operations: forward kinematics calculation to determine the wheel contact point position and suspension joint angle; Jacobian matrix calculation to obtain the Jacobian matrix of the stability domain; and mass matrix query to provide inertial parameters.

[0015] The simulation step size is 0.005 seconds, corresponding to a simulation frequency of 200Hz; the attitude channel, suspension channel, wheel channel, and center of gravity channel are integrated using a semi-implicit Euler method; the steering channel uses a position control method, directly writing the target steering angle into the simulation state vector without using differential equation integration.

[0016] The process also includes a statistical comparison and verification step: on a scenario set covering multiple standardized test conditions, the whole-body control (WBC) mode, traditional control mode, and passive mode are automatically run in parallel; standardized performance indicators are calculated, including the load transfer ratio (LTR); and an automatic runaway detection mechanism is introduced, which includes six runaway indicators, including sideslip and rollover indicators, and sets trigger and recovery thresholds for each indicator to form hysteresis protection.

[0017] Beneficial effects 1. This invention enables the transfer of full-body robot control to an automobile chassis.

[0018] 2. This invention divides the vehicle's degrees of freedom into five independent integration channels and employs a bidirectional suspension coupling force transmission mechanism based on the least squares method. While utilizing the automatic calculation capabilities of the multibody dynamics library, it fundamentally eliminates the non-physical cross-coupling introduced by the unified integration of traditional articulated body algorithms, achieving real-time vehicle dynamics simulation that combines physical accuracy and computational efficiency. 3. The statistical comparison and verification framework automatically runs WBC control, traditional control and passive mode on a standard scenario set covering multiple working conditions, calculates standardized performance indicators, and introduces an automatic runaway detection mechanism that includes six indicators such as sideslip and rollover and sets trigger and recovery threshold lag.

[0019] 4. The system integrates the Pinocchio multibody dynamics library, but its usage differs fundamentally from traditional methods. In this invention, the Pinocchio library is not used for performing ABA forward simulations, but only for forward kinematics (FK) calculations to determine wheel contact point locations and suspension joint angles; Jacobian matrix calculations to obtain the stability domain Jacobian matrix; and mass matrix lookups to provide inertial parameters. This fundamentally avoids the non-physical cross-coupling caused by inverting the unified ABA mass matrix.

[0020] 5. The simulation method of this invention divides the 16 degrees of freedom of the vehicle into five independent integration channels: attitude channel (3 degrees of freedom), suspension channel (4 degrees of freedom), wheel channel (4 degrees of freedom), steering channel (2 degrees of freedom), and center of gravity channel (3 degrees of freedom). Each channel operates independently, using the integration method best suited to its physical characteristics. Attached Figure Description

[0021] Figure 1 This is a flowchart of a method according to an embodiment of the present invention. Detailed Implementation

[0022] The present invention will be further described below with reference to the accompanying drawings and specific embodiments.

[0023] This invention proposes a decoupled multi-degree-of-freedom dynamics simulation method for automotive chassis systems, implemented within a "HybridForwardSimulator" system. This system aims to address the non-physical coupling issues inherent in the traditional ABA (Articulated Body Algorithm) unified integration method, achieving physically accurate and computationally efficient simulations.

[0024] Furthermore, the system integrates the Pinocchio multibody dynamics library, but its usage differs fundamentally from traditional methods. In this invention, the Pinocchio library is not used to perform ABA forward simulations, but only for: (1) forward kinematics (FK) calculations to determine wheel contact point locations and suspension joint angles; (2) Jacobian matrix calculations to obtain the stability domain Jacobian matrix; and (3) mass matrix lookups to provide inertial parameters. This fundamentally avoids the non-physical cross-coupling caused by the inversion of the unified ABA mass matrix. Figure 1 This is a flowchart of a method according to an embodiment of the present invention.

[0025] The simulation method of this invention divides the 16 degrees of freedom of the vehicle into five independent integration channels: attitude channel (3 degrees of freedom), suspension channel (4 degrees of freedom), wheel channel (4 degrees of freedom), steering channel (2 degrees of freedom), and center of gravity channel (3 degrees of freedom). Each channel operates independently, using the integration method best suited to its physical characteristics, and interacts through coupling forces. The detailed design of the five independent integration channels is as follows: The attitude channel is used to integrate the three rotational degrees of freedom of the vehicle body: roll (φ), pitch (θ), and yaw (ψ). This channel uses Euler's rigid body equations with inertia product (Ixz) for dynamic integration and employs quaternion methods to avoid singularities.

[0026] The synthesis of external moments is key to this process. The total synthesized external moments (τx, τy, τz) include the direct acting moments and the coupled moments from the suspension system (Mx is the roll moment, and My is the pitch moment), as shown in the following equation: τx = Mx + (Iyy - Izz) * ωy * ωz τy = My + (Izz - Ixx) * ωz * ωx τz = Mz + (Ixx - Iyy) * ωx * ωy Where ωx, ωy, and ωz represent the roll, pitch, and yaw angular velocities, respectively; Ixx, Iyy, and Izz represent the principal inertia; and Ixz represents the inertia product. The above equations, expressed in matrix form, are: the inertia matrix I multiplied by the angular acceleration vector equals the vector on the right (the sum of the gyro terms and the external torque). The inertia matrix I is [Ixx, 0, Ixz; 0, Iyy, 0; Ixz, 0, Izz]. The angular accelerations of each axis are obtained by solving the 2×2 coupled subsystem (dwx and dwz are cross-coupled through Ixz, and dwy is solved independently) using Cramer's rule. This equation directly incorporates the suspension coupling torque into the attitude integral channel, rather than transmitting the coupling force solely through the mass matrix as in the traditional ABA method.

[0027] Attitude integration is performed using quaternions. The vehicle body attitude is represented by the quaternion q = [q0, q1, q2, q3]. This means that the angular velocity ω = [ωx, ωy, ωz] Update as needed. The transformation relationship from angular velocity to quaternion derivative is as follows: dq / dt = 1 / 2 * Ω(ω) * q Where Ω(ω) is a skew-symmetric matrix composed of angular velocity components. The quaternions are updated using a first-order explicit integration method: q(t+dt) = q(t) + dq / dt * dt To preserve the properties of a unit quaternion, normalization is required after integration.

[0028] The suspension channel is used to independently integrate four suspension compression parameters (zFL, zFR, zRL, zRR). Each suspension uses a spring-damped model with the following dynamic equations: mw * dvsusp / dt = Fz_tire + (-ks * (z - z0) - cs * vsusp) + Fcoupling- mw * g Where mw is the unsprung mass (wheel-end mass, e.g., 40 kg), vsusp is the suspension compression velocity, ks is the suspension stiffness, cs is the suspension damping coefficient, z is the current compression, z0 is the static equilibrium position, Fcoupling is the force from the coupling force channel, Fz_tire is the tire normal force, and g is the acceleration due to gravity (9.81 m / s²). This channel uses a semi-implicit Euler method for integration to ensure numerical stability.

[0029] The angular velocities (ωFL, ωFR, ωRL, ωRR) of the four wheels are independently integrated along the wheel path. The dynamic equation is: Iw * dω / dt = Tdrive - Fx * R Where Iw is the wheel's moment of inertia, Tdrive is the driving torque from the controller, Fx is the tire-ground longitudinal force, and R is the effective wheel radius. This channel combines semi-implicit integrals with a tire slip ratio model for calculation: The slip ratio s = (ω * R - vx) / max(|vx|, v_ref), where v_ref is the low-speed regularization reference velocity (e.g., 0.5 m / s) to avoid the denominator approaching zero at low speeds, which could lead to numerical instability. Longitudinal force Fx = ftire(s, Fz) Where vx is the longitudinal velocity at the wheel center and Fz is the vertical force of the tire.

[0030] The steering channel directly controls the steering angles of the two front wheels (δFL, δFR) using a position control method: the target steering angle δcmd is directly written into the simulation state vector without using differential equation integration. This design eliminates the need for a time constant parameter in the steering channel, simplifies the control interface, and allows direct interface with the controller's steering output.

[0031] (Note: The original first-order steering actuator dynamics model has been replaced with the above position control method. The following paragraphs are historical descriptions.) The center-of-mass channel is used to integrate the three translational degrees of freedom (x, y, z) of the vehicle's center of mass, which obey Newton's second law: m * dv / dt = Ftotal Where m is the total vehicle mass (using sprung mass m_chassis for longitudinal and lateral movement, and total mass m_total = m_chassis + 4 × mw for vertical movement), v is the center-of-mass velocity vector, and Ftotal is the resultant force acting on the center of mass (including longitudinal tire force, lateral tire force, air resistance, etc.). This channel uses a semi-implicit Euler method for position and velocity updates.

[0032] The specific braking condition is as follows: the braking deceleration ax generates a pitching moment My = m * ax * h_cg (where m is the vehicle mass and h_cg is the center of gravity height). This moment is distributed to the four suspensions using the least squares method, and the coupling force distribution formula is: ΔFz_i = (My / Sxx) * (x_i - x ), where Sxx = Σ(x_i - x )² represents the sum of squared deviations of the longitudinal coordinates of the suspension, and x_i represents the longitudinal coordinate of the i-th suspension. The mean This method guarantees three physical constraints: ΣΔFz_i=0 (pure couple, no net vertical force generated), ΣΔFz_i*x_i=My (accurately restoring the pitch moment), and the force of each wheel is distributed according to the wheelbase ratio. The coupling force is transmitted to the suspension channel on the one hand, changing the suspension displacement; on the other hand, the pitch moment My generated by it is transmitted to the attitude channel, directly changing the pitch rate of the vehicle body, forming the correct physical feedback for the "nodding" effect.

[0033] The roll condition is as follows: the lateral acceleration ay generates a roll moment Mx = m * ay * h_cg. The coupling force distribution formula is: ΔFz_i += (Mx / Syy) * (y_i - y ), where Syy = Σ(y_i - y )² represents the sum of squared deviations of the lateral coordinates of the suspension, and y_i represents the lateral coordinate of the i-th suspension. The mean This coupling force is also transmitted bidirectionally to the suspension channel and attitude channel, accurately simulating vehicle body roll.

[0034] The simulation integration timing is as follows: within a typical simulation step (e.g., dt = 0.005s, i.e., 200Hz), the decoupled multi-degree-of-freedom hybrid simulator performs the following steps: Step 1: State Reading and Coupling Force Calculation. Read the vehicle state (generalized coordinates q and generalized velocity v) at the current simulation moment. Calculate the coupling force Fcoupling generated by the pitch moment My and roll moment Mx based on the current state.

[0035] Step 2: Ground force calculation. Based on the current vehicle state (such as suspension displacement, wheel speed, and vehicle posture) and road model, calculate the interaction force between each tire and the ground, including normal force Fz, longitudinal force Fx, and lateral force Fy.

[0036] Step 3: Parallel integration across five channels. Using the concept of parallel computation, the state of each of the five channels is updated independently (the steering channel is for position control and is directly assigned a value): The attitude channel updates the angular velocity according to Euler's rigid body equations, and updates the vehicle attitude using the quaternion method based on the angular velocity.

[0037] The suspension channel updates the compression amount and compression speed of the four suspensions based on the spring-damping model and the received coupling force.

[0038] The wheel channel updates the angular velocity of the four wheels based on the driving torque and the longitudinal force of the tires.

[0039] The steering channel updates the steering angle of the two front wheels according to the steering command.

[0040] The center of mass channel updates the velocity and position of the center of mass based on the total force.

[0041] Step 4: State assembly and output. Assemble the new states obtained by integrating the five channels to form the complete 16-DOF vehicle state (VehicleState) for the next simulation moment, and output the state for use by external modules (such as controllers and data loggers).

[0042] Through the above methods, this invention provides a physically accurate and computationally efficient decoupled multi-degree-of-freedom dynamics simulation method. Five independent integration channels eliminate the non-physical cross-coupling generated by the inverse mass matrix in the ABA method, and in particular, solve the problem of inaccurate pitch / roll angle simulation in traditional simulations. Meanwhile, the appropriate use of the Pinocchio library and the channel-parallel integration strategy ensure the system's customizability and computational efficiency.

[0043] The symbols and abbreviations involved in this invention are explained below: I. Explanation of Symbols 1. Mass and Gravity m—Vehicle mass, unit: kg m_chassis — Sprout mass (the mass of the vehicle body and frame supported by the suspension), unit: kg mw—Unsprung mass (the mass below the suspension, including wheels, hubs, brake discs, etc.), unit: kg m_total—Total vehicle mass, m_total = m_chassis + 4×mw, unit: kg g—acceleration due to gravity, valued at 9.81 m / s² h_cg—Vehicle center of gravity height, i.e., the vertical distance from the center of gravity to the ground, unit: m. R—Effective radius of the wheel (rolling radius), unit: m Ftotal—The vector of the resultant force acting on the vehicle's center of mass, in N. v—Vehicle center of gravity velocity vector, unit: m / s 2. Moment of inertia and rotational moment of inertia Ixx—The moment of inertia (principal inertia) of the vehicle about its longitudinal axis (x-axis), also known as the roll moment of inertia, unit: kg·m² Iyy—The moment of inertia (principal moment of inertia) of the vehicle about its lateral axis (y-axis), also known as pitch moment of inertia, unit: kg·m² Izz—The moment of inertia (principal moment of inertia) of the vehicle about its vertical axis (z-axis), also known as the yaw moment of inertia, measured in kg·m². Ixz—The product of inertia (cross inertia) of the vehicle about the xz plane, reflecting the asymmetry of the mass distribution about the xz plane, unit: kg·m² I is a 3×3 inertia matrix, containing the principal inertia and the off-diagonal elements of the inertia product, in the form [Ixx, 0, Ixz; 0, Iyy, 0; Ixz, 0, Izz]. Iw—Moment of inertia of the wheel about its own axis of rotation, unit: kg·m² 3. Angle and angular velocity φ — Roll angle, the angle of rotation of the vehicle body about its longitudinal axis, unit: rad θ — Vehicle pitch angle, the angle of rotation of the vehicle body about its lateral axis, unit: rad ψ—Yaw angle, the angle of rotation of the vehicle body about its vertical axis, unit: rad ωx — vehicle body roll rate, i.e., dφ / dt, unit: rad / s ωy—the pitch velocity of the vehicle body, i.e., dθ / dt, unit: rad / s ωz — the yaw rate of the vehicle body, i.e., dψ / dt, unit: rad / s ω—wheel angular velocity, unit: rad / s δ—Wheel steering angle, unit: rad δcmd — Target steering angle (from controller command), unit: rad δFL, δFR — Steering angles of the left and right front wheels, in rad. 4. Force and Torque Mx—Roll moment, the moment acting on the lateral axis of the vehicle body, causing the vehicle body to roll. Unit: N·m My—Pitch moment, the torque acting on the longitudinal axis of the vehicle body, causing pitch motion (such as braking "nodding"), unit: N·m Mz—Yaw moment, the torque acting on the vertical axis of the vehicle body, causing the vehicle to yaw. Unit: N·m τx — Total external moment of the resultant roll axis, equal to the sum of the direct acting moment and the suspension coupling moment Mx, unit: N·m τy—The total external moment of the resultant pitch axis, equal to the sum of the direct acting moment and the suspension coupling moment My, in N·m. τz — The total external torque of the resultant yaw axis, equal to the sum of the directly acting torque and the external yaw torque Mz, in N·m. Fx — Tire-ground longitudinal force, tangential force along the direction of wheel travel, unit: N Fy—Tire-ground lateral force, a lateral force perpendicular to the direction of wheel travel, unit: N. Fz — Tire-ground normal force (vertical load), unit: N Fz_tire — Tire normal force, i.e., the vertical supporting force exerted on the tire by the ground, unit: N. Fcoupling—coupling force, an additional vertical force distributed to each suspension by the pitch and roll moments using the least squares method, measured in N. ΔFz_i — The coupling force component of the i-th suspension, in N Tdrive—Drive torque, the driving torque applied to the wheels by the controller, unit: N·m 5. Suspension parameters z—Current suspension compression (displacement), unit: m z0—Suspension static balance position (design compression), unit: m vsusp — Suspension compression speed, i.e., dz / dt, unit: m / s ks—Suspension spring stiffness coefficient, unit: N / m cs—Suspension damping coefficient, unit: N·s / m zFL, zFR, zRL, zRR — Front left, front right, rear left, and rear right suspension compression, unit: m 6. Tires and Skid s—Tire longitudinal slip ratio, dimensionless, reflects the degree of wheel slippage relative to the ground. v_ref—Low-speed regularization reference velocity, used to prevent the denominator of the slip ratio formula from approaching zero at low speeds; typical value is 0.5 m / s. vx — longitudinal velocity at the wheel center (in the direction of vehicle movement), unit: m / s ftire(s, Fz) — Tire force model function, takes slip ratio s and vertical load Fz as input, and outputs longitudinal force Fx. Supports linear models and Pacejka Magic Formula models. 7. Acceleration ax—Longitudinal acceleration of the vehicle (braking deceleration is negative), unit: m / s² ay—lateral acceleration of the vehicle, unit: m / s² 8. Coupling force calculation x_i — The longitudinal coordinate of the i-th suspension connection point relative to the center of mass, in meters. y_i — The lateral coordinate of the i-th suspension connection point relative to the center of mass, in meters. x —The mean of the four longitudinal coordinates of the suspension, x = Σx_i / 4 y —The mean of the four suspension lateral coordinates, y = Σy_i / 4 Sxx—Sum of squared deviations of the longitudinal coordinates of the suspension, Sxx = Σ(x_i - x )², used for least squares coupling force allocation Syy—The sum of squared deviations of the suspension's lateral coordinates, Syy = Σ(y_i - y )², used for least squares coupling force allocation 9. Quaternions q — attitude quaternion vector, q = [q0, q1, q2, q3] Used to represent three-dimensional rotations without singularity q0 — The scalar part (real part) of the quaternion, also known as cos(θ / 2). q1, q2, q3 — Quaternion vector part (imaginary part), related to the rotation axis and half-angle. Ω(ω)—from angular velocity ω = [ωx, ωy, ωz] The constructed oblique symmetric matrix is ​​used for the quaternion kinematic equation dq / dt = ½Ω(ω)q 10. Simulation and Integration dt—Simulation integration step size, with a value of 0.005 seconds (corresponding to a simulation frequency of 200Hz). dω / dt — Angular acceleration vector, i.e., the time derivative of angular velocity, unit: rad / s² dv / dt — Linear acceleration vector, i.e., the time derivative of velocity, unit: m / s² II. Explanation of Abbreviations WBC—Whole Body Control, is a multi-degree-of-freedom coordinated control method derived from robotics, which is applied in this invention to the coordinated control of multiple actuators in an automotive chassis. ABA—Articulated Body Algorithm, a recursive algorithm for solving the forward dynamics of multibody systems. DOF—Degree of Freedom, degree of freedom FK—Forward Kinematics HIL—Hardware-in-the-Loop—hardware-in-the-loop testing LTR—Load Transfer Ratio, a metric used to assess the risk of vehicle rollover / rollover. Pacejka Magic Formula—an empirical tire mechanics model proposed by HB Pacejka—fits the nonlinear relationship between tire force and slip ratio using a sine function and is widely used in vehicle dynamics simulation. Pinocchio is an open-source library for rigid body multibody dynamics calculations, providing efficient functions for calculating forward kinematics, Jacobian matrices, and mass matrices. Cramer's rule—a method for solving linear equations using determinants—is used in this invention to solve a 2×2 coupled angular acceleration subsystem containing the Ixz inertia product. The Least Squares Method, a mathematical optimization method, is used in this invention to optimally distribute external torques to each suspension connection point while satisfying pure couple constraints (zero net force). The semi-implicit Euler method, also known as the symmetric Euler method, is a numerical integration method that updates the velocity first and then the position. Compared with the explicit Euler method, it has better numerical stability and energy conservation properties.

[0044] This invention includes, but is not limited to, the above embodiments. Any equivalent substitutions or partial improvements made under the spirit and principles of this invention shall be considered within the scope of protection of this invention.

Claims

1. A decoupled multi-degree-of-freedom dynamic simulation method for an automobile chassis system, characterized in that, include: The vehicle's 16 degrees of freedom are divided into five independent integration channels: a 3-DOF attitude channel, a 4-DOF suspension channel, a 4-DOF wheel channel, a 2-DOF steering channel, and a 3-DOF center of mass channel. The attitude channel uses Euler's rigid body equations with inertia product Ixz and quaternion integration to independently integrate the vehicle's pitch, roll, and yaw angles. The suspension channel uses a spring-damping model to independently integrate the compression of the four suspension components. The wheel channel uses a tire slip ratio model to independently integrate the angular velocities of the four wheels. The suspension coupling forces generated by pitch and roll moments are calculated and simultaneously transmitted to both the suspension and attitude channels, forming a two-way physical feedback closed loop. Within the same simulation step, the five independent integration channels are integrated in parallel, with each channel using an integration method matched to its physical characteristics.

2. The method according to claim 1, characterized in that, During the attitude channel integration process, the synthesized total external moments τx, τy, and τz are calculated using the following formulas, directly incorporating the roll moment Mx and pitch moment My from the suspension system: τx = Mx + (Iyy - Izz) * ωy * ωz, τy = My + (Izz - Ixx) * ωz * ωx, τz = Mz + (Ixx - Iyy) * ωx * ωy, where ωx, ωy, and ωz are the roll, pitch, and yaw angular velocities, respectively; Ixx, Iyy, and Izz are the principal inertia; and Ixz is the inertia product. The above equations are expressed in matrix form as the inertia matrix I multiplied by the angular acceleration vector equals the vector on the right. The angular accelerations of each axis are obtained by solving the 2×2 coupled subsystem (dwx and dwz are cross-coupled through Ixz, and dwy is solved independently) using Cramer's rule.

3. The method according to claim 1, characterized in that, The dynamic equation of the spring-damped model of the suspension channel is: mw * dvsusp / dt = Fz_tire + (-ks * (z - z0) - cs * vsusp) + Fcoupling -mw * g, where mw is the unsprung mass (wheel-end mass), vsusp is the suspension compression velocity, ks is the suspension stiffness, cs is the suspension damping coefficient, z is the current compression, z0 is the static equilibrium position, and Fcoupling is the force from the coupling force channel.

4. The method according to claim 1, characterized in that, The dynamic equation of the wheel channel is: Iw * dω / dt = Tdrive - Fx * R, where Iw is the wheel moment of inertia, ω is the wheel angular velocity, Tdrive is the driving torque from the controller, Fx is the tire-ground longitudinal force, and R is the effective radius of the wheel; the tire slip ratio model is: slip ratio s = (ω * R - vx) / max(|vx|, v_ref), where v_ref is the low-speed regularized reference speed, longitudinal force Fx = ftire(s, Fz), where ftire is the tire force model, supporting linear and Pacejka Magic Formula models, where vx is the wheel center longitudinal velocity, and Fz is the tire vertical force.

5. The method according to claim 1, characterized in that, In step (e), under braking conditions, the pitching moment My is generated by the braking deceleration ax. The suspension coupling force distribution formula generated by this moment is as follows: The coupling force is calculated from the external moment using the least squares method to ensure zero net force constraint: For the pitching moment My, the coupling force of the i-th suspension is ΔFz_i = (My / Sxx) * (x_i - x ), where Sxx is the sum of squared deviations of the longitudinal coordinates of the suspension, x_i is the i-th longitudinal coordinate of the suspension, and x The mean is used; this method guarantees ΣΔFz_i=0 (pure couple) and ΣΔFz_i*x_i=My; for the roll moment Mx, the coupling force of the i-th suspension is ΔFz_i += (Mx / Syy) * (y_i - ), where Syy is the sum of squared deviations of the suspension lateral coordinates, and y_i is the i-th suspension lateral coordinate. The value is the average; the external torque is generated by acceleration: during braking, My = m*ax*h_cg, and during tilting, Mx = m*ay*h_cg, where h_cg is the height of the center of mass.

6. The method according to claim 1, characterized in that, The method also integrates the multibody dynamics library Pinocchio, but the Pinocchio library is not used to perform forward simulation of the articulated body algorithm ABA, but only to perform the following operations: forward kinematics calculation to determine the wheel contact point position and suspension joint angle; Jacobian matrix calculation to obtain the Jacobian matrix of the stability domain; And mass matrix lookup to provide inertia parameters.

7. The method according to claim 1, characterized in that, The simulation step size is 0.005 seconds, corresponding to a simulation frequency of 200Hz; the attitude channel, suspension channel, wheel channel, and center of gravity channel are integrated using a semi-implicit Euler method; the steering channel uses a position control method, directly writing the target steering angle into the simulation state vector without using differential equation integration.

8. The method according to any one of claims 1-7, characterized in that, It also includes a statistical comparison and verification step: on a scenario set covering multiple standardized test conditions, the whole-body control WBC mode, traditional control mode and passive mode are automatically run in parallel; standardized performance indicators are calculated, including load transfer ratio (LTR); an automatic runaway detection mechanism is introduced, which includes six runaway indicators, including sideslip and rollover indicators, and trigger thresholds and recovery thresholds are set for each indicator to form hysteresis protection.

Images