Dettagli della rivista
Formato
Rivista
eISSN
2444-8656
Prima pubblicazione
01 Jan 2016
Frequenza di pubblicazione
2 volte all'anno
Lingue
Inglese
Accesso libero

# Research on Dynamics of Flexible Multibody System with Deployable Antenna Based on Static Lagrangian Function

###### Accettato: 25 Apr 2022
Dettagli della rivista
Formato
Rivista
eISSN
2444-8656
Prima pubblicazione
01 Jan 2016
Frequenza di pubblicazione
2 volte all'anno
Lingue
Inglese
Modeling and analysis

Generally speaking, there are three dynamic modeling methods for flexible multi-body systems: The first is Newton-Euler method, which requires the momentum theorem to clarify the corresponding dynamic equations of all separated bodies, and then use constraints to remove the corresponding constraints. This method meets the requirement that there is a rigid body in the system. Secondly, it refers to the virtual displacement method, which should combine the virtual displacement principle to obtain the first Lagrangian equation, or combine the variational principle to construct the second Lagrangian equation, so as to obtain the dynamic model. Finally, it refers to other typical methods such as IGC, AIGC, NSD, etc.[1]

Definition 1. Euler parameters are used for analysis and Lagrange equation is selected to construct dynamic control equations of flexible multi-body system, namely differential equations (DAEs). The specific formula is as follows: ${M(q,t)q¨+GT(q,t)λ=F(q,q,t)g(q,t)=0$ \left\{{\matrix{{M\left({q,t} \right)\ddot q + {G^T}\left({q,t} \right)\lambda = F\left({q,q,t} \right)} \hfill \cr {g\left({q,t} \right) = 0} \hfill \cr}} \right.

Theorem 1. IIn the above formula, qRnq Represents the generalized coordinates of the system, M(g, t) ∈ Rnq×nq The generalized mass matrix representing positive definite symmetry of the system, g(q, t) : Rnq × RRnq Represents the left part of the constraint equation of generalized coordinates, λRn×q Stands for Lagrange times subvector, $G(q,t)=∂∂qg(q,t)∈Rnq×nq$ G\left({q,t} \right) = {\partial \over {\partial q}}g\left({q,t} \right) \in {R^{{n_q} \times {n_q}}}

Is the constraint matrix, F(q, q, t): Rnq × Rnq × RRnq Represents the generalized force of the system, which contains two contents: one is to refer to the generalized external force, the other is to refer to the velocity quadratic term. The dimension of this kind of equation exceeds the freedom degree of the system, and the actual calculation quantity is large. When the system is in a more exotic position, the actual value is more difficult to calculate.[2]

Proposition 2 In the study, Likins et al. proposed the use of mixed coordinates to represent the miniaturization of flexible multi-body systems. According to the conceptual analysis of relative and absolute coordinates of multi-rigid body system dynamics, it is mainly divided into two forms, one refers to the relative coordinates and elastic coordinates of floating coordinate system, and the other refers to the absolute coordinates and elastic coordinates of floating coordinate system. The specific formula is as follows: $q=(pTaT)T$ q = {\left({{p^T}{a^T}} \right)^T}

In the above formula, P represents relative or absolute coordinates in the floating coordinate system, and A represents deformation coordinates.

Modeling and numerical simulation are the core content of the study of dynamics of flexible multi-body system. It is very difficult to quickly construct the corresponding motion equations of the system, especially to form the numerical solution equations required by computer operation. Due to the large scale of solving the coefficients and variables of the dynamics equation of flexible multi-body system, the actual derivation process is complicated, so the symbolic calculus method is used to deal with the modeling of the system by using finite element method.[3]

Numerical analysis of DAEs

Lemma 3 After the establishment of the dynamic model of flexible multi-body system, it is necessary to study the dynamic change of the system with the combination of numerical simulation, and the study of numerical method is directly related to the dynamic control equations of the model. In the field of multi-body dynamics, this paper mainly studies the first kind of Lagrange equation, the specific form is as follows[4]: ${ddt(∂T∂q˙l)T−(∂T∂ql)T+Gq′T λ=Qii=1,…,nqg(q, t)=0$ \left\{{\matrix{{{d \over {dt}}{{\left({{{\partial T} \over {\partial {{\dot q}^l}}}} \right)}^T} - {{\left({{{\partial T} \over {\partial {q^l}}}} \right)}^T} + G_{{q^{'}}}^T\,\lambda = {Q^i}i = 1, \ldots,{n_q}} \hfill \cr {g\left({q,\,t} \right) = 0} \hfill \cr}} \right.

Conjecture 5. In the above formula, the unknown quantity contained in the first formula is not independent, and the actual quantity is consistent with the number of equations. The Lagrange multiplier λ={λ1, λ2 … λnλ} is also an unknown quantity that needs to be determined, and the actual number is consistent with the complete constraint equations. At this point, the constraint equations obtained are algebraic equations, which are also the core content of the numerical method studied in this paper.[5]

There are many methods to transform algebraic mixed equations into numerical models, such as augmented method, reduced method and mixed method. Based on the analysis of the algebraic equations obtained in the above research, it can be seen that: ${G(q,t)=(∂gi∂qj)nq×nqi=1,…,nλ, j=1, …,nqrank G(q, t)=nλ(nλ \left\{{\matrix{{G\left({q,t} \right) = {{\left({{{\partial {g_i}} \over {\partial {q_j}}}} \right)}_{{n_q} \times {n_q}}}i = 1, \ldots,{n_\lambda},\,j = 1,\, \ldots,{n_q}} \hfill \cr {rank\,G\left({q,\,t} \right) = {n_\lambda}\left({{n_\lambda} < {n_q}} \right)} \hfill \cr}} \right.

Example 6. Corresponding initial conditions are: $q0=q(0)q˙0=q˙(0)w q0=q(0)q˙0=q˙(0)w$ {q_0} = q\left(0 \right){\dot q_0} = \dot q{\left(0 \right)_w}\,{q_0} = q\left(0 \right){\dot q_0} = \dot q{\left(0 \right)_w}

Taking the derivative of t using the second expression of the algebraic equations, we can get: $G(q,t)⋅q˙+g1(g,t)=0G(q,t)⋅q¨+Φ(q,q˙,t)=0$ \matrix{{G\left({q,t} \right) \cdot \dot q + {g_1}\left({g,t} \right) = 0} \hfill \cr {G\left({q,t} \right) \cdot \ddot q + \Phi \left({q,\dot q,t} \right) = 0} \hfill \cr}

Among them: $η=−Φ(q,q˙,t)=−(G(q,t)⋅q˙)⋅q˙−2Gt(q,t)⋅q˙−gu(q,t)$ \eta = - \Phi \left({q,\dot q,t} \right) = - \left({G\left({q,t} \right) \cdot \dot q} \right) \cdot \dot q - 2{G_t}\left({q,t} \right) \cdot \dot q - {g_u}\left({q,t} \right)

The compatible conditions are as follows: $g(q0,0)=0G(q0,0)q˙0=−gt(q0,0)λ0=λ(q0,q˙0,0)$ \matrix{{g\left({{q_0},0} \right) = 0} \hfill \cr {G\left({{q_0},0} \right){{\dot q}_0} = - {g_t}\left({{q_0},0} \right)} \hfill \cr {{\lambda _0} = \lambda \left({{q_0},{{\dot q}_0},0} \right)} \hfill \cr}

According to an antenna in the flexible multi-body system dynamics model, the control equation of characteristic analysis, comparison on the relatively mature after the numerical simulation method, the constraint default direct correction method, the precise integration method, the nonlinear programming method acquired two, on the basis of continue to analyze the solving process, actual problems and computational cost, In the end, the direct default correction method of the augmented method is selected for solving analysis. The specific process is as follows:[6]

Note 7. The feature of the augmented method is to integrate the dynamics equation and acceleration constraint equation, so the generalized coordinates and Lagrange multiplier are regarded as concrete variables, and the equations as shown below can be obtained: $[M(q,t)GT(q,t)G(q,t)][q¨λ]=[F(q,q˙,t)η]$ \left[{\matrix{{M\left({q,t} \right){G^T}\left({q,t} \right)} \hfill \cr {G\left({q,t} \right)} \hfill \cr}} \right]\left[{\matrix{{\ddot q} \hfill \cr \lambda \hfill \cr}} \right] = \left[{\matrix{{F\left({q,\dot q,t} \right)} \hfill \cr \eta \hfill \cr}} \right]

The above formula conforms to $η=−Φ(q,q˙,t)=−(G(q,t)⋅q˙)⋅q˙−2Gi(q,t)⋅q˙−gu(q,t)$ \eta = - \Phi \left({q,\dot q,t} \right) = - \left({G\left({q,t} \right) \cdot \dot q} \right) \cdot \dot q - 2{G_i}\left({q,t} \right) \cdot \dot q - {g_u}\left({q,t} \right)

This condition. For nq+nλ unknown variables q and λ, the above equations belong to closed algebraic equations, and the actual solvability will be affected by the non-singularity of the coefficient matrix. In order to better analyze the application principle of direct default correction method, we can first define the method of constraint stability. In Baumgarte (Baumgarte) constrained default stability method, the following equations can be obtained because of the large error accumulation: ${g(q,t)=ε1g˙(q,t)=ε2g¨(q,t)=ε3$ \left\{{\matrix{{g\left({q,t} \right) = {\varepsilon _1}} \hfill \cr {\dot g\left({q,t} \right) = {\varepsilon _2}} \hfill \cr {\ddot g\left({q,t} \right) = {\varepsilon _3}} \hfill \cr}} \right.

Open Problem 8. The corresponding dynamic equation can be transformed into: $[M(q,t)GT(q,t)G(q,t)0][q¨λ]=[F(q,q˙,t)η*]$ \left[{\matrix{{M\left({q,t} \right){G^T}\left({q,t} \right)} \hfill \cr {G\left({q,t} \right)0} \hfill \cr}} \right]\left[{\matrix{{\ddot q} \hfill \cr \lambda \hfill \cr}} \right] = \left[{\matrix{{F\left({q,\dot q,t} \right)} \hfill \cr {{\eta ^*}} \hfill \cr}} \right]

The above formula conforms to $G(q,t)q¨=η*=−2α g˙(q,t)−β2g(q,t)+η$ G(q,t)\ddot q = {\eta ^*} = - 2\alpha \,\dot g(q,t) - {\beta ^2}g(q,t) + \eta This condition.

Generally speaking, α, β=5 ~ 50 should be set as the critical error resistance state under the condition of α=β, and the actual response stability is faster. The Lagrange multiplier can be obtained by combining the number of dynamic equations of flexible multi-body system and the actual derivative results: $λ=(GM−1GT)−1(GM−1F+Φ)$ \lambda = {\left({G{M^{- 1}}{G^T}} \right)^{- 1}}\left({G{M^{- 1}}F + \Phi} \right)

After sorting, it can be obtained: $Mq¨=F*$ M\ddot q = {F^*}

The above formula conforms to F* = FGT(GM−1GT)−1 (GM−1F + Φ) This condition.

On the one hand, coordinate default should be corrected. The numerical solution of the coordinate at time t=t (I) is assumed to be $q^(i)$ {\hat q^{\left(i \right)}} , Then affected by the error, the actual displacement constraint equation has a certain default phenomenon, and the equation as follows can be obtained: $g^(i)def__g(q^(i), t(i))≠0$ {\hat g^{\left(i \right)}}\underline{\underline {def}} g\left({{{\hat q}^{\left(i \right)}},\,{t^{\left(i \right)}}} \right) \ne 0

Example 7. When the breach of displacement constraint equation exceeds the specified range, numerical solution is needed $q^(i)$ {\hat q^{\left(i \right)}} Make corrections. Assuming that the correction term is, then it can be obtained: $q(i)=q^(i)+δq(i)$ {q^{\left(i \right)}} = {\hat q^{\left(i \right)}} + \delta {q^{\left(i \right)}}

Taylor expansion of the position constraint equation can obtain: $g(i)=g(q^(i), t(i))+g^q(i)⋅δq(i)=0$ {g^{\left(i \right)}} = g\left({{{\hat q}^{\left(i \right)}},\,{t^{\left(i \right)}}} \right) + \hat g_q^{\left(i \right)} \cdot \delta {q^{\left(i \right)}} = 0

Among them: $g^q(i)=gq(q^(i), t(i))=G(q^(i), t(i))$ \hat g_q^{\left(i \right)} = {g_q}\left({{{\hat q}^{\left(i \right)}},\,{t^{\left(i \right)}}} \right) = G\left({{{\hat q}^{\left(i \right)}},\,{t^{\left(i \right)}}} \right)

After sorting, it can be obtained: $G(q^(i), t(i))⋅δq(i)=−g^q(i)$ G\left({{{\hat q}^{\left(i \right)}},\,{t^{\left(i \right)}}} \right) \cdot \delta {q^{\left(i \right)}} = - \hat g_q^{\left(i \right)}

In general, because the number of generalized coordinates exceeds the number of constraint equations in the system, in other words G(q, t) ∈ Rnq × nq Is full row rank. Combined with the generalized force theory analysis of the matrix, because the constraint equation is not independent, so there is (GCT) −1, then the generalized inverse G+ (q, t) equation of the Jacobi matrix G (q, t) is: G+ (q, t) = GT(GGT−1

At this point, the minimal norm solution is obtained, as follows: $δq(i)=−G+(q^(i), t(i))g^(i)=−GT(q^(i), t(i))[G(q^(i), t(i))GT(q^(i), t(i))]−1⋅g(q^(i), t(i))$ \matrix{{\delta {q^{\left(i \right)}} = - {G^ +}\left({{{\hat q}^{\left(i \right)}},\,{t^{\left(i \right)}}} \right){{\hat g}^{\left(i \right)}} = - {G^T}\left({{{\hat q}^{\left(i \right)}},\,{t^{\left(i \right)}}} \right)} \hfill \cr {{{\left[{G\left({{{\hat q}^{\left(i \right)}},\,{t^{\left(i \right)}}} \right){G^T}\left({{{\hat q}^{\left(i \right)}},\,{t^{\left(i \right)}}} \right)} \right]}^{- 1}} \cdot g\left({{{\hat q}^{\left(i \right)}},\,{t^{\left(i \right)}}} \right)} \hfill \cr}

Minimal norm solution δq(i) As the least influential correction during the study of the dynamics equation, the generalized coordinates can be obtained which are closest to the real motion of the system, which is very important for numerical simulation research.

On the other hand, speed default is corrected. Combined with the above analysis, the hypothesis is consistent $g˙^(i)=G(q(i), t(i))⋅q˙^(i)+gt(q(i),t(i))=0$ {{\hat {\dot g}}^{\left(i \right)}} = G\left({{q^{\left(i \right)}},\,{t^{\left(i \right)}}} \right) \cdot {\hat {\dot q}^{\left(i \right)}} + {g_t}\left({{q^{\left(i \right)}},{t^{\left(i \right)}}} \right) = 0

For this condition, the modified item is δq(i), Then we can get: $q˙(i)=q˙^(i)+δq˙(i)$ {{\dot q}^{\left(i \right)}} = {{\hat {\dot q}}^{\left(i \right)}} + \delta {\dot q^{\left(i \right)}}

The constraint equation of the corresponding velocity is: $g˙(1)=G(q(i), t(i))⋅(q˙^(i)+δq˙(i))+gt(q(i),t(i))=0$ {{\dot g}^{\left(1 \right)}} = G\left({{q^{\left(i \right)}},\,{t^{\left(i \right)}}} \right) \cdot \left({{{\hat {\dot q}}^{\left(i \right)}} + \delta {{\dot q}^{\left(i \right)}}} \right) + {g_t}\left({{q^{\left(i \right)}},{t^{\left(i \right)}}} \right) = 0

By simplifying the conditions, we can get: $g˙(i)=g˙^(i)+G(q(i), t(i))δq˙(i)=0$ {\dot g^{\left(i \right)}} = {\hat {\dot g}^{\left(i \right)}} + G\left({{q^{\left(i \right)}},\,{t^{\left(i \right)}}} \right)\delta {\dot q^{\left(i \right)}} = 0

Thus, the equation of velocity correction value can be obtained as follows: $G(q(i),t(i))δq˙(i)=−g˙^(i)$ G\left({{q^{\left(i \right)}},{t^{\left(i \right)}}} \right)\delta {\dot q^{\left(i \right)}} = - {\hat {\dot g}^{\left(i \right)}}

The minimal norm solution of the above equation is studied by using the computational method of nonlinear programming method, and the specific results are shown as follows: $δq˙(i)=−G(i)T(G(i)G(i)T)−1g˙^(i)$ \delta {\dot q^{\left(i \right)}} = - {G^{\left(i \right)T}}{\left({{G^{\left(i \right)}}{G^{\left(i \right)T}}} \right)^{- 1}}{\hat {\dot g}^{\left(i \right)}}

In the study of the basic numerical methods of dynamic equations, the following points are specific: First, Treanor method. Assume that the unknown function value at tj is yij (I = 0,1... N-1), then calculate the unknown function value yi of tj+1=t+h, j+1, the formula is: $yi,j+1=yij+Δyi(i=0,1,…, n−1)$ {y_{i,j + 1}} = {y_{ij}} + \Delta {y_i}\left({i = 0,1, \ldots,\,n - 1} \right)

Where, when PI ≤0 is met, the following equation can be obtained: $Δyi=h6[qi(1)+2(qi(2)+qi(3))qi(4)]$ \Delta {y_i} = {h \over 6}\left[{q_i^{\left(1 \right)} + 2\left({q_i^{\left(2 \right)} + q_i^{\left(3 \right)}} \right)q_i^{\left(4 \right)}} \right]

In accord with pi&gt; When the condition is 0, we can get: $Δyi=h{qi(1)ri(2)+[−3(qi(1)+piwi(1))+2(qi(2)+piwi(2))+2(qi(3)+piwi(3))−(qi(4)+piwi(4))]4[(qi(1)+piwi(1))−(qi(2)+piwi(2))−(qi(3)+piwi(3))+(qi(4)+piwi(4))ri(4)]}$ \Delta {y_i} = h\left\{{\matrix{{q_i^{\left(1 \right)}r_i^{\left(2 \right)} + \left[{- 3\left({q_i^{\left(1 \right)} + {p_i}w_i^{\left(1 \right)}} \right) + 2\left({q_i^{\left(2 \right)} + {p_i}w_i^{\left(2 \right)}} \right) + 2\left({q_i^{\left(3 \right)} + {p_i}w_i^{\left(3 \right)}} \right) - \left({q_i^{\left(4 \right)} + {p_i}w_i^{\left(4 \right)}} \right)} \right]} \hfill \cr {4\left[{\left({q_i^{\left(1 \right)} + {p_i}w_i^{\left(1 \right)}} \right) - \left({q_i^{\left(2 \right)} + {p_i}w_i^{\left(2 \right)}} \right) - \left({q_i^{\left(3 \right)} + {p_i}w_i^{\left(3 \right)}} \right) + \left({q_i^{\left(4 \right)} + {p_i}w_i^{\left(4 \right)}} \right)r_i^{\left(4 \right)}} \right]} \hfill \cr}} \right\}

In the above formula: $pi=(qi(3)−qi(2))wi(3)−wi(2)(i=0,1,…,n−1)ri(1)=e−p1h, ri(2)=ri(1)−1−p1h, ri(3)=ri(2)−1−p1h, ri(4)=ri(3)−0.5−p1h(i=0,1,…, n−1) wi(1)=yijqi(1)=fi(ti, w0(1), w1(1),…,wn−1(1))wi(2)=wi(1)+h2qi(1)qi(2)=fi(tj+h2, w0(2), w1(2), …,wn−1(2))wi(3)=wi(1)+h2qi(2)qi(3)=fi(tj+h2, w0(3), w1(3), …,wn−1(3))wi(4)=wi(1)+qiqi(4)=fi(tj+h, w0(4), w1(4), …,wn−1(4))(i=0,1,…,n−1)$ \matrix{{{p_i} = {{\left({q_i^{\left(3 \right)} - q_i^{\left(2 \right)}} \right)} \over {w_i^{\left(3 \right)} - w_i^{\left(2 \right)}}}\left({i = 0,1, \ldots,n - 1} \right)} \hfill \cr {r_i^{\left(1 \right)} = {e^{- p{_1}h}},\,r_i^{\left(2 \right)} = {{r_i^{\left(1 \right)} - 1} \over {- {p_1}h}},\,r_i^{\left(3 \right)} = {{r_i^{\left(2 \right)} - 1} \over {- {p_1}h}},\,r_i^{\left(4 \right)} = {{r_i^{\left(3 \right)} - 0.5} \over {- {p_1}h}}\left({i = 0,1, \ldots,\,n - 1} \right)\,} \hfill \cr {w_i^{\left(1 \right)} = {y_{ij}}q_i^{\left(1 \right)} = {f_i}\left({{t_i},\,w_0^{\left(1 \right)},\,w_1^{\left(1 \right)}, \ldots,w_{n - 1}^{\left(1 \right)}} \right)} \hfill \cr {w_i^{\left(2 \right)} = w_i^{\left(1 \right)} + {h \over 2}q_i^{\left(1 \right)}q_i^{\left(2 \right)} = {f_i}\left({{t_j} + {h \over 2},\,w_0^{\left(2 \right)},\,w_1^{\left(2 \right)},\, \ldots,w_{n - 1}^{\left(2 \right)}} \right)} \hfill \cr {w_i^{\left(3 \right)} = w_i^{\left(1 \right)} + {h \over 2}q_i^{\left(2 \right)}q_i^{\left(3 \right)} = {f_i}\left({{t_j} + {h \over 2},\,w_0^{\left(3 \right)},\,w_1^{\left(3 \right)},\, \ldots,w_{n - 1}^{\left(3 \right)}} \right)} \hfill \cr {w_i^{\left(4 \right)} = w_i^{\left(1 \right)} + {q_i}q_i^{\left(4 \right)} = {f_i}\left({{t_j} + h,\,w_0^{\left(4 \right)},\,w_1^{\left(4 \right)},\, \ldots,w_{n - 1}^{\left(4 \right)}} \right)\left({i = 0,1, \ldots,n - 1} \right)} \hfill \cr}

Among them, qi should be defined in combination with the following formula: $qi={qi(3), pi<02(qi(3)−qi(1))ri(3)+(qi(1)−qi(2))ri(2)+qi(2),pi>0(i=0,1,…,n−1)$ {q_i} = \left\{{\matrix{{q_i^{\left(3 \right)},\,{p_i} < 0} \hfill \cr {2\left({q_i^{\left(3 \right)} - q_i^{\left(1 \right)}} \right)r_i^{\left(3 \right)} + \left({q_i^{\left(1 \right)} - q_i^{\left(2 \right)}} \right)r_i^{\left(2 \right)} + q_i^{\left(2 \right)},{p_i} > 0} \hfill \cr} \left({i = 0,1, \ldots,n - 1} \right)} \right.

Example 8. This method makes it easier to design programs. For general non-stiff equations, this method can be transformed into the fourth order Runge-Kutta method if the condition of PI ≤0 is met. And for the Stiff equation, it has the advantage of integral stability.[8]

Second, Gill's law. Kiel's formula for integrating from t to t+h is: $k0i=hfi(t,y0(0), y1(0),…yn−1(0))k1i=hfi(t+h2, y0(1), y1(1), …,yn−1(1))k2i=hfi(t+h2, y0(2), y1(2), …,yn−1(2))k3i=hfi(t+h, y0(3), y1(3), …,yn−1(3)) (i=0,1,…, n−1)$ \matrix{{{k_{0i}} = h{f_i}\left({t,y_0^{\left(0 \right)},\,y_1^{\left(0 \right)}, \ldots y_{n - 1}^{\left(0 \right)}} \right)} \hfill \cr {{k_{1i}} = h{f_i}\left({t + {h \over 2},\,y_0^{\left(1 \right)},\,y_1^{\left(1 \right)},\, \ldots,y_{n - 1}^{\left(1 \right)}} \right)} \hfill \cr {{k_{2i}} = h{f_i}\left({t + {h \over 2},\,y_0^{\left(2 \right)},\,y_1^{\left(2 \right)},\, \ldots,y_{n - 1}^{\left(2 \right)}} \right)} \hfill \cr {{k_{3i}} = h{f_i}\left({t + h,\,y_0^{\left(3 \right)},\,y_1^{\left(3 \right)},\, \ldots,y_{n - 1}^{\left(3 \right)}} \right)\,\left({i = 0,1, \ldots,\,n - 1} \right)} \hfill \cr}

Among them: $yi(1)=yi(0)+12(koi−2qi(0))yi(2)=yi(1)+(1−12)(k1i−qi(1))yi(3)=yi(2)+(1−12)(k2i−qi(2))yi(4)=yi(3)+16(k3i−2qi(3))(i=0,1,…, n−1)$ \matrix{{y_i^{\left(1 \right)} = y_i^{\left(0 \right)} + {1 \over 2}\left({{k_{oi}} - 2q_i^{\left(0 \right)}} \right)} \hfill \cr {y_i^{\left(2 \right)} = y_i^{\left(1 \right)} + \left({1 - \sqrt {{1 \over 2}}} \right)\left({{k_{1i}} - q_i^{\left(1 \right)}} \right)} \hfill \cr {y_i^{\left(3 \right)} = y_i^{\left(2 \right)} + \left({1 - \sqrt {{1 \over 2}}} \right)\left({{k_{2i}} - q_i^{\left(2 \right)}} \right)} \hfill \cr {y_i^{\left(4 \right)} = y_i^{\left(3 \right)} + {1 \over 6}\left({{k_{3i}} - 2q_i^{\left(3 \right)}} \right)\left({i = 0,1, \ldots,\,n - 1} \right)} \hfill \cr}

In the above formula: $qi(1)=qi(0)+3[12(k0i−2qi(0))]−12k0iqi(2)=qi(1)+3[(1−12)(k1i−qi(1))]−(1−12)k1iqi(3)=qi(2)+3[(1+12)(k2i−qi(2))]−(1+12)k2iqi(4)=qi(3)+3[16(k3i−2qi(3))]−12k3i(i=0,1,…, n−1)$ \matrix{{q_i^{\left(1 \right)} = q_i^{\left(0 \right)} + 3\left[{{1 \over 2}\left({{k_{0i}} - 2q_i^{\left(0 \right)}} \right)} \right] - {1 \over 2}{k_{0i}}} \hfill \cr {q_i^{\left(2 \right)} = q_i^{\left(1 \right)} + 3\left[{\left({1 - \sqrt {{1 \over 2}}} \right)\left({{k_{1i}} - q_i^{\left(1 \right)}} \right)} \right] - \left({1 - \sqrt {{1 \over 2}}} \right){k_{1i}}} \hfill \cr {q_i^{\left(3 \right)} = q_i^{\left(2 \right)} + 3\left[{\left({1 + \sqrt {{1 \over 2}}} \right)\left({{k_{2i}} - q_i^{\left(2 \right)}} \right)} \right] - \left({1 + \sqrt {{1 \over 2}}} \right){k_{2i}}} \hfill \cr {q_i^{\left(4 \right)} = q_i^{\left(3 \right)} + 3\left[{{1 \over 6}\left({{k_{3i}} - 2q_i^{\left(3 \right)}} \right)} \right] - {1 \over 2}{k_{3i}}\left({i = 0,1, \ldots,\,n - 1} \right)} \hfill \cr}

Where, yi (o) represents the position function value yi (t) at point T, yi (4) represents the unknown function value yi (t+ H) at point T + H after one step of integration, qi (O) represents that the initial value is 0, and qi (4) is regarded as the next qi (O) at each step of integration. In practice, this method can offset the rounding error accumulated in each step, effectively improve the accuracy of integration application, and has strong numerical stability for ill-posed dynamic equations. At the same time, the equation calculation has a lower stage error, compared with other methods can improve the integration speed, the actual calculation is also effectively controlled, the overall program has been improved. Third, the Gear method. The general expression formula is: $∑j=0kαjyn+j=hβkfn+kβk≠0, αk=1$ \sum\limits_{j = 0}^k {{\alpha _j}{y_{n + j}} = h{\beta _k}{f_{n + k}}{\beta _k} \ne 0,\,{\alpha _k} = 1}

Using Taylor expansion, the order is p=k, and the coefficient αj (j= 0,1 …) can be calculated under the condition that k= 1–6. K-1), βk, and the corresponding formula can be obtained: $k=1, yn+1=yn+hfn+1(The backward Euler method)k=2, yn+2=13(4yn−1−yn+2hfn+2)k=3, yn+3=111(18yn+2−9yn+1+2yn+6hfn+3)k=4, yn+4=125(48yn+3−36yn+2+16yn+1−3yn+12hfn+4)k=5, yn+5=1137(300yn+4−300yn+3+200yn+2−75yn+1+2yn+60hfn+5)k=6, yn+6=1147(360yn+5−450yn+4+400yn+3−255yn+2+72yn+1−10yn+60hfn+6)$ \matrix{{k = 1,\,{y_{n + 1}} = {y_n} + h{f_{n + 1}}(\rm The\,backward\,Euler\,method)} \hfill \cr {k = 2,\,{y_{n + 2}} = {1 \over 3}\left({4{y_{n - 1}} - {y_n} + 2h{f_{n + 2}}} \right)} \hfill \cr {k = 3,\,{y_{n + 3}} = {1 \over {11}}\left({18{y_{n + 2}} - 9{y_{n + 1}} + 2{y_n} + 6h{f_{n + 3}}} \right)} \hfill \cr {k = 4,\,{y_{n + 4}} = {1 \over {25}}\left({48{y_{n + 3}} - 36{y_{n + 2}} + 16{y_{n + 1}} - 3{y_n} + 12h{f_{n + 4}}} \right)} \hfill \cr {k = 5,\,{y_{n + 5}} = {1 \over {137}}\left({300{y_{n + 4}} - 300{y_{n + 3}} + 200{y_{n + 2}} - 75{y_{n + 1}} + 2{y_n} + 60h{f_{n + 5}}} \right)} \hfill \cr {k = 6,\,{y_{n + 6}} = {1 \over {147}}\left({360{y_{n + 5}} - 450{y_{n + 4}} + 400{y_{n + 3}} - 255{y_{n + 2}} + 72{y_{n + 1}} - 10{y_n} + 60h{f_{n + 6}}} \right)} \hfill \cr}

When K= 1,2, the method is a-stable, and when K= 3 ~ 6, the method is A (α) -stable and rigid stable. The results of A (α) -stable α Max and rigid stable parameters Dmin and θ Max are shown in the table 1 below:

Stability parameter analysis

K 1,2 3 4 5 6
αmax A - stable, 90° 86°54′ 73°14′ 51°50′ 18°47′
Dmin 0 0.1 0.7 2.4 6.1
Dmax 0.75 0.75 0.75 0.5

There are many countermeasures to improve this method, as shown below: $∑j=0kαjyn+j=h(βkfn+k+βk−1fn+k−1) αk≠0$ \sum\limits_{j = 0}^k {{\alpha _j}{y_{n + j}} = h\left({{\beta _k}{f_{n + k}} + {\beta _{k - 1}}{f_{n + k - 1}}} \right)\,{\alpha _k} \ne 0}

In the above formula, assuming βk is not equal to 0, and meets βk+β K-1 =1 condition, then we still need to calculate the order P =k in the formula, and βk=β is regarded as a free parameter, using Taylor expansion and k order compatible conditions, accurately calculate the corresponding expression formula.

Assuming β=20, the following formula can be obtained: $k=2, yn+2=141(80yn+1−39yn+40hfn+2−38hfn+1)k=3, yn+3=1182(417 yn+2−294 yn+1+59 yn+120hfn+3−114hfn+2)$ \matrix{{k = 2,\,{y_{n + 2}} = {1 \over {41}}\left({80{y_{n + 1}} - 39{y_n} + 40h{f_{n + 2}} - 38h{f_{n + 1}}} \right)} \hfill \cr {k = 3,\,{y_{n + 3}} = {1 \over {182}}\left({417\,{y_{n + 2}} - 294\,{y_{n + 1}} + 59\,{y_n} + 120h{f_{n + 3}} - 114h{f_{n + 2}}} \right)} \hfill \cr}

Compared with the Gear method of the same order, the actual stability region of these formulas is larger, and the corresponding results are shown in the following table:

Stability region comparison results

K 2 3 4 5 6 7
Improve gill method A - stable 89°55′ 85°32′ 73°2′ 51°23′ 18°32′
Jill method A - stable 86°54′ 73°13′ 51°50′ 18°47′ unstable

It can be seen that the disadvantage of this method is that the order of the equation keeps rising, the initial value is difficult to judge, and the multi-step time discrete has a comprehensive influence on the constraint equation, resulting in the drift of the solution. In other words, the rigidity should not be too high, and the Gear method will lose effect when the integral variable appears high-frequency oscillation, because the Gear method does not have A stability region.

Dynamic simulation program design and analysis of flexible multi-body system

The whole software system design includes the finite element analysis, structural design, the content such as geometry modeling, optimization design module, the deployment dynamics analysis module is the process in control, and to exchange data information of an operation, the corresponding numerical simulation module is master data information during the process control need to know. It can be seen that numerical simulation module is the basic requirement to ensure that the process control module can play an effective role. Flow chart 1 of the system platform constructed for deployable antenna in this paper is shown as follows:

The corresponding design system flow chart 2 is:

Flow chart 3 of the numerical calculation module is:

Key Contents

First, we have to deal with differential equations. In order to calculate the differential equation more quickly, it is necessary to deal with the equation obtained in the above research effectively. Because the equations are nQ-dimensional second order ordinary differential equations, so x1 = q, x2 = q This condition can be obtained: $x˙=[x2M−1F*]=f(x,t)$ \dot x = \left[{\matrix{{{x_2}} \hfill \cr {{M^{- 1}}{F^*}} \hfill \cr}} \right] = f\left({x,t} \right)

Among them, $x=[x1x2]$ x = \left[{\matrix{{{x_1}} \hfill \cr {{x_2}} \hfill \cr}} \right]

Assuming t=0, the initial formula of the system is: $x(0)=(q0Tq˙0T)T$ x\left(0 \right) = {\left({q_0^T\dot q_0^T} \right)^T}

Then the above equation of dynamics can be numerically integrated.

Secondly, we need to deal with numerical simulation methods. In this paper, a numerical integration method library is constructed, which contains a variety of different methods for effective selection during simulation. At present, in addition to several methods proposed in this paper, runge-Kutta method with fixed or variable step size and Gear method with variable step size and variable order are also involved. In addition, input more parameters, such as simulation start time, end time, etc.

Finally, the kinetic equation must be calculated and analyzed. On the one hand, for the multi-body system with closed loop, its dynamics equation is differential mixed equation, which needs to be contracted to obtain the pure differential form of the system dynamics equation, which is consistent with the dynamics equation of the tree multi-body system in form, and can be treated and solved in the same way. On the other hand, when the joints in the system have active drive, the simulation should focus on the justice mixing problem.[9,10]

Simulation Analysis

On the one hand, the study of the dynamics of rigid crank is mainly to verify the accuracy of the numerical simulation model. The Gill method of fixed step length is used to study the movement of rigid crank under the action of constant driving torque. The cross section of homogeneous crank is rectangular and shows a trend of horizontal rotation. At this point, the model parameters of the rigid crank involve the arm length L reached 0.75 meters, the cross-sectional area A was 60×15×10 square meters, and the mass M reached 1.4 kg. The driving torque is 62N•m, the crank is in horizontal position at the initial moment, and the initial conditions are as follows: $t=0, q(0)=0, q˙(0)=0$ t = 0,\,q\left(0 \right) = 0,\,\dot q\left(0 \right) = 0

Among them, the left figure is based on the results obtained from the theoretical research of traditional fixed-axis rotation mechanics, and the right figure is the results obtained from the theoretical simulation analysis of multi-body system dynamics in this paper. Combined with the model parameters, the actual angular acceleration of the crank is 236.19rad/s2. Therefore, it can be found that, on the one hand, the rigid model obtained according to the dynamic model is valid; on the other hand, the dynamic response obtained by the program is consistent with the traditional calculation, which further indicates the correctness of the design model studied in this paper.

On the other hand, the dynamics of two - link flexible manipulator is simulated. The corresponding model is shown in Figure 6 below:

Assuming that both links are in horizontal position at the initial moment and the actual deformation is 0, the initial conditions are as follows: $t=0, q(0)=0, q˙(0)=0$ t = 0,\,q\left(0 \right) = 0,\,\dot q\left(0 \right) = 0

Gill method with variable step size was used for analysis. The integral step size was H = LE-4, the simulation accuracy was EPS = le-L2, the simulation time was tend=1.2sec, and the calculation time was 2min56sec. The left and right graphs in FIG. 7 show the relation curves between angular velocity and time of two connecting rods at angles θ1 and θ2 respectively.

The Treanor method with fixed step length is used for analysis. Except that the calculation time is 47sec, other conditions are consistent. The relation between the actual angular velocity and time of the two connecting rods is shown in Figure 8 as follows:

Meanwhile, the above two methods are used to calculate the simulation results of h= LE-5 and H = LE-6 integral step sizes respectively. The actual calculation speed is shown in Table 3 below:

Speed comparison

The integral step Le-4 Le-5 Le-6
Variable step size Gill 2min56sec 48min6sec 2h43min53sec
Fixed step length Treanor 47sec 8min21sec 11min20sec

Combined with the analysis in Table 3 above, it can be seen that Gill method takes a long time to calculate a small integral step. However, Treanor method has no clear requirements on the integral step length, and the actual settlement time is very close. Compared with the corresponding dynamic change graphs, Gill's method has stronger stability for numerical solutions than Treanor's method, with higher actual accuracy. With the increase of calculation time, the angular velocity of the two links will continue to rise, which may cause dynamic hardening and increase the probability of data simulation failure.[11,12]

Conclusion

To sum up, the flexible multi-lift system dynamics studied in this paper integrates multiple disciplines. The actual development not only supplements and summarizes the original discipline, but also expands the thinking space of researchers and scholars, and promotes the construction of practical numerical simulation models to be more diversified. At present, the research work of flexible multibody system dynamics has obtained many research results from modeling scheme, symbolic calculus, tempering effect, numerical simulation and so on. Although related research is still a new subject, with the continuous exploration and data accumulation of domestic and foreign researchers, it is inevitable to start from this in the future development of deployable antenna application optimization. In addition, the future research direction of flexible multi-body dynamics will also change. For example, because the current high-precision integration algorithm cannot meet the requirements of implementation simulation, researchers should continue to discuss the appropriate simulation algorithm. In the calculation strategy, we can choose a new modeling method and computing system, such as recursive model, so as to make use of the parallel processing and vector processing capabilities of multi-processing computing, and conduct in-depth research on the simulation of relevant models.

#### Speed comparison

The integral step Le-4 Le-5 Le-6
Variable step size Gill 2min56sec 48min6sec 2h43min53sec
Fixed step length Treanor 47sec 8min21sec 11min20sec

#### Stability parameter analysis

K 1,2 3 4 5 6
αmax A - stable, 90° 86°54′ 73°14′ 51°50′ 18°47′
Dmin 0 0.1 0.7 2.4 6.1
Dmax 0.75 0.75 0.75 0.5

#### Stability region comparison results

K 2 3 4 5 6 7
Improve gill method A - stable 89°55′ 85°32′ 73°2′ 51°23′ 18°32′
Jill method A - stable 86°54′ 73°13′ 51°50′ 18°47′ unstable

Yiqun Zhang. Integrated Design of Structure and Control for Flexible Space Deployable Antenna [D]. Xidian University. ZhangYiqun Integrated Design of Structure and Control for Flexible Space Deployable Antenna [D] Xidian University Search in Google Scholar

Genyong Wu, Xingsuo He, P. Frank Pai. Dynamic Simulation of Large Deformation Flexible Multi-body System Based on Geometric Accurate Modeling Theory [C]// Chinese Congress on Mechanics. 2013. WuGenyong HeXingsuo PaiP. Frank Dynamic Simulation of Large Deformation Flexible Multi-body System Based on Geometric Accurate Modeling Theory [C] Chinese Congress on Mechanics 2013 Search in Google Scholar

Meizhi Li. Study on dynamics of flexible multi - bar system. Huazhong University of Science and Technology, 2004. MeizhiLi Study on dynamics of flexible multi - bar system Huazhong University of Science and Technology 2004 Search in Google Scholar

Dawei Zhu, Ping Zhu, Jiancheng Miao. Comparative Study on Dynamics of Rigid body and Flexible Body Based on Lagrange Method [J]. Computer Simulation, 2008, 025(001):314–319. ZhuDawei ZhuPing MiaoJiancheng Comparative Study on Dynamics of Rigid body and Flexible Body Based on Lagrange Method [J] Computer Simulation 2008 025 001 314 319 Search in Google Scholar

Hongxin Li. Analysis and Experimental Research on Deployable Planar Antenna Mechanism deployable Process [D]. Harbin Institute of Technology, 2015. LiHongxin Analysis and Experimental Research on Deployable Planar Antenna Mechanism deployable Process [D] Harbin Institute of Technology 2015 Search in Google Scholar

Pei Li, Cheng LIU, Qiang TIAN, et al. Dynamics of Deployable Antenna Structures with Large Flexible ring Trusses [C]// 2014 Conference on Deployable Spatial Structures. LiPei LIUCheng TIANQiang Dynamics of Deployable Antenna Structures with Large Flexible ring Trusses [C] 2014 Conference on Deployable Spatial Structures Search in Google Scholar

Shaoze YAN, Yuming GUAN, Jinwei FAN, et al. Dynamic Modeling of Flexible Multi-body System Based on Small Deformation [J]. Journal of Hebei University of Technology, 1999, 28(1):5. YANShaoze GUANYuming FANJinwei Dynamic Modeling of Flexible Multi-body System Based on Small Deformation [J] Journal of Hebei University of Technology 1999 28 1 5 Search in Google Scholar

Zhengfeng Bai, Yang Zhao, Hao Tian. Journal of Vibration and Shock, 2009, 28(6):4. BaiZhengfeng ZhaoYang TianHao Journal of Vibration and Shock 2009 28 6 4 Search in Google Scholar

Dazhi Cao, Hongfu Qiang, Gexue Ren. Parallel Computation of Flexible Multi-body System Dynamics Based on OpenMP and Pardiso [J]. Journal of Tsinghua University: Science & Technology, 2012, 52(11):7. CaoDazhi QiangHongfu RenGexue Parallel Computation of Flexible Multi-body System Dynamics Based on OpenMP and Pardiso [J] Journal of Tsinghua University: Science & Technology 2012 52 11 7 Search in Google Scholar

Selvi, M. Salai Mathi and Rajendran, L.. “Application of modified wavelet and homotopy perturbation methods to nonlinear oscillation problems” Applied Mathematics and Nonlinear Sciences, vol.4, no.2, 2019, pp.351–364. https://doi.org/10.2478/AMNS.2019.2.00030 SelviM. MathiSalai RajendranL. “Application of modified wavelet and homotopy perturbation methods to nonlinear oscillation problems” Applied Mathematics and Nonlinear Sciences 4 2 2019 351 364 https://doi.org/10.2478/AMNS.2019.2.00030 10.2478/AMNS.2019.2.00030 Search in Google Scholar

Hassan, Sk. Sarif, Reddy, Moole Parameswar and Rout, Ranjeet Kumar. “Dynamics of the Modified n-Degree Lorenz System” Applied Mathematics and Nonlinear Sciences, vol.4, no.2, 2019, pp.315–330. https://doi.org/10.2478/AMNS.2019.2.00028. HassanSk. SarifReddy ParameswarMoole RoutRanjeet Kumar “Dynamics of the Modified n-Degree Lorenz System” Applied Mathematics and Nonlinear Sciences 4 2 2019 315 330 https://doi.org/10.2478/AMNS.2019.2.00028. 10.2478/AMNS.2019.2.00028 Search in Google Scholar

Recuero A M, Escalona, José L. Dynamics of the coupled railway vehicle–flexible track system with irregularities using a multibody approach with moving modes[J]. Vehicle System Dynamics, 2014, 52(1):45–67. RecueroA M EscalonaJosé L. Dynamics of the coupled railway vehicle–flexible track system with irregularities using a multibody approach with moving modes[J] Vehicle System Dynamics 2014 52 1 45 67 10.1080/00423114.2013.857030 Search in Google Scholar

Articoli consigliati da Trend MD