A derivative-free algorithm for nonlinear equations and its applications in multibody dynamics

Abstract

Local floating coordinate system is used to represent the deployment motion of each rigid and flexible body of multibody system dynamics. Normal substructure modes are employed to describe the flexibility of a flexible body. Constraint equations establish the linkage between different bodies, part of them to specify positions and the others to specify orientations. System's governing equations are then derived using generalized coordinates by Lagrange methods. The resulting differential-algebraic equations are transformed to algebraic equations using backward differential formula corrector method, thus highly coupled nonlinear equations are obtained. However, Jacobian matrix of the nonlinear equations is hard to calculate, and then a quasi-Newton method based on Broyden–Fletcher–Goldfarb–Shanno update approach for the solution of the nonlinear equations is proposed. And a suitable line search approach is combined with the Broyden–Fletcher–Goldfarb–Shanno method to improve its efficiency. Some numerical results are reported to show efficiency of the proposed method. Afterwards, the Broyden–Fletcher–Goldfarb–Shanno method is integrated into multibody dynamics method. A rigid multibody case and a rigid-flex multibody case are further studied to show the efficiency of the proposed multibody solver.

Introduction

The multibody dynamic equations are described by differential-algebraic equations (DAEs) and the time integration method for DAE is a key issue in computational multibody dynamics. The accuracy, stability and efficiency of the method have direct influence on accuracy in multibody simulations, but unified solution method of DAE is not mature yet. Until now, many methods have been developed for DAE solver in multibody dynamics. ¹ One way to solve the DAE is to transfer the DAE to ordinary differential equations (ODE) by differentiating constraint equations. Then time-stepping methods of ODE, such as Newmark method, HHT-α method, generalized-α method, ${\theta }_1$ method, etc. can be used in structural dynamic simulations.^2–4 However, constraints may be violated during numerical integration of ODE. ⁵ And coefficients in stabilization of constraint equations are set empirically.

Another way to solve the DAE is to transfer the DAE to nonlinear equations. ⁶ Actually, in DAE terminology, dynamic system with holonomic constraint is of index 3.^7,8 Gear–Gupta–Leimkuhler formulation (index 2 formulation) ⁹ is often used to transfer the second order system to a first order system, thus positions and velocities of system are constrained simultaneously. The approach of Gear is useful because the multibody formalism is not modified during the solution process of the DAE essentially and it is suitable for stiff equations. For the application of general purpose, the DAE equations with Gear–Gupta–Leimkuhler formulation can be written in residual form ¹⁰

$$\begin{array}{c}\begin{array}{c}F(y,\stackrel{·}{y},t)=0\\ C(y,t)=0\end{array}\end{array}$$

(1)

where y and $\stackrel{·}{y}$ are generalized coordinates and derivative of generalized coordinates of the multibody system respectively. Then the previous equation (1) can be calculated by k-step backward differential formula (BDF) method and nonlinear equations will be obtained. Generally, nonlinear equations have the unified form

$$\begin{array}{c}A(x)=0,x\in R^n\end{array}$$

(2)

where x is an unknown variable and $A:R^n\to R^n$ is a continuously differentiable nonlinear operator. This implicit time integration method requires solving the nonlinear equations in each time step. Now nonlinear equations are solved by Newton ¹¹ and Newton-like methods iteratively in multibody dynamics. The backbone of multibody dynamics solver is specialized in numerical solution methods that are tailed to the structure of multibody dynamic equations. ¹² And system Jacobian matrix should be calculated in Newton method.

However, a drawback of these methods is that the numerical calculation of Jacobian matrices with respect to the subsystem variables is very time-consuming. Especially when the physics or engineering applications are related to multidiscipline, the nonlinear equations are usually numerically coupled and system Jacobian matrix may be too complex to calculate.^13,14 Thus quasi-Newton method with derivative-free will be an alternative one and the method is regarded as one of the most efficient methods of nonlinear equations. Since quasi-Newton method was first introduced by Broyden, ¹⁵ there have been significant progresses in the theoretical numerical analysis. Originally, Broyden and Broyden-like methods are local convergent,^16,17 thus initial value should be chosen carefully. A systematic and comprehensive review about local convergence analysis can be found by Dennis and More. ¹⁸ At that time, line search methods require calculation of derivatives, and it seems that the line search methods are inappropriate for global convergence analysis of quasi-Newton method. To overcome this shortcoming, an early method of globally derivative-free convergent nonlinear equation solver was presented by Griewank ¹⁹ with a derivate-free line search method. It is shown that if the Jacobian matrices are uniform nonsingular and Lipschitz continuously, then Griewank method converges to the unique solution. However, the linear search was not well defined all the time as pointed out in Li and Fukushima. ²⁰ Then another derivative-free line search method was proposed by Fukushima to resolve this difficulty. In order to obtain a larger step size, a derivative-free method was proposed by Youjun. ²¹ Derivative-free line search methods require backtracking, and a disadvantage of backtracking is that step size may be too small, which will lead to a round error. Then, some filter methods are also brought into line search methods.^22,23

The derivative-free method proposed by Youjun is simple and has a high computational efficiency. In the method, nonlinear equation (2) is solved by constructing an iterative algorithm, and linearized equations are used to represent the original equations

$$\begin{array}{c}L(x_k)=J_{k}x+b_k=0\end{array}$$

(3)

where k is the $\mathit{kth}$ iterative step, $J_k$ is usually the non-singular Jacobian matrix, and $b_k\in R^n$ . However, its calculating amount is large, the linearized equation (3) should be solved in each step and calculation cost of Broyden method is of order O( $n^3$ ). To reduce the calculation cost, inverse BFGS iteration method has been used to reduce calculation cost of order O( $n^3$ ) to order O( $n^2$ ). ²³ By combing the BFGS method with an exact line search or some special inexact line search, it was shown that a global convergence can be obtained. ²⁴ Yuan et al.^25,26 have proved the global convergence of a BFGS method for both convex functions and nonconvex functions under a modified weak Wolfe–Powell line search. A novel modified BFGS update method was developed for nonlinear equations ²⁷ and D. Chandra Sekhar applied the BFGS method to the aeroelastic problem of helicopter. ²⁸ Thus, inverse BFGS iteration method can combine the derivative-free line search method ²¹ to improve computational efficiency.

The purpose of this paper is to propose a new solver with derivative-free for multibody dynamics. An inverse BFGS method based on a derivative-free line search is developed and we bring it into DAE solver to simulate multibody dynamics. Thus, a new multibody dynamics solution method without Jacobian matrices calculation is established. In the method, multibody system dynamic equations are transferred to nonlinear equations by a predictor-corrector method, where predictor is using Taylor expansion method and corrector is using BDF method. Then the BFGS method is used to solve the nonlinear equations to overcome the problem of Jacobian matrix calculation.

A derivative-free BFGS method

Youjun ²¹ proposed a derivative-free line search based on the approximate norm descent condition below

$$\begin{array}{c}\Vert A(x_k+\alpha d_k)\Vert -\Vert A(x_k)\Vert \le -\sigma \Vert A(x_k+\alpha d_k)-A(x_k)\Vert +{\varepsilon }_k\end{array}$$

(4)

where $\sigma \in (0,1)$ and $\{{\varepsilon }_k\}$ is a positive sequence satisfying the following inequality for some constant $\varepsilon >0$

$$\begin{array}{c}\sum _{k=0}^{\infty }{\varepsilon }_k<\varepsilon <\infty \end{array}$$

(5)

In the previous equation, $\alpha$ is found by backtracking, and the right side of equation (4) goes to a positive constant, while the left side goes to zero. Thus, the iteration is well defined and a large step size of $\alpha$ can be obtained in the backtracking. The updating formula of estimation of Jacobian matrix proposed by Youjun is

$$\begin{array}{c}{J}_{k+1}=J_k+{\theta }_k\frac{(y_k-J_k{s}_k)s_k^T}{{\Vert s_k}^2\Vert }\end{array}$$

(6)

$$\begin{array}{c}{J}_{k}d+A(x_k)=0\end{array}$$

(7)

where $s_k=x_{k+1}-x_k$ , $y_k=A(x_{k+1})-A(x_k)$ , ${\theta }_k={\eta }^k$ .

However, it may spend a lot of computing cost on solving the linear equations of equation (7) especially when equation number becomes large. Instead, estimate the inverse of dynamic system Jacobian matrix will be an alternative way. ²³ In this paper, an inverse of Jacobian matrix updating formula (9) is combined with the line search method of formula (4). In the numerical experience section, the computing results are compared with the algorithm of paper Deng and Liu. ²¹

Now the algorithm of this paper is given as follows:

ALGORITHM 1.

Step 0: Initialization.

 Choose an initial point $x_0\in R^n$ and a nonsingular matrix $B_0\in R^{n\times n}$ , Calculate $A(x_0)$ . Given a positive sequence $\{{\varepsilon }_k\}$ satisfying equation (4). Choose constants $\sigma$ , $\eta \text{and}\rho \in (0,1)$ . Select the final accuracy tolerance $\delta >0$ . Let $k:=0$ .

Step 1: Test for optimality.

 If $\Vert A(x_k)\Vert <\delta$ then stop, else calculate $d_k$ by the equation

$$d=-B_{k}A(x_k)$$

(8)

for $d_k$ . Let ${\alpha }_k:=1$ .

Step 2: Backtracking.

 If $\Vert A(x_k+\alpha d_k)\Vert -\Vert A(x_k)\Vert =-\sigma \Vert A(x_k+\alpha d_k)-A(x_k)\Vert +{\varepsilon }_k$ , then let $x_{k+1}:=x_k+{\alpha }_k{d}_k$ , else let ${\alpha }_k:=\rho {\alpha }_k$ and back to Step 2.

Step 3: Update.

Step 3.1. Let $s_k=x_{k+1}-x_k$ , $y_k=A(x_{k+1})-A(x_k)$ , ${\theta }_k:=1$ .

Step 3.2

$$\begin{array}{c}{B}_{k+1}=B_k+{\theta }_k(-\frac{B_k{s}_k{s}_k^T{B}_k}{s_k^T{B}_k{s}_k}+\frac{y_k{y}_k^T}{y_k^T{s}_k})\end{array}$$

(9)

Step 3.3 If det ( $B_{k+1}$ ) $\ne 0$ then set $k:=k+1$ and go to Step 1, else set ${\theta }_k:=\eta {\theta }_k$ and go to Step 3.2.

Multibody dynamics

Kinemics analysis

Rigid body motions are represented by the motions of floating frame of reference, and the elastic motions are represented by modal shapes. The total generalized coordinates of the ith flexi body are $q_i=[x_i,y_i,z_i,{\varphi }_i,{\theta }_i,{\psi }_i,{\varsigma }_i^1,\dots ,{\varsigma }_i^k{]}^T$ , where $[x_i,y_i,z_i{]}^T$ are locations of a body, $[{\varphi }_i,{\theta }_i,{\psi }_i{]}^T$ are Euler angles, ${\varsigma }_i=[{\varsigma }_i^1,\dots ,{\varsigma }_i^k{]}^T$ are generalized coordinates of the flexible modal shapes. ²⁹ When a body moves with elastic deformations, shown as Figure 1, point $P\text{'}s$ location of body i can be written as

$$r_i=r_{i0}+A_i(u_{i0}+u_{\mathit{if}})$$

(10)

where $r_{i0}$ is the location of ith body center under the global Cartesian coordinate, $A_i$ is the Euler rotation matrix of the floating frame, $u_{i0}$ is location under the floating frame. $u_{\mathit{if}}=\sum _{p=1,k}{\Psi }_i^p{\varsigma }_i^p$ represents the elastic deformations, while ${\Psi }_i=[{\Psi }_i^1,\dots ,{\Psi }_i^k]$ is modal shape of the ith flexible body. Velocity vector of the point P can be obtained by differentiating equation (10)

$$\begin{array}{c}{\stackrel{·}{r}}_i={\stackrel{·}{r}}_{i0}+{\stackrel{·}{A}}_i(u_{i0}+u_{\mathit{if}})+A_i{\Psi }_i{\varsigma }_i\end{array}$$

(11)

Figure 1.

Flex body in the reference and deformed configurations.

A body may suffer loads during the process of movement, all loads are acted in form of node forces under the global coordinate. The loads are divided into two categories, the force $F_i=[f_{\mathit{ix}},f_{\mathit{iy}},f_{\mathit{iz}}{]}^T$ , and the Toque $T_i=[T_{\mathit{ix}},T_{\mathit{iy}},T_{\mathit{iz}}{]}^T$ . All loads are transferred to the generalized force thorough the principle of virtual work

$$\begin{array}{c}{Q}_{\mathit{ie}}=[F_i,G_i^T(T_i+h_i\times F_i),{\Psi }_i^T{\text{A}}_i^T{F}_i+{\Psi }_i^{*T}{\text{A}}_i^T{T}_i{]}^T\end{array}$$

(12)

where $G_i$ is a transfer matrix, $h_i$ is location vector, ${\Psi }_i^{*}$ is rotational mode at acting point. ²⁹

Constraints

Different parts of the multibody system are connected by kinds of joints. Joints are modeled by constraint equations, which are composed of location constraints and vector constraints. ³⁰ Assume a joint connects ith body and jth body at point m and point n, respectively. Location constraints are constructed at point m and point n, while z axis of one joint parallels to the other as shown in Figure 2. The location constraint equations are

$$\begin{array}{c}{r}_{\mathit{im}}-r_{\mathit{jn}}=0\end{array}$$

(13)

where $r_{\mathit{im}}$ , $r_{\mathit{jn}}$ are locations of point m and point n respectively under the global coordinate. The vector constraint equations are constructed to meet the relative angular requirements of the two bodies as shown in equation (14).

$$v_m^x{v}_n^z=0$$

(14)

Figure 2.

Joint connection.

Figure 3.

Multibody method based on BFGS method. BFGS: Broyden–Fletcher–Goldfarb–Shanno. ²⁷

A fixed joint and a rotational motion joint are modeled in this paper. ³⁰ The equations of fixed joint are shown in equation (15), while the equations of rotational motion joint are shown in equation (16).

$$\begin{array}{c}{C}_{\mathit{fixed}}=[\begin{array}{c}{r}_{\mathit{im}}-r_{\mathit{jn}}=0\\ v_m^x{v}_n^z=0\\ v_m^y{v}_n^z=0\\ v_m^y{v}_n^x=0\end{array}]\end{array}$$

(15)

$$\begin{array}{c}{C}_{\mathit{driving}}=[\begin{array}{c}{r}_{\mathit{im}}-r_{\mathit{jn}}=0\\ v_m^x{v}_n^z=0\\ v_m^y{v}_n^z=0\\ v_m^y{v}_n^x=\text{sin}(\alpha t)\end{array}]\end{array}$$

(16)

where α is the angle of rotational motion.

Governing equations

Assume a multi body system is composed of u rigid bodies, v flexible bodies and k constraints, the system freedom is n. The generalized coordinates are $q=[q_1,\dots ,q_n{]}^T$ , while constraint equations are $C_i(q,t)=0,i=1,\dots ,s$ . Because the generalized coordinates are not independent, thus Lagrange multiplier method is used to transfer the constraint equations to the kinematic equations as follows ³⁰

$$Q_I+Q_e-\mathit{Kq}-C_q^T\lambda =0$$

(17)

where $Q_I=-\int {\mathit{no}}_v\rho (\frac{\partial r}{\partial q}{)}^T\stackrel{\mathrm{··}}{r}\mathit{dv}$ is the generalized inertial term, K is the stiffness matrix, $Q_e$ is the generalized force. λ is the Lagrange multiplier, $C_q$ is the partial derivative of constraint equations to the generalized coordinates q, ρ is the density, $\mathit{dv}$ is the finite volume of a body. The DAEs of the multi body system are obtained by combining equations (15) to (17) ³¹

$$\begin{array}{c}\{\begin{array}{c}M\stackrel{\mathrm{··}}{q}+\mathit{Kq}+C_q^T\lambda =Q_e+Q_v\\ C(q,t)=0\end{array}\end{array}$$

(18)

where M and $Q_v$ are the mass matrix and velocity part of generalized inertial force of the system.

Numerical integration scheme based on Broyden method

Let $u=\stackrel{·}{q}$ , the differential equations of DAEs (18) are rewritten

$$\begin{array}{c}\{\begin{array}{c}\stackrel{·}{q}=u\\ M\stackrel{·}{u}+\mathit{Kq}+C_q^T\lambda =Q_e+Q_v\end{array}\end{array}$$

(19)

Let $\phi =[q,u,\lambda {]}^T$ and then equation (19) can be rewritten in a uniform formula ³²

$$\begin{array}{c}\{\begin{array}{c}F(\phi ,\stackrel{·}{\phi },t)=0\\ C(\phi ,t)=0\end{array}\end{array}$$

(20)

When solving equation (20) at time step $t+\Delta t$ , the initial value ${\phi }_{t+\Delta t}$ is predicted by 2th order Taylor predictor

$$\begin{array}{c}{\phi }_{t+\Delta t}={\phi }_t+{\stackrel{·}{\phi }}_t\Delta t+\frac{1}{2}{\stackrel{\mathrm{··}}{\phi }}_t\Delta t^2\end{array}$$

(21)

Then equation (20) are corrected by 4 steps BDF formula

$$\begin{array}{c}\frac{12}{25}\mathit{hf}({\phi }_{n+4},t_{n+4})={\phi }_{n+4}-\frac{48}{25}{\phi }_{n+3}+\frac{36}{25}{\phi }_{n+2}\\ -\frac{16}{25}{\phi }_{n+1}+\frac{3}{25}{\phi }_n\end{array}$$

(22)

Thus, nonlinear equations are obtained. And inverse BFGS method with Jacobean-free of ALGORITHM 1 is used to solve the equations.

When at initial time step, initial velocities ${\stackrel{·}{q}}_0$ of system are solved using minimal optimization

$$\begin{array}{c}\{\begin{array}{c}\mathit{Minimize}:L=\frac{1}{2}{\stackrel{·}{q}}^T\stackrel{·}{q}\\ \mathit{Subjectto}:C_q\stackrel{·}{q}+C_t=0\end{array}\end{array}$$

(23)

Thus, linear equations are obtained

$$\begin{array}{c}[\begin{array}{cc}I& \hfill C_q^T\hfill \\ C_q& \hfill 0\hfill \end{array}][\begin{array}{c}\stackrel{·}{q}\\ \lambda \end{array}]=[\begin{array}{c}0\\ C_t\end{array}\end{array}]$$

(24)

Initial velocities $\stackrel{·}{q}$ can be solved and velocity part $Q_v$ can be obtained. The initial ${\stackrel{\mathrm{··}}{q}}_0$ and $\lambda$ can be calculated by solving equation (18) based on the inverse BFGS method.

The algorithm is given as follows and the calculation flow diagram is shown in Figure 3:

ALGORITHM 2:

Step 0: Initialize parameters.

The multibody model is built and the initial $q_0$ is set. Initial velocity ${\stackrel{·}{q}}_0$ is obtained by solving the linear equations (24). Set default value of BFGS algorithm 1: $\sigma$  = 0.5, ${\gamma }_{\phi }$  = 0.5, $\rho =0.2$ , ${\varepsilon }_k=0.5/{\text{k}}^2$ , $\eta$  = 0.8, $\delta ={10}^{-10}$ . Set $B_0=I$ and $\lambda =0$ , where I is the identity matrix.

Step 1: Multibody and constraint modeling.

Transfer the $q_0$ and velocity ${\stackrel{·}{q}}_0$ of system to each body and connector. Update position, mass matrix and stiffness matrix of each body. Update location and vector constraint equations of each connector.

Step 2: DAEs, Taylor predictor and BDF corrector.

Assemble ODE of each body and algebraic equations of each connector and DAE of system is obtained. The initial values at each time step are predicted by the 2th order Taylor predictor of equation (21). DAE are corrected using BDF formula and nonlinear equations are obtained.

Step 3: BFGS system.

Transfer freedom of multibody system to nonlinear equations and independent variables of nonlinear equations are composed of system freedom. Thus, ALGORITHM 1 is used to solve the nonlinear equations.

Step 4: Solve the nonlinear equations until a convergence is obtained. Transfer solution of nonlinear equations to freedom of multibody system. If solution is at the end then output results, else go to step 1 and calculate for the next time step.

Numerical experiences and discussion

In this section, to verify the solver of nonlinear equations, some nonlinear equations case is studied. Then the method of nonlinear equations is integrated to multibody dynamics, a multi-rigid body case and a rigid-flex body case are studied.

Numerical experience of nonlinear equations

We wrote a FORTRAN implementation of the Algorithm 1 and tested the following set of test problems from paper. ²¹ Calculations results of Algorithm 1 were compared with those of algorithm in Deng and Liu. ²¹ It was summarized in Tables 1 and 2, in which the following notation is used.

Table 1.

Numerical comparison.

No.(number of variables)/ Algorithm		I	II
1(1)	BKT	0	2
	NI	28	6
	NT
	NORM	6.18648974E-010	3.69436467E-0012
2(3)	BKT	22	24
	NI	19	17
	NORM	4.77143839E-10	2.38951784E-0012
3(2)	BKT	2	1
	NI	27	31
	NORM	2.05218333E-10	3.24256509E-0011
4(3)	BKT	0	0
	NI	3	2
	NORM	8.94069263E-11	8.94069726E-0011
5(3)	BKT	6	7
	NI	14	12
	NORM	2.77388197E-11	2.28055025E-0014
6(6)	BKT	889	1651
	NI	211	333
	NORM	2.81484198E-10	9.81004324E-0012
7(3)	BKT	17	12
	NI	19	18
	NORM	1.06859940E-10	1.20985694E-13
8(4)	BKT	40	29
	NI	27	23
	NORM	4.13369440E-10	3.52400505E-0011
9(2)	BKT	13	11
	NI	19	14
	NORM	5.83911549E-12	4.82542583E-11
10(3)	BKT	4	4
	NI	14	12
	NORM	1.62014066E-11	1.07978781E-0012
11(2)	BKT	12	3
	NI	17	11
	NORM	1.50045571E-10	2.21079707E-0012
12(2)	BKT	915	199
	NI	327	114
	NORM	3.77680670E-11	4.07477480E-0011
13(2)	BKT	9	9
	NI	16	17
	NORM	3.06880163E-10	2.44920164E-0013
14(2)	BKT	2	2
	NI	10	9
	NORM	1.25580131E-10	2.59750844E-0011
15(3)	BKT	22	11
	NI	22	14
	NORM	2.38350824E-10	2.08273690E-0012
16(2)	BKT	8	9
	NI	15	12
	NORM	1.54408415E-11	8.79997619E-0012
17(3)	BKT	15	9
	NI	15	12
	NORM	4.21043151E-10	2.89976083E-0011

Table 2.

The solution.

No./ Algorithm		I	II
1(1)	X(1)	−1.56994424461096	−0.52359877559706
2(3)	X(1)	0.99999999999626	0.99999999999996
	X(2)	1.00000000004168	1.00000000000024
	X(3)	1.00000000001225	1.00000000000002
3(2)	X(1)	−0.64332284412465	−0.64367630882303
	X(2)	−2.21367347966790	−2.21491060590535
4(3)	X(1)	0.000000000000000	0.00000000000000
	X(2)	0.100000005215406	0.10000000521540
	X(3)	1.00000000000000	1.00000000000000
5(3)	X(1)	0.908926369284422	0.90892636926707
	X(2)	1.08560001168726	1.08560001167168
	X(3)	0.682147261469357	0.68214726146584
6(6)	X(1)	4.98959204978069	−4.39918312687848
	X(2)	−13.0455080951923	−9.82147462555338
	X(3)	10.2101761241605	16.49336143134665
	X(4)	−4.39918312687750	4.98959204978883
	X(5)	9.02808129598120	9.90391544159937
	X(6)	−8.63937979736478	−11.78097245096115
7(3)	X(1)	1.99607232807405	1.67654644242118
	X(2)	1.10965799691528	1.43125749183231
	X(3)	2.37032108681996	2.72624304638332
8(4)	X(1)	0.999999999968179	1.00000000000430
	X(2)	8.489888375268451E−012	1.75401157458567E−12
	X(3)	−1.00000000004407	−0.99999999999357
	X(4)	0.500000000035616	0.49999999999627
9(2)	X(1)	−0.707106830988472	−8.77940724716718E−12
	X(2)	1.50000007042472	0.99999999995471
10(3)	X(1)	−0.767760562423625	−0.76776056243563
	X(2)	7.533067640423938E−002	7.53306764041538
	X(3)	−1.16215118587028	−1.16215118585959
11(2)	X(1)	4.02858655786816	0.13927499920091
	X(2)	1.96721290655287	4.28294321065657
12(2)	X(1)	1.00000000000741	1.00000000000822
	X(2)	1.00000000001476	1.00000000001640
13(2)	X(1)	1.33333334611025	1.33333331302455
	X(2)	0.577350279992796	0.57735025160186
14(2)	X(1)	0.998606944112507	0.99860694409400
	X(2)	−0.105530492267675	−0.10553049229841
15(3)	X(1)	0.258405740156031	0.25840574019806
	X(2)	−0.251681148001244	−0.25168114803984
	X(3)	−1.86988307312539	−1.86988307315382
16(2)	X(1)	1.06734608580823	1.06734608580459
	X(2)	0.139227666899697	0.13922766688039
17(3)	X(1)	0.499999999993126	0.49999999999946
	X(2)	−5.322180602400565E−009	−7.08341966140051E−0009
	X(3)	−0.523598780931556	−0.52359881418775

I: Algorithm 1 of this paper, II: Algorithm 1 in Deng and Liu. ²¹

BKG: Backtracking steps, NI: Number of iterations.

NORM: Norm of $A(x_k)$ of the last iterative step.

x(i): ith component if independent variable x.

Equations are listed as follows

$$\text{cos}(3x_1)=0$$

(25)

$$\begin{array}{c}\begin{array}{c}2x_1{x}_2-x_2+5x_3-6=0\\ x_2^2+3x_1-10x_2{x}_3+6=0\\ x_1+2x_2-5x_3+2=0\end{array}\end{array}$$

(26)

$$\begin{array}{c}\begin{array}{c}0.7\text{sin}(x_1)-0.2\text{cos}(x_2)+0.3=0\\ 0.7\text{cos}(x_1)-0.2\text{sin}(x_2)-0.4=0\end{array}\end{array}$$

(27)

$$\begin{array}{c}\begin{array}{c}{x}_1+\text{cos}(x_1{x}_2{x}_3)-1=0\\ {(1-x_1)}^{0.25}+x_2+0.05x_3^2-0.15x_3-1=0\\ -x_1^2-0.1x_2^2+0.01x_2+x_3-1=0\end{array}\end{array}$$

(28)

$$\begin{array}{c}\begin{array}{c}12x_1-x_2^2-4x_3-7=0\\ x_1^2+10x_2-x_3-11=0\\ x_2^2+10x_3-8=0\end{array}\end{array}$$

(29)

$$\begin{array}{c}\begin{array}{c}\begin{array}{c}{x}_1^2\text{cos}(x_2)\text{cos}(x_3)+x_4^2\text{cos}(x_5)\text{cos}(x_6)+3.0=0\\ x_1^2\text{cos}(x_2)\text{sin}(x_3)+x_4^2\text{cos}(x_5)\text{sin}(x_6)+3.0=0\\ -x_1^2\text{sin}(x_2)-x_4^2\text{sin}(x_5)-4.0=0\end{array}\\ \begin{array}{c}{x}_1\text{cos}(x_2)\text{cos}(x_3)+x_4\text{cos}(x_5)\text{cos}(x_6)+6.0=0\\ x_1\text{cos}(x_2)\text{sin}(x_3)+x_4\text{cos}(x_5)\text{sin}(x_6)+6.0=0\\ -x_1\text{sin}(x_2)-x_4\text{sin}(x_5)-4.0=0\end{array}\end{array}\end{array}$$

(30)

$$\begin{array}{c}\begin{array}{c}{x}_1^2{x}_2^2+x_1^2+x_2-10=0\\ x_2^2{x}_3^2+x_2^2+x_3-20=0\\ x_1^2{x}_3^2+x_3^2+x_1-30=0\end{array}\end{array}$$

(31)

$$\begin{array}{c}\begin{array}{c}\begin{array}{c}2x_1+6x_2-6x_3{x}_4-5=0\\ {(x_1+x_2)}^2-4x_3-10x_4=0\end{array}\\ \begin{array}{c}2.5x_1+2x_2-2x_3-3x_4-4=0\\ x_2{x}_2+x_3+x_4+0.5=0\end{array}\end{array}\end{array}$$

(32)

$$\begin{array}{c}\begin{array}{c}{x}_1^2-x_2+1=0\\ x_1-\text{cos}(\pi x_2/2)=0\end{array}\end{array}$$

(9)

$$\begin{array}{c}\begin{array}{c}{x}_1^3+e^{x_1}+2x_2+x_3+1=0\\ -x_1+x_2+x_2^3+2e^{x_2}-3=0\\ -2x_2+x_3+e^{x_3}+1=0\end{array}\end{array}$$

(33)

$$\begin{array}{c}\begin{array}{c}{e}^{x_1^2+7x_2-30}-1=0\\ x_1^2+7x_2-30=0\end{array}\end{array}$$

(11)

$$\begin{array}{c}\begin{array}{c}200(x_2-x_1^2)=0\\ -2.0(x_1(200(x_2-x_1^2)-1)+1)=0\end{array}\end{array}$$

(34)

$$\begin{array}{c}\begin{array}{c}\frac{1}{2}{x}_1[\text{sin}(\frac{1}{2}\pi x_1)]-x_2=0\\ x_2^2-x_1+1=0\end{array}\end{array}$$

(35)

$$\begin{array}{c}\begin{array}{c}4x_1^2+x_2^2-4=0\\ x_1+x_2-\text{sin}(x_1-x_2)=0\end{array}\end{array}$$

(36)

$$\begin{array}{c}\begin{array}{c}{x}_1^2+x_2^2-x_3-2=0\\ x_1+5x_2+1=0\\ x_1{x}_3-2x_1+1=0\end{array}\end{array}$$

(37)

$$\begin{array}{c}\begin{array}{c}{x}_1^2-x_2-1=0\\ {(x_1-2)}^2+{(x_2-0.5)}^2-1=0\end{array}\end{array}$$

(38)

$$\begin{array}{c}\begin{array}{c}3x_1-\text{cos}(x_2{x}_3)-\frac{1}{2}=0\\ x_1^2-81{(x_2+0.1)}^2+\text{sin}(x_3)+1.06=0\\ e^{-x_1{x}_2}-20x_3+\frac{1}{3}(10\pi -3)=0\end{array}\end{array}$$

(39)

The parameters of these algorithms are the same as in paper Deng and Liu ²¹ as follows, $\sigma$  = 0.5, ${\gamma }_{\phi }$  = 0.5, $\rho$  = 0.2, ${\varepsilon }_k=0.5/{\text{k}}^2$ , $\eta$  = 0.8, $\delta ={10}^{-10}$ , $x_0=(0,0,\dots ,0{)}^T$ . $B_0=I$ .

From Table 1, we see that some of the equations do not satisfy uniform nonsingular assumptions; however, algorithm of this paper still generates converge sequences. Norm of $A(x_k)$ of the last iterative step of this paper is smaller than Deng and Liu ²¹ in each nonlinear problem. Table 2 shows the value of the final iterative point.

To check the total calculating steps, total steps TS of each nonlinear equation are calculated by adding backtracking steps and number of iterations. Then BKT, NI, TS of this paper are compared with those of Deng and Liu, ²¹ the results are summarized in Table 3. It is shown that there are 8 nonlinear problems whose backtracking steps of this method are smaller than Deng and Liu, ²¹ 14 nonlinear problems whose number of iterations are smaller, 13 nonlinear problems whose total steps are smaller. All the nonlinear equations are solved accurately, the inverse BFGS method of this paper has been verified.

Table 3.

Calculating steps comparison.

Step/Comparison	Smaller	Equal	Larger
BKT	8	4	5
NI	14	0	3
TS	13	1	3

Numerical experience of multi-rigid body

A swing problem of two rigid bodies connected by a driving hinge and a rotating hinge is studied. Rigid body A with steel material is connected to the ground by a driving hinge at the left side, as shown in Figure 4. The hinge m drives body A rotating about z axis with the speed of 60 °/s. Rigid body B is connected to the right side of body A by a rotating hinge n. And body B can rotate freely about the rotating shaft. The mass of body A is 336.14 kg, moments of inertia $I_x$ , $I_y$ , $I_z$ around the center of mass are 1.38 kgm², 131.3 kgm² and 130.5 kgm², respectively. The length, width and thickness of body A are 2.0 m, 0.2 m and 0.1 m respectively. Materials and sizes of rigid body B are the same with body A.

Figure 4.

Rigid body system.

Initial positions of the rigid body system are shown in Figure 4, where body A and body B are placed horizontally. The dynamic responses are computed with the method of this paper, and then the results are compared with those of Software Adams. Dynamic responses of displacement of body A are shown in (a) and (b) of Figure 5, while dynamic responses of displacement of body B are shown in (c) and (d) of Figure 5.

Figure 5.

Dynamic responses of displacement of rigid body systems. (a) mass center of body A in x direction, (b) mass center of body A in y direction, (c) mass center of body B in x direction, (d) mass center of body B in y direction.

The dynamic problem of the two multi-rigid body is solved using the DAE method of this paper in which both kinetic equations of the two bodies and constraint equations of the two joints are involved. Dynamic responses of velocity of body B are shown in (c) and (d) in Figure 6, and the rotating angle responses are compared in Figure 7. It is shown that the results of these paper fit well with those of software Adams. The validities of hinge models, rigid models and the predictor-corrector method of this paper are verified.

Figure 6.

Dynamic responses of velocity of rigid body B. (a) mass center of body B in x direction, (b) mass center of body B in y direction.

Figure 7.

Dynamic responses of rotating angle of rigid body B.

Numerical experience of rigid-flexi body

Elastic blades of wind turbine NH1500 ³³ rotate about rotational axis of hub when wind turbine operates. Nacelle is simplified as a rigid body with mass of 6000 kg, hub is simplified as beams to connect the three blades. Finite element model of turbine rotor is construed, modes of rotor are analyzed. Frequencies of flap-wise mode and edge-wise mode are presented in Table 4. The corresponding modal shapes are shown in Figures 8 to 10.

Table 4.

Frequencies of NH1500 rotor.

Modal numer	Modal type	Frequency/Hz
First order of flap-wise mode	1. Antisymmetry	1.1609
	2. Antisymmetry	1.1609
	3. Symmetry	1.1609
Second order of flap-wise mode	4. Antisymmetry	3.102
	5. Antisymmetry	3.102
	6. Symmetry	3.1023
Third order of flap-wise mode	7. Antisymmetry	6.0375
	8. Antisymmetry	6.0375
	9. Symmetry	6.0383
First order of edge-wise mode	10. Antisymmetry	1.4754
	11. Antisymmetry	1.4757
	12. Symmetry	1.4757
Second order of edge-wise mode	13. Antisymmetry	4.7204
	14. Antisymmetry	4.7213
	15. Symmetry	4.7213

Figure 8.

First order of flap-wise modal shapes of wind rotor. (a) antisymmetric, (b) antisymmetric, (c) symmetric.

Figure 9.

First order of edge-wise modal shapes of wind rotor. (a) antisymmetric, (b) antisymmetric, (c) symmetric.

Figure 10.

Second order of flap-wise modal shapes of wind rotor. (a) antisymmetric, (b) antisymmetric, (c) Symmetric.

Rotor-nacelle system is constructed by elastic blades, a rigid hub, a rigid nacelle and a generator. The nacelle is fixed to the ground at bottom, the generator is simplified as a driving motion with rotation speed of 17.2 r/min, as shown in Figure 11.

Figure 11.

Rotor-nacelle system.

Nacelle and rotor are modeled as a rigid body and flexible bodies respectively by multibody dynamics. A driving hinge is built with a speed of 17.2 r/min to represent the rotation with constant speed. Meanwhile, the nacelle is locked at bottom and thus a nacelle-rotor system is constrained. Blades are subjected to the combined action of inertia forces, elastic forces and aerodynamic forces in the process of rotation, where the inertial forces mainly come from gravity, Coriolis force and centrifugal force. Aerodynamic forces are calculated by free vortex wake method and then applied to blade structure at each blade section by multipoint method. ³⁴ Structural dynamic responses are calculated by multibody method of this paper.

On the other hand, to verify the method of this paper, the rotor coordinate system is simplified as inertial system because of the low rotational speed. Thus, the body motions are simplified as vibration process excited by gravity forces, elastic forces and aerodynamic forces in the inertial coordinate system. Then structural dynamic responses are calculated using modal method and compared with the results of multibody method of this paper. During the calculation, aerodynamic forces are calculated in advance and applied to blade structure afterwards. Vibration responses at 40.1 m of a blade in rotor coordinate system by the two methods are compared as shown in Figure 12. It is shown that responses of this paper are coincident with those of modal method generally, thus multibody method of this paper is verified. It is implied that modal method can be used to calculated structural responses approximately. Dynamic responses of generalized coordinates of rotor are shown in Figure 13. It is also implied that dynamic responses of symmetric models can represent the averaged static deformations and dynamic responses of antisymmetric models represent the vibration around averaged static deformations. However, rotor coordinate system is not inertial coordinate system actually. There are deviations in the calculated results of the two methods, and Coriolis force and centrifugal force should be considered. The method of this paper is suitable for sophisticated dynamics simulation of large-scaled wind turbine blade.

Figure 12.

Displacemental responses at a blade tip. (a) flap-wise displacement, (b) edge-wise displacement.

Figure 13.

Dynamic responses of generalized coordinates of the turbine rotor by multibody method.

Conclusion

A solver with derivative-free for multibody dynamics is developed based on a derivative-free line search. In the method, a line search method is combined with inverse BFGS method, a nonlinear equations method with derivative-free is established. The inverse of Jacobian matrix is estimated instead of solving the linearized equations of this paper, it is suitable for problems whose Jacobian matrix is hard to get. A series of nonlinear equations are solved using the nonlinear equation method, its practicability is verified. Furthermore, the nonlinear equation method is integrated into multibody dynamics method, thus a new multibody dynamics solver is proposed. Obviously, the Jacobian matrix of dynamic system need not be calculated. The proposed multibody method can be used to those problems whose system Jacobian matrix cannot be obtained. Also, a multi-rigid-body system and a rigid-flexible-body system are simulated. The practicability of proposed multibody dynamics method is also verified. It is also implied that dynamic responses of flexible turbine rotor can be divided into two parts. Symmetric models can represent the averaged static deformation part and antisymmetric models can represent the vibration part around averaged static deformations.