Vincent Monchiet
email: vincent.monchiet@ed.univ-lille1.fr

Eric Charkaluk
email: eric.charkaluk@univ-lille1.fr

Djimedo Kondo
email: kondo@univ-lille1.fr

An improvement of Gurson-type models of porous materials by using Eshelby-like trial velocity fields

Keywords: Computational solid mechanics, Ductile porous metals, Generalized Gurson criterion, Prolate and oblate voids, Eshelby velocity fields, Limit analysis Mécanique des solides numériques, Métaux ductiles poreux, Critère de Gurson généralisé, Cavités allongées et aplaties, Champ de vitesse d'Eshelby

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Since three decades, the modelling of the behavior of ductile porous media has been the subject of important researches in non linear mechanics of materials. In his pioneering work, Gurson [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criterion and flow rules for porous ductile media[END_REF] developed a limit analysis approach of a hollow sphere. The plastic matrix is assumed to obey the von Mises criterion: f(σ ) = σ eq -σ 0 0; with: σ eq = 3 2 σ : σ

where σ denotes the microscopic stress field, σ its deviatoric part and σ eq the microscopic von Mises equivalent stress. σ 0 represents the yield stress in tension. More specifically, Gurson obtained a macroscopic criterion which, in the case of spherical voids, reads:

Σ 2 eq σ 2 0 + 2f cosh 3Σ h 2σ 0 -1 -f 2 = 0( 2 )
Σ h denotes the hydrostatic stress and f the porosity. Σ eq is the macroscopic von Mises equivalent stress. It has been demonstrated (see for instance [START_REF] Perrin | Contribution à l'Étude Théorique et Numérique de la Rupture Ductile des Métaux[END_REF]), that the yield surface defined by ( 2) constitutes an upper bound for Hashin's wellknown composite spheres assemblage and gives the exact result for purely hydrostatic macroscopic loading. Later, using variational techniques, Ponte-Castañeda [START_REF] Ponte-Castañeda | The effective mechanical properties of nonlinear isotropic composites[END_REF] and Suquet [START_REF] Suquet | On bounds for the overall potential of power law materials containing voids with an arbitrary shape[END_REF] obtained a rigorous nonlinear Hashin-Shtrikhman upper bound which, for spherical voids, takes the following form1 :

1 + 2 3 f Σ 2 eq σ 2 0 + 9 4 fΣ 2 h -(1 -f) 2 = 0( 3 )
An important observation is that the Gurson model (Eq. ( 2)) violates this upper bound for low values of the stress triaxiality T = Σ h /Σ eq . However, its predominance over the Hashin-Shtrikhman bound is still observed for high stress triaxialities. A possible method to improve the predictions of the original limit analysis approach of Gurson and in fact the subsequent models consists in considering refined trial velocity fields (see for instance [START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF] and [START_REF] Garajeu | A micromechanical approach of damage in viscoplastic materials by evolution in size, shape and distribution of voids[END_REF]). Still, due to the limitation of trial velocity fields which have been explored, few significant improvements have been obtained in the past studies. Therefore, the main objective of the present Note is to develop a limit analysis approach based on Eshelby-like velocity fields and to derive new expression of the yield function. The calculations will be performed in the general case of prolate and oblate voids (Sections 2 and 3). In this way, it is expected that the new approach will also provide an improved version and a generalization of the results obtained by [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric prolate ellipsoidal cavities[END_REF] and [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric oblate ellipsoidal cavities[END_REF]. In order to provide simple illustrations of the obtained approximate criterion, some specific cases such as spherical or cylindrical voids will be examined (Section 4).

Basic concepts and methodology

The studied cell

As in [START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF], consider a spheroidal (axisymmetric) prolate or oblate cavity with semi-axes a 1 (along e 3 ), and b 1 (along e 1 and e 2 ) embedded in a cell which has the shape of a confocal spheroid with the semi-axes a 2 (along e 3 ), and b 2 (along e 1 and e 2 ). a 1 >b 1 corresponds to a prolate cavity while b 1 >a 1 is associated to an oblate one. Let us denote c the focal distance and e 1 the eccentricity defined by: with L λ = a 2 sin 2 (β) + b 2 cos 2 (β), θ ∈[0, 2π ], β ∈[0,π] and e ρ = cos(θ )e 1 + sin(θ )e 2 .

c = a 2 1 -b 2 1 ; e 1 = c a 1 (prolate); c = b 2 1 -a 2 1 ; e 1 = c b 1 (oblate) (4)

The Eshelby-like velocity field

The velocity field, v, in the matrix, is classically decomposed into a uniform velocity field, A.x, and an heterogeneous velocity field, v E , as follows: v = A.x + v E .Forv E , we consider the exterior point Eshelby solution (see [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problem[END_REF]) adapted here to an incompressible viscous fluid containing a spheroidal inclusion. For convenience, this solution (see also [START_REF] Mura | Micromechanics of Defects in Solids[END_REF]) can be put in the form:

v E = r=6 r=1 v r d * r (6)
where

d * r = d * : Q r /(Q r : Q r ).
Tensor d * is an eigenstrain in the inhomogeneity and Q r are defined by:

Q 1 = 1; Q 2 = 1 -3e 3 ⊗ e 3 ; Q 3 = e 2 ⊗ e 2 -e 1 ⊗ e 1 ; Q 4 = e 1 ⊗ e 2 + e 2 ⊗ e 1 Q 5 = e 1 ⊗ e 3 + e 3 ⊗ e 1 ;
Q 6 = e 2 ⊗ e 3 + e 3 ⊗ e 2 [START_REF] Garajeu | A micromechanical approach of damage in viscoplastic materials by evolution in size, shape and distribution of voids[END_REF] with the property Q r : Q s = 0ifr = s. The quantities d * r are then related to the components of d * by:

d * 1 = 1 3 (d * 11 + d * 22 + d * 33 ); d * 2 = 1 3 d * 11 + d * 22 2 -d * 33 d * 3 = d * 22 -d * 11 2 ; d * 4 = d * 12 ; d * 5 = d * 13 ; d * 6 = d * 23 (8)
The velocity fields v r for r = 1, 6 are given in Appendix A, in the spheroidal frame. Note that v 1 and v 2 are two axisymmetric velocity fields; they are independent of the coordinate θ and the component v θ = 0. The first field v 1 is the one used by [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric prolate ellipsoidal cavities[END_REF] and by [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric oblate ellipsoidal cavities[END_REF] for the determination of the macroscopic yield function of a plastic matrix containing prolate and oblate voids respectively. Due to some limitations of v 1 , it has been proposed in [START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF] to incorporate a supplementary velocity field which does not appear in the Eshelby velocity field.

Garajeu and Suquet [START_REF] Garajeu | A micromechanical approach of damage in viscoplastic materials by evolution in size, shape and distribution of voids[END_REF] have also derived an expression for the macroscopic yield function of the porous material in the case of prolate cavities by using a truncated expression of the velocity fields v 1 and v 2 which corresponds to f 1 (λ) = g 1 (λ) = 0 (see Appendix A). It is also interesting to mention that v 1 and v 2 are contained in the general axisymmetric velocity fields proposed by Lee and Mear [START_REF] Lee | Axisymmetric deformation of power-law solids containing a dilute concentration of aligned spheroidal voids[END_REF]. However, the other fields v r for r = 3, 6, considered in this study, are nonaxisymmetric and are not contained in the Lee and Mear fields.

The microscopic plastic strain rate is defined by

d = A + r=6 r=1 d r d * r in which d r =∇ s v r .
The microscopic dissipation reads π(d) = σ 0 d eq where d eq = 2 3 d : d is the equivalent plastic strain rate defined by:

d 2 eq = A 2 eq + 2 r=6 r=1 d * r A : d r + r=6 r=1 s=6 s=1 d * r d * s d r : d s (9)
Since the Eshelby-like velocity field introduced by ( 6) does not comply with the uniform strain on the cell boundary, we propose to use the classical average rule which relates the macroscopic strain rate tensor D to d. This allows us to relate d * to D by:

D = 1 |Ω| Ω d dV = A + f S(e 2 ) : d * (10)
where |Ω|=4πa 2 b2 2 /3 denotes the volume of the studied cell (matrix + void) and S(e) is the Eshelby tensor corresponding to a spheroidal cavity of eccentricity e embedded in an incompressible medium; the components of S(e) are expressed in Appendix B. By combining then (B.1) and [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problem[END_REF], one obtains:

d * 1 = D h f ; A 33 = D 33 (1 -α 2 ) -(D 11 + D 22 )α 2 + εf 2α 2 a 2 2 + (α 2 -1)b 2 2 2c 2 d * 2 A 22 -A 11 = D 22 -D 11 -εf 2a 2 2 + 3b 2 2 (α 2 -1) 4c 2 d * 3 ; A 12 = D 12 -εf 2a 2 2 + 3b 2 2 (α 2 -1) 8c 2 d * 4 ( 11
)
A 13 = D 13 -εf a 2 2 + b 2 2 4c 2 (1 -3α 2 )d * 5 ; A 23 = D 23 -εf a 2 2 + b 2 2 4c 2 (1 -3α 2 )d * 6
in which D h 1 represents the spherical part of D and the convention ε = 1 for prolate voids and ε =-1 for oblate voids is adopted. α 2 = α(e 2 ) where α(e) is function of eccentricity e and is given by (A.2).

It is readily seen that the complete velocity field is defined by 11 parameters, the components of A and the scalar d * r for r = 1, 6. Condition [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problem[END_REF], detailed in [START_REF] Mura | Micromechanics of Defects in Solids[END_REF] together with (A.2), gives six relations between the different parameters. Still, their remains five unknown parameters which have to be determined. 2 These are the d * r for r = 2, 6 or equivalently the components of d * (the deviatoric part of d * ). It is interesting to note that, by putting d * r = 0forr = 2, 6, and choosing A axisymmetric in [START_REF] Mura | Micromechanics of Defects in Solids[END_REF], the homogeneous boundary strain rate conditions obtained by [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric prolate ellipsoidal cavities[END_REF][START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric oblate ellipsoidal cavities[END_REF] are recovered.

The macroscopic dissipation Π(D) and the minimization principle

Due to the dependence of the velocity field d with the following unknown kinematic parameters d * , let us now introduce Π(D, d * ) defined by:

Π D, d * = σ 0 |Ω| Ω-ω d eq dV (12)
where ω denotes the volume of the void. It is readily understood that the searched expression of the macroscopic dissipation, Π(D), is derived from a minimization procedure of Π(D, d * ) with respect to d * :

Π(D) = min d * Π D, d * (13)
The yield surface, related to the macroscopic dissipation, is then assumed to be given by:

Σ = ∂Π ∂D (14)
In fact such consideration can be interpreted as the condition of the necessary coherence of the energy definition at the two scales. Due to the difficulty of integrating [START_REF] Lee | Axisymmetric deformation of power-law solids containing a dilute concentration of aligned spheroidal voids[END_REF], some approximations are needed either for prolate and oblate voids in order to determine the macroscopic dissipation and then the macroscopic yield function.

Approximate expression of Π(D, d * ) and determination of the macroscopic yield locus

We now aim to derive the expression of Π(D,d * r ) by using the same type of approximations as in [START_REF] Monchiet | Approximate yield criteria for anisotropic metals with non spherical voids[END_REF]. Let us first recall the expression of Π(D,d * r ) in the spheroidal frame.

Π D, d * = σ 0 |Ω| λ=λ 2 λ=λ 1 β=π β=0 θ =2π θ =0 d eq bL 2 λ sin β dλ dβ dθ (15)
The following approximations are made:

A1: d eq is replaced by a mean value { d 2 eq E } 1/2 on each confocal ellipsoid E.

It follows that (15) takes the following form:

Π D, d * = σ 0 a 2 b 2 2 λ=λ 2 λ=λ 1 d 2 eq E 1/2 b 2a 2 + b 2 dλ ( 16
)
where d 2 eq E is such that:

d 2 eq E = 3 4π(2a 2 + b 2 ) β=π β=0 θ =2π θ =0 d 2 eq L 2 λ sin(β) dβ dθ (17)
Then, from the definition ( 9), one has:

d 2 eq E = A 2 eq + 2 r=6 r=1 d * r A : W r (e) + r=6 r=1 s=6 s=1 d * r d * s P rs (e) (18)
in which W r (e) = 2 3 d r E and P rs (e) = 2 3 d r : d s E . Note that W r (e) and P rs (e) also depend on e 1 . Let us now introduce the following variables x and y defined by:

x = a 1 b 2 1 ab 2 ; y = χa 1 b 2 1 c 3 + χab 2 (19)
with χ = 3 4 π 2 + 32 3 . As a unified notation, let us introduce u such that u = x for prolate cavities and u = y for oblate cavities. A2: Following [START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF], the expressions of W r (e) and P rs (e) are replaced by w r u 2 and p rs u 2 respectively where w r and p rs are constant.

As in [START_REF] Monchiet | Approximate yield criteria for anisotropic metals with non spherical voids[END_REF], parameters p rs and tensors w r are respectively determined as the mean values of P rs (e)/u 2 and W r (e)/u 2 along the interval [u 1 ,u 2 ]:

p rs = 1 u 2 -u 1 u 2 u 1 P rs (e) du u 2 ; w r = 1 u 2 -u 1 u 2 u 1 W r (e) du u 2 ; for r = 1, 6 (20)
u 1 = 1 for a prolate cavity and

u 1 = χa 1 b 2 1 /(c 3 + χa 1 b 2 1
) for an oblate one. This approximation is consistent only if the obtained values for p rs and w r show low variations according to e 1 and e 2 . The explicit determination of p rs and w r is performed by using Maple software; the resulting expressions are given in Appendix C. The quantity d 2 eq E can be rewritten as follows: where W and P are two fourth order tensors defined in Appendix C. Eq. ( 21) can also be expressed as follow:

d 2 eq E = A 2 eq + 2
d 2 eq E = A 2 eq -A : W : P -1 : W T : Au 2 + d * + A : W : P -1 : P : d * + P -1 : W T : A u 2 (22)
where W T represent the transpose of W such that W T ij kl = W klij .

A3: Term A : W : P -1 : W T : Au 2 is neglected.

Although this term is not neglected in [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric oblate ellipsoidal cavities[END_REF], [START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF] and [START_REF] Monchiet | Approximate yield criteria for anisotropic metals with non spherical voids[END_REF], its effect on the macroscopic criterion is generally weak3 and motivates the approximation A3. Let us then introduce the following change of variable:

B 2 = B : P : B; with: B = d * + P -1 : W T : A (23)
Π(D, d * ) (Eq. ( 16)) takes then the following form:

Π D, d * = σ 0 f u 1 u 2 A 2 eq + B 2 u 2 1/2 du u 2 = σ 0 f B arcsinh uB A eq - A 2 eq + u 2 B 2 u u 1 u 2 (24)
As previously indicated, the macroscopic yield locus is defined by [START_REF] Ponte-Castañeda | Nonlinear composites[END_REF] in which Π(D) corresponds to the minimum of Π(D, d * ) with respect to d * , d * 1 being completely identified and related to D h (see [START_REF] Mura | Micromechanics of Defects in Solids[END_REF]). The macroscopic criterion is then obtained by resolution of:

Σ = ∂ Π(D, d * ) ∂D with: ∂ Π(D, d * ) ∂ d * = 0 (25)
An analytical solution to this system, given in Appendix D, leads to the following generalized approximate criterion corresponding to prolate and oblate voids4 :

Σ 2 eq σ 2 0 + 2(1 + g)(f + g)cosh Σ P σ 0 -(1 + g) 2 -(f + g) 2 = 0 ( 26
)
where g is conventionally taken as 0 for a prolate void and is given by g = c 2 /(χa 2 b 2 2 ) for an oblate one. Σ P and Σ eq which enters in the criterion (26) are defined by:

Σ 2 P = Σ : S(e 2 ) : P -1 : S T (e 2 ) : Σ; Σ 2 eq = 3 2 Σ : Σ; with: Σ = Σ -W : P -1 : S(e 2 ) : Σ (27)
where the components of P -1 are given in Appendix C. Eqs. ( 26) and ( 27) constitute the most important results of the present Note. Their new features and advantage over the existing Gurson-type criteria of porous media (see [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric prolate ellipsoidal cavities[END_REF] and [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric oblate ellipsoidal cavities[END_REF]) come mainly from the terms introduced by the components of d * . In particular, it is interesting to point out that, due to Σ P , the new criterion introduces a coupling between pure shear stresses and the porosity f which does not exist in the existing criteria. Note also that the nonlinear Hashin-Shtrikman bound exhibits a coupling of this type (see [START_REF] Ponte-Castañeda | Nonlinear composites[END_REF]).

The particular cases of spherical and cylindrical voids

Let us consider, as a first particular illustration, the case of the spherical cavity; this is obtained by putting e 1 → 0 and e 2 → 0. It follows that W = 0, and then Σ eq = Σ eq . Coefficients p rs are p 11 = 4, p 12 = 0, p 22 = 3p 33 = 3p 55 = 24/25, and then Σ 2 P = 9Σ 2 h /4 + 2Σ 2 eq /3. The resulting macroscopic yield locus5 is given by:

Σ 2 eq σ 2 0 + 2f cosh 1 σ 0 9 4 Σ 2 h + 2 3 Σ 2 eq -1 -f 2 = 0 (28)
In the particular loading case Σ eq = 0, the proposed criterion (28) reduces to the Gurson one which gives the exact solution under purely hydrostatic pressure. Eq. ( 28) provides also the strongly remarkable property that the deviatoric equivalent stress Σ eq enters with the mean stress Σ h in the cosh term which controls the cavity growth process.

Consider now low values of Σ h ; it can be verified that in this case cosh{ 1

σ 0 9 4 Σ 2 h + 2 3 Σ 2 eq }≃1 + 9 8 Σ 2 h σ 2 0 + 1 3 Σ 2 eq σ 2 0 .
Reporting this result into (28), one recovers the Hashin-Shtrikman bound (see Eq. ( 3)). The predictions of the yield surface by the different models, in the case of spherical voids and for two different porosities, are represented on Fig. 1 and confirm the good performance of the proposed model. It must be recalled that, unlike the studies [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criterion and flow rules for porous ductile media[END_REF] and [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric prolate ellipsoidal cavities[END_REF][START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric oblate ellipsoidal cavities[END_REF], the results obtained in the present study are based on a nonkinematical velocity field defined by the condition [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problem[END_REF]; it is this choice which allows the connection with the nonlinear homogenization bounds [START_REF] Ponte-Castañeda | Nonlinear composites[END_REF].

Let us consider now, as a second illustration, the case of a cylindrical cavity; then e 1 → 1, e 2 → 1, W = 0, then Σ eq = Σ eq . Coefficients p rs are p 11 = 3, p 12 = p 22 = 0 and p rr = 1/3f o rr = 3, 6, then Σ 2 P = (2Σ 2 11 + 2Σ 2 22 + 4Σ 2 12 + 4Σ 2 13 + 4Σ 2 23 )/3. The macroscopic yield surface is described by:

Σ 2 eq σ 2 0 + 2f cosh 1 σ 0 3 2 Σ 2 11 + Σ 2 22 + 3Σ 2 12 + 3Σ 2 13 + 3Σ 2 23 -1 -f 2 = 0 (29)
As for spherical voids, similar comments can be done for cylindrical cavities: the exact solution under purely hydrostatic pressure is recovered; for low values of Σ h /σ 0 the criterion leads to the approximation:

Σ 2 eq σ 2 0 + 3f 2 1 σ 2 0 Σ 2 11 + Σ 2 22 + 2Σ 2 12 + 2Σ 2 13 + 2Σ 2 23 -(1 -f) 2 = 0 (30)
which coincides with the result established by Suquet [START_REF] Suquet | On bounds for the overall potential of power law materials containing voids with an arbitrary shape[END_REF].

with:

f 1 (λ) = 1 -3α; g 1 (λ) =- 3 4ab 2a 2 α + b 2 (α -1) ; k 2 = ε 2a 2 1 + b 2 1 2c 2 α(e) = ab 2 c 3 arctanh c a - b 2 c 2 (prolate),α ( e ) =- ab 2 c 3 arctan c a + b 2 c 2 (oblate) (A.2)
The two next velocity fields are given by: with:

v 3 = a 1 b 2 1 bL λ f 3 (λ) sin 2 (β)e λ + 1
f 3 (λ) = 1 + k 3 b 2 (3α + 1)b 2 -2a 2 ; g 3 (λ) = k 3 ab 2a 2 + 3b 2 (α -1) ; k 3 = ε b 2 1 4c 2
The last velocity fields are:

v 5 = a 1 b 2 1 bL λ f 5 (λ) sin(2β)cos(θ )e λ + a 1 b 2 1 bL λ g 5 (λ) 1 + k 5 cos(2β) cos(θ )e β - a 1 b 2 1 b 2 G 5 (λ)(1 + k 5 ) cos(β) sin(θ )e θ v 6 = a 1 b 2 1 bL λ f 5 (λ) sin(2β)sin(θ )e λ + a 1 b 2 1 bL λ g 5 (λ) 1 + k 5 cos(2β) sin(θ )e β + a 1 b 2 1 b 2 g 5 (λ)(1 + k 5 ) cos(β) cos(θ )e θ (A .4)
with:

f 5 (λ) = 1 4ab (3α + 1)(1 + k 5 )a 2 + 3(1 -α)(1 -k 5 )b 2 ; g 5 (λ) = 1 2 (3α -
S 3333 = ε c 2 3αa 2 -b 2 ; S 1133 (e) = S 2233 (e) = (1 -3α) εa 2 1 2c 2 S 3311 (e) = S 3322 (e) = (1 -3α) εb 2 1 2c 2 ; S 2323 (e) = S 1313 (e) = ε(1 -3α) a 2 1 + b 2 1 4c 2 (B.1)
Let us recall that α is function of eccentricity e and is given by (A.2).

Appendix C. Expression of tensor W and P

Q r : P : Q s = p rs for r = 1, 6; with:

p rs = (1 + g)(f + g) f(1 -f) H rs (e 1 ) -fH rs (e 2 ) (C.1)
where the nonnull H rs (e) are:

H 11 = 3(1 + 3α)(1 -α); H 12 = 3(1 -3α)(1 -α -β); H 22 = 3(3α + 3β -1)(1 -α -β) H 33 = 1 12 (1 + 3α -β)(3 -3α + β); H 55 = 1 3 (1 -3α -2β)(3α -3 + 2β) (C.2)
Tensor W is given by:

W = 2 3 (1 + g)(f + g) f(1 -f) S(e 2 ) -S(e 1 ) (C.3)
The nonzero components of tensor P are: The aim of the following demonstration is to put (D.1) in the form:

P 1111 = P 2222 = 1
Σ eq = ∂ Π ∂A eq ; Σ P = 1 f ∂ Π ∂B (D.2)
where Σ eq and Σ P have to be determined. Considering (24), Σ eq and Σ P , defined by (D.2) read: with J = 1 3 1 ⊗ 1 and K = I -J, I being the symmetric fourth-order identity tensor and 1 the second-order identity tensor.

Σ eq = ∂ Π ∂A eq = σ 0 f 1 + f 2 ξ 2 f -1 + ξ 2
For the determination of Σ eq and Σ P , it is convenient to expressed Σ : S(e 2 ) from the first relation into (D.6), or more precisely its deviatoric part and its hydrostatic part Σ : S(e 2 ). Using the following properties 1 : S(e

  It is convenient to introduce the spheroidal coordinates characterized by λ, β, θ, and defined in the cylindrical frame (coordinates ρ,θ,z)b yρ = b sin β, z = a cos β. The iso-λ surface defines confocal spheroids, with semi-axes a = c cosh(λ), b = c sinh(λ) and eccentricity e = c/a, for prolate. An oblate spheroid is associated to semi-axes a = c sinh(λ), b = c cosh(λ) and eccentricity c/b. The unit vectors of the new base are: e λ = 1 L λ a sin(β)e ρ + b cos(β)e z ; e β = 1 L λ b cos(β)e ρ -a sin(β)e z ; e θ = e θ (5)

s u 2

 2 = A 2 eq +[2A : W : d * + d * : P : d * ]u 2 (21)

Fig. 1 .

 1 Fig. 1. Yield surface of the porous material-case of spherical void: comparison of the proposed model (28), the Gurson yield locus and the Ponte-Castañeda model. (a) Porosity f = 0.01; (b) porosity f = 0.3.

Fig. 1 .

 1 Fig. 1. Surface de charge des matériaux poreux-cas d'une cavité sphérique : comparaison du modèle proposé (28) avec la surface de Gurson et avec le modèle de Ponte-Castaneda. (a) Porosité f = 0,01 ; (b) porosité f = 0,3.

2 g 3 b 2 1 b 2 4 = a 1 b 2 1 bL λ f 3 2 g 3 a 1 b 2 1 b 2

 3124132312 (λ) sin(2β)e β cos(2θ)-a 1 g 3 (λ) sin(β) sin(2θ)e θ v (λ) sin 2 (β)e λ + 1 (λ) sin(2β)e β sin(2θ)+ g 3 (λ) sin(β) cos(2θ)e θ (A.3)

2 3 S

 23 Appendix B. Components of S(e) for a spheroidal inclusion with axis along e 1111 (e) = S 2222 (e) = 3S 1122 (e) = 3S 2211 (e) = 3S 1212 (e) = 3ε 2a 2 + 3(α -1)b 2 8c 2

(C. 4 )

 4 from which are easily obtained the components of P -1 Appendix D. Determination of the macroscopic yield locus Let us introduce the following change of variables (D, d * ) ⇒ (A,B) in (25):

eq σ 0 2 + 2 ( 1

 221 0 arcsinh(f ξ ) -arcsinh(ξ ) (D.[START_REF] Ponte-Castañeda | The effective mechanical properties of nonlinear isotropic composites[END_REF] with ξ = B/A eq . Eliminating ξ in relations (D.3), we find:Σ + g)(f + g)cosh Σ P σ 0 -(1 + g) 2 -(g + f) 2 = 0( D . 4 )Recalling that the expression of tensor A is given by[START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problem[END_REF] and the one of tensor B by (23), one has:A = Df S(e 2 ) : d * ; B = d * + P -1 : W T : Df S(e 2 ) : d * ; with: d * = ∂A : I -S(e 2 ) : J + 1 fB ∂ Π ∂B B : P : J + f P -1 : W T : I -S(e 2 ) : J ∂ Π ∂A : S(e 2 ) : K = 1 fB ∂ Π ∂B B : P : Kf P -1 : W T : S(e 2 ) : K (D.6)

  2 ) : 1 = 3 and J : S(e 2 ) : K = 0, one obtains: ⎧ ⎨ ⎩ Σ : S(e 2 ) : 1 = 1 fB ∂ Π ∂B B : P : 1 Σ : S(e 2 ) : K = ∂ Π ∂A : S(e 2 ) : K + 1 B ∂ Π ∂B B : W T : S(e 2 ) : K (D.7)

Similar results for fluid saturated porous media with a Drucker-Prager matrix can be found in[START_REF] Dormieux | Microporomechanics[END_REF].

Note that the components of A could also be chosen as unknowns but this does not change the final results.

However, note that the consideration of these terms does not rise an important computational difficulty in the determination of the macroscopic yield locus but gives more complex expressions.

Obviously, it will be verified that this result generalizes also the Gurson criterion.

Note that this result extends to fluid saturated porous media by replacing Σ h by Σ h + P where P is the fluid pressure.

Acknowledgements

The authors are grateful to Professor J.-B. Leblond for fruitful suggestions on Note [13] which are partly at the origin of the present study.

Appendix A. Eshelby velocity field

The exterior point Eshelby solution associated to a spheroidal inclusion in an isotropic matrix, can be found in [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problem[END_REF] or in the textbook [START_REF] Mura | Micromechanics of Defects in Solids[END_REF]; the corresponding velocity field is given as a combination of the following velocities: which leads to the expression of Σ P (27). We now come to the determination of Σ eq . Considering the deviatoric part of Σ from (D.6), and noticing that W T : K = W T ,wehave:

Replacing B by its expression deduced from (D.9), one obtains:


Hassan Adeyoola

Advanced Persistent Threat: Detection and Defence

Keywords: cyber-criminal, organization, Advanced Persistent Threat, undercover activities (keywords)

as the growth and popularity of technology has become simultaneous ascend in both impacts and numbers of cyber criminals thanks to the web. For many years, the organization has strived in ways of preventing any attacks from cyber-criminal with advanced techniques. Cybercriminals and intruders are developing a more advanced way to breach the security surface of an organization. Advanced Persistent Threats are also known as APT are new and a lot more sophisticated version for multistep attack scenarios that are known and are targeted just to achieve a goal most commonly undercover activities. this report, there will cover everything I know that tells us about APT with more word and brief explanations

I. INTRODUCTION

Since the invention of the internet, we have been relying on networks and servers to help function our lives of 1.4 billion people (and counting) on the planet earth. Unfortunately, in the world we live in many malicious activities are going on, some report the lately and their damage has a massive effect on their victim's lives. From stealing information to masquerading as someone authentic, these are cybercrimes, and they are a danger in our daily lives and this report, I will go profoundly on a topic that correlates to this type. I will adhere to research that I have found on the internet about this topic. This will also include the common types of mitigating cyber-attacks which means the detections and defense that are used. What is also included in this report will show how some sites progressions on extracting data malware that cybercriminals use including their location just to show insight on how the world is fighting back against cybercriminals.

II. WHAT IS APT

"An advanced persistent threat (APT) is a broad term used to describe an attack campaign in which an intruder, or team of intruders, establishes an illicit, long-term presence on a network in order to mine highly sensitive data." [START_REF] Berros | 5 Most Common Types Of APT Attack Vectors | XM Cyber[END_REF].

In the cyber security business, an APT is s broad term that is used to describe an attack campaign especially when an intruder or a group of intruders, set up an unlawful longterm bearing on a network for the sake of mining highly important data.

III. TYPES OF APT ATTACKS

The most crucial thing to comprehend about APTs is their flexibility. Although they are capable of infiltrating attacks of tremendous and sophistication, they can also indulge in a very basic attack. Rivals are concerned about the ROI also known as the return of investment and the efficiency just like everyone in an organization and it is quite often that a single attack is successful. Below are some of the most common vectors for the attacks and they are:

A. Phishing

This is known to be a type of social engineering which mostly consists of stealing user data, which includes the user's credit card numbers and login credentials [START_REF] Ghafir | A Survey on Botnet Command and Control Traffic Detection[END_REF]. It happens when an attacker, masquerades themselves as a trusted or authentic entity, then deceives the victim into opening a text, email, or instant message.

B. DNS modifications

Also known as dns hijacking. It consists of manipulating transactions and user will be deceived since they are not aware of the activity of their servers especially during an internet session. Since it is malicious exploitation, this means users are being redirected with the help of DNS servers that are rouge and never change the IPS address of the user that is being redirected of their servers especially during an internet session.

C. Zero days attacks

Which is also described as Day Zero, is an attack that involves the exploitation of a possibly a serious software security weakness in which the developers or vendors are not aware of [START_REF] Ghafir | Advanced Persistent Threat and Spear Phishing Emails[END_REF][START_REF]What Is APT (Advanced Persistent Threat)[END_REF]. This means that the vendors or security software developer must quickly resolve the weakness immediately when it has been discovered as this could limit the threat to those who uses the software.

D. Supply chain attacks

A supply chain attack also known as either third party or value chain attack, it happens when an intruder penetrates your system through a provider or an outside partner that can have access to an organisation system [START_REF] Ghafir | Advanced Persistent Threat Attack Detection: An Overview[END_REF]. This is known to have significantly change enterprises' attack surface for the past few years, since there have been more service providers and suppliers grabbing their hands-on important data more than ever before E. Ransomware First, ransomware is a type of mischievous software that is designed to stop the access to a computer files or even the computer until the victim paid the sum of money. Most ransomware alternatives are encrypting files on the computer affected, that will make them inaccessible, and demand a ransom payment to restore access to the computer.

F. Pirated software Also known as software piracy, this is when an unapproved copying of a software that has been purchased [START_REF] Ghafir | DNS query failure and algorithmically generated domain-flux detection[END_REF][START_REF] Eyers | Sabotage Of Critical Infrastructure A Growing Risk, Says Former Chief Spy[END_REF]. This means when you buy the software from a company it does not mean you are the owner of that software; all it means is that you are a licensed user.

IV. RISK AND IMPACTS

The victims of these types of APT that are especially the main target are being carefully researched and illegally monitored. They are mostly large industrial business or local networks. These are known to be mostly vulnerable and the impact it could are known to be following and they are :

A. Rational property theft

This known to include patents or the secrets for trading. The impact goes deeper according to this website it says that these sorts of situations keep chiefs up around evening time all things considered: Intellectual property (IP) which is one of the the essential part of the 21st-century organization, a basic engine driving improvement, intensity, including the development of organizations and the economy "Intellectual property can constitute more than 80 percent of a single company's value today.1 It's no surprise, then, that thievesarmed with means, motive, and opportunity-are in hot pursuit." [START_REF] Hammoudeh | Continuous Monitoring in Mission-Critical Applications Using the Internet of Things and Blockchain[END_REF][START_REF] Ghafir | DNS traffic analysis for malicious domains detection[END_REF] which briefly shows the fundamental of what the attacker wants to materialize.

B. Sabotaged critical organisation infrastructures

That means deletion of a database. As a result of this here is a brief explanation of the impact, according to the source I got from a website which says "a bigger threat than foreign states targeting political party emails was "covert sabotage, by chucking a cyber spanner, as it were, into the now cyber-dependent critical infrastructure that we have". This could include the financial system or energy utilities." [START_REF] Ghafir | DNS query failure and algorithmically generated domain-flux detection[END_REF][START_REF] Eyers | Sabotage Of Critical Infrastructure A Growing Risk, Says Former Chief Spy[END_REF].

C. Sites fully taken over

Gangs have taken steps to offer taken information to intruders; exploit taken information to assault casualties' colleagues; and plug casualties' "messy mysteries" based on reasonable research which can be viewed in websites especially local or government ones for instance "Some attackers took advantage of COVID-19 to coax people into opening malicious emails and attachments, while other ransomware groups agreed to an ad-hoc ceasefire on healthcare vendors." [START_REF] Novinson | The 11 Biggest Ransomware Attacks Of 2020 (So Far)[END_REF] In this scenario it could be that the attacker could be asking for ransom or else they could not release the site for the company until there is an agreement that needs to be done and that could cost the business a lot of money and time especially if the business is a marketing and financial industry.

D. Classified information compromised

Which includes employee of an organisation or a private user data a great example for this is when a person in the name of Harold T. Martin III, who was a worker at the N.S.A. cooperation personalized accessed the co-operation hacking unit, confessed his criminal activity, two years after his arrest in what may be described as the biggest breach of classified information in history. "F.B.I. agents who swarmed his modest home south of Baltimore in 2016 found stacks of documents and electronic storage devices stashed in his car, his home and even a garden shed." [START_REF] Ghafir | Blacklist-based Malicious IP Traffic Detection[END_REF].

V. DECTECTION AND DEFENCE

In the computing industry especially in the cyber security environment, there are requirements in how things should be before its being utilized. This case we must solve the APT and have a fully functioning plan to detect and defend users from APT. Which means it must have a multi-faceted approach on the role of network admins, individual users, and providers of security here are the activities during an APT session

A. Traffic monitoring

Now the best practice of preventing the instalment of stolen data extraction or something we call backdoors which is a malware that is known for negating normal authentication procedure for accessing a system. We are encouraged to use monitoring egress and ingress traffic, this because it does not only stop malware and stolen data extraction, inspecting the traffic inside the desired network can also aid in alerting security personnel to any strange behavior that may be a malicious activity.

There are examples of this types of process and one of this is a WAF which is a Web Application Firewall, and they are been installed into your network traffic that is inside the perimeter of your network filter traffic which are connected to your web application servers; therefore, it is protecting one of your weak attack sides. The other features that are included in a WAF is that it can clear out application layer attacks. This includes attack like SQL injections, Remote File inclusion also known as RFI as they are known for commonly used on a phase of an APT infiltration. Network firewalls which are also known to be a service for network monitoring are different in this matter. This is because they are able to provide a rough view of how users are communicating within the network while aiding in finding an internal traffic abnormality (strange or unusual large data transfer, abnormal logins). The end could indicate an APT attack that is being performed. It is possible to monitor system honeypots or file sharing access

In conclusion the use of incoming traffic monitoring system could be helpful in terms of detection and the removal of malwares like backdoor shells. This can be discovered just by the interception of remote request from the operators or users [START_REF] Raza | An IoT and Business Processes Based Approach for the Monitoring and Control of High Value-Added Manufacturing Processes[END_REF].

B. Access control

As for criminals, the largest and the most vulnerable soft spot in your security perimeter is your employees especially the ones in organizations'. As common as it may be, this is potentially why intruders can view network users since it can be an easy entrance to penetrate your defence, as they are hold withing the organisations' security perimeter.

Here below are the likely mistakes or deliberates that these target or employee do that could open a gateway for intruders and they are:  The carelessness of a user, especially those who fail to comply to the network security policies. This could grant access to unauthorized users and could lead to potential threats or attacks  Users who negotiate with the intruder to have the network access by them since the rouge employee is able to access them, leaving them open for planned threat by him/her and the intruder  Malicious insiders who deliberately took the advantage of their identity just to give the intruder access

The way to mitigate this type of problem is to develop an effective control that will require a thorough review of everyone in the organisations' especially the information that have been accessed. A great example is classifying data on a need-to-know basis, this known to be affective since its aids in blocking the intruder's ability to be able to infiltrate login credentials all from a low-level employee member, using it to gain access to classified materials. The other things that are required in defence are to secure key network access point with a two-factor authentication also known as 2FA. In this process it requires users to use a 2 nd form of authentication when they want to see the organisations' secretive materials, this usually done by sending a verification code to the users communicating devices, mobile phones, tablet and many more [START_REF] Ghafir | Tor-based malware and Tor connection detection[END_REF][START_REF] Ghafir | A Survey on Intrusion Detection and Prevention Systems[END_REF]. This proven to being able to prevent intruders masking as an authentic user to being able to explore the organisations' network.

C. Domain whitelisting and Application

The method of whitelisting consists of controlling the domains that can be accessed from the organisation or your network, this includes applications which can be installed by users. This is an effective method as is well known for being able to reduce the accomplishment rate of APT attack all by reducing the attack surfaces that are available. It is known that this type of security measure is not that deceiving, on the other hand it has been know that domains that are trusted can be compromised. Majority of security experts knows that malicious or strange files often arrive as a software that seems legitimate. To add into this, software that have old product version are vulnerable to compromisation and exploitation.

In order to have an efficient whitelisting, it is recommended to have a strict update policy. It must be enforced to reassure that your users or the organisations' employees are always running with the latest version of any application provided on the list Target areas: APT41 has straightforwardly focused on associations in at any rate 14 nations going back to as right on time as 2012. The gathering's surveillance crusades have focused on medical care, telecoms, and the cutting-edge area, and have verifiably included taking protected innovation. Their digital wrongdoing interruptions are generally evident among computer game industry focusing on, including the control of virtual monetary forms, and endeavored sending of ransomware. APT41 operations against advanced education, travel administrations, and news/media firms give some sign that the gathering additionally tracks people and directs observation. "Overview: APT41 is a prolific cyber threat group that carries out Chinese state-sponsored espionage activity in addition to financially motivated activity potentially outside of state control.

D. APT diagram

Related malware: APT41 has been noticed utilizing in any event 46 diverse code families and apparatuses.

Assault vectors: APT41 frequently depends on lance phishing messages with connections, for example, aggregated HTML (.chm) records to at first trade off their casualties. Once in a casualty association, APT41 can use more complex TTPs and send extra malware [START_REF] Novinson | The 11 Biggest Ransomware Attacks Of 2020 (So Far)[END_REF]. For instance, in a mission running close to 12 months, "APT41 compromised hundreds of systems and used close to 150 unique pieces of malware including backdoors, credential stealers, keyloggers, and rootkits. APT41 has also deployed rootkits and Master Boot Record (MBR) bootkits on a limited basis to hide their malware and maintain persistence on select victim systems [START_REF] Shane | Contractor Arrested In Biggest Breach Of U.S. Secrets Pleads Guilty[END_REF]." Suspected attribution: China Target areas: Numerous, including the government along with the transnational financial company, aerospace, and defense organizations, including high tech, construction along with the manufacturing, communications, media, and the insurance. "APT30 Suspected attribution: China Target sectors: Members of the Association of Southeast Asian Nations (ASEAN) Overview: APT30 is noted not just for sustained activity over a long period of time but also for being able to successfully modify and adjusting the source code to sustain the same tools, methods, and structure ever since at least the year of 2005. Indication reveals that the group arranges targets, possibly works in shifts in a collaborative environment and builds malware from a coherent development plan. The group has had the capacity to infect air-gapped networks and this has been happening according to the report, since 2005.

All of the associated malware includes: SPACESHIP, SHIPSHAPE, FLASHFLOOD "Attack vectors: APT30 uses a suite of tools that includes downloaders, backdoors, a central controller, and several components designed to infect removable drives and cross air-gapped networks to steal data. APT30 frequently registers its own DNS domains for malware CnC activities [START_REF] Ghafir | Disguised Executable Files in Spear-Phishing Emails: Detecting the Point of Entry in Advanced Persistent Threat[END_REF]." CONCLUSION Throughout the research and the actions that have been taken to stop intruders, in the security industry. This defense strategy which includes the mentioned techniques clearly showed crucial importance on how companies always come up with plans to reduce the risks of problems that could take place on their surface. Furthermore, the strict implementation of the company including education employees or users mostly focuses on the psychology of humans. The ability to distribute a workforce on requirements needed to not fall into the trap of social engineering. Those strategies mentioned are not technical and they deal with APT attacks taking the advantage of blind trusted individuals who have emails Personal or private accounts and attachments.

In the future, as companies or organizations are obtaining ideas and more defense systems to stop attacks like the APT, they should have an advanced team that is solely based on the APT attacks since it is a traditional attack, and it could be in a system unknown and undetected. On the other hand, some employees cannot be trusted in making sure they can mitigate attacks like APT since they could work with the intruder. Organizations should have a way of finding these rough employees as it will reduce the risk of APT forming or even worse changing the industry reputations. Local users should always make sure that if the software they are using should be updated luckily most software will remind the users when it needs updating. Furthermore, there have been many cases of APT and there will be more cases rising and it is based on preparations that can be used to detections and defense to reduce the incident of APT.

Fig. 1 .

 1 Fig. 1. Example of an APT attack (figure caption) REAL LIFE CASES OF APT In some areas of APT there are ways of measuring the attack and according to a research I found a site called fire eye, it clearly goes through a many case. here are the examples


Muhao Chen
email: muhaochen@tamu.edu

Jun Chen
email: jchen@qnlm.ac

Manoranjan Majji
email: mmajji@tamu.edu

Eminent Professor, Robert E Skelton
email: bobskelton@tamu.edu

Tensegrity system dynamics in fluids

Keywords: Tensegrity systems, Multibody dynamics, Nonlinear dynamics, Flexible structures, Fluid structure interaction

come   L'archive ouverte pluridisciplinaire

Introduction

Tensegrity, named by Buckminister Fuller [START_REF] Fuller | Synergetics: explorations in the geometry of thinking (Estate of R. Buckminster Fuller[END_REF] for the art form created by Ioganson (1921) and Snelson (1948) [START_REF] Lalvani | [END_REF]. It is a stable network of compressive members (bars) and tensile members (strings). If no compressive members are touching each other in the network, it is called a "class-1" structure. If there are a maximum of k compressive members connecting at one node, it is called a "class-k" structure [START_REF] Skelton | Tensegrity systems[END_REF]. From the definition, we know that the most fundamental property of the tensegrity structure is that all the structural members are axially loaded. The overall 3D structure is made from one-dimensional structure members. Thus, the structure members can be placed along the load path to make the best of each element [START_REF] Yuan | [END_REF]5], and the total mass of the structure can be greatly reduced [6]. And, since there is no material bending, the uncertainty of the onedimensional material movement also brings more accurate models [START_REF] Ma | Design and analysis of deployable clustered tensegrity cable domes[END_REF].

After years' study, many more properties and benefits of tensegrity systems have been discovered and verified by researchers: 1). Tensegrity in biological structures. For nano-scale evidence, Wang et al. found that living cell structure shape control is consistent with a tensegrity model [START_REF] Wang | [END_REF]. Seeman showed that tensegrity structures from DNA could be used to make nanomaterials [9]. For micro/macro-scale findings, Simmons et al. presented that the molecular structure of the spider fiber has both compressive and tensile members [10].

Scarr studied the significant features of a human's elbow as a tensegrity structure [11]. 2). Lightweight structure designs. For example, Fraternali et al. designed lightweight masonry domes and vaults by a tensegrity approach [12]. Ma et al. designed a lightweight tensegrity cantilever structure [START_REF] Ma | Composite Structures p[END_REF]. Wang et al. presented minimal mass algorithms for the design of active tensegrity structures [START_REF] Wang | [END_REF]. 3). Soft robot designs. For example, Caluwaerts et al. developed design and analysis tools for tensegrity robots and verified them by a spherical six-bar tensegrity ball [15]. Chen et al. derived a nonlinear shape control law for any class-k structures [16]. Kim et al. designed cable-driven spherical tensegrity robots for rolling locomotion [START_REF] Kim | Soft robotics[END_REF]. Sabelhaus et al. presented model-predictive control for tensegrity spine robots [START_REF] Sabelhaus | [END_REF]. Zheng et al. evaluated the robustness of the six-strut robot [19]. Baines et al. implemented the soft membrane to drive the tensegrity robot [20]. 4). Energy absorber/harvester. For example, Pajunen et al. analyzed the impact response of a 3D-printable tensegrity lander [21]. Baines et al. used pneumatic membrane actuators and a generalized path planning algorithm for rolling tensegrity robots [20].

There are a few studies have been conducted on tensegrity dynamics. For example, Sulton et al. derived linearized equations of motion for tensegrity models around arbitrary equilibrium configurations [22]. In 2005, Skelton [START_REF] Skelton | IUTAM symposium on vibration control of nonlinear mechanisms and structures[END_REF] presented one of the simplest dynamics forms for class-1 structure by using non-minimal coordinates and assumed the compressive elements to have no inertia about the longitudinal axis. Later, Wroldsen added wave forces to the dynamics of class 1 tensegrity structures [START_REF] Wroldsen | Proceedings of SPIE, the International Society for Optical Engineering[END_REF]. Cheong et al. extended the non-minimal dynamics to class-k tensegrity systems in a second-order matrix form [START_REF] Cheong | [END_REF]. Goyal et al. presented a compact matrix form of tensegrity dynamics by including massive strings, an analytical solution of Lagrange multipliers for class-k tensegrity structures, and bar length correction algorithms [26], a corresponding general software for modeling of any tensegrity structures can be found in [27]. Ma et al. formulated the tensegrity dynamics based on the finite element method by allowing the structure members to have elastic or plastic deformation [28]. Yang and Sultan investigated the deployment of foldable tensegrity-membrane systems [29]. Kan et al. derived tensegrity dynamics with clustered cables [30]. Multiphysics (thermoelasticity, fluidstructure interaction, piezoelectricity, magnetoelasticity, biomechanics, and intelligent systems) has become a popular topic in the past fifteen years [31]. And there are many research on fluid-structure interaction (FSI) problems, i.e., dynamics modeling for parafoil dynamics [32], study on the vibrations of cylinders [33], plate flut-ter and cavity [34], and beam/shell elements of a flapping wing [35]. However, in all these dynamics models, none of them considered fluid forces for tensegrity systems, which limited our ability to analyze fluid-based tensegrity structures. Applications include aquaculture, robotic fish, airplanes, etc. Fluid-structure interaction (FSI) simulation in the computational fluid dynamics (CFD) field is an option to do the analysis, but on the one hand, it requires skillful modeling experience in fluid dynamics and a lot of computational resources. On the other hand, the analytical fluid force model, i.e., the fluid force is proportional to the square of the relative velocity between the tensegrity element and fluid, has been successfully applied in the studies of aquatic animal swimming [36,37]. The accuracy of the analytical fluid force model was validated by the experimental data [38]. Thus, our study here provides a dynamics model to investigate the behaviors of tensegrity structure in the presence of fluid, also suitable for control design.

This paper is structured as follows: Section 2 derives the translation and rotation equation of motion of one single rod. Section 3 models the fluid particles interacting with the rod and the corresponding fluid force expressions on the rod. Section 4 stacks all the equation of motion for each rod to form the tensegrity network and develops a compact vector form of class-1 and class-k structures. Section 5 presents a compact matrix form of class-1 and class-k structures. Section 6 discusses the results of the dynamics by a 3D prism example. Section 7 gives the conclusions.

2 Dynamics of a single rod

Rotational dynamics

Let us consider a single rod immersed in an infinite volume of fluid field. The length of the rod is ||b|| = l, where b is a vector along the rod, shown in Fig. 1. Vector r gives the position of the rod center of mass O. We separate the forces on the rod into two parts: fluid force (f 0 and τ 0 ) and other external forces (f 1 and f 2 applied at the two nodes of the rod).

To describe the position of a rigid body, vectors are usually coordinatized in two different reference frames: inertia frame and body frame. Let the vectrices E = e 1 e 2 e 3 and B = b 1 b 2 b 3 denote the dextral set of unit vectors e i which are inertially fixed and the body-fixed dextral set of b i fixed in the coordinates of the rigid body [26]. The angular velocities of frame B relative to frame E is ω = Eω E = Bω B , where

ω B = ω B 1 ω B 2 ω B 3
T . Coordinate transformation by a unitary direction cosine matrix Θ satisfies:

B = EΘ, Θ T Θ = I. (1)
Take the first time derivative of Θ T Θ = I, one can get:

Θ T Θ = ωB =   0 -ω B 3 ω B 2 ω B 3 0 -ω B 1 -ω B 2 ω B 1 0   , (2)
which is skew-symmetric matrix. Then, the time derivative of B can be obtained,

Ḃ = E Θ = BΘ T Θ = B ωB . (3)
ḃ = Ḃb B + B ḃB = Ḃb B = B ωB b B . (4)
Then, since we have ωb = -bω and bb = -l 2 I +bb T , one can obtain:

b × ḃ = (Bb B ) × (B ωB b B ) = B bB ωB b B (5) = -B bB bB ω B = -B( bB ) 2 ω B (6) = B(b B T b B I -b B b B T )ω B (7) = B( b B 2 I -b B b B T )ω B (8) = B   l 2 I -   0 0 0 0 0 0 0 0 l 2     ω B (9) = B I 2 0 0 0 l 2 ω B = Bl 2   ω B 1 ω B 2 0   = l 2 ω b . ( 10
)
The angular momentum of bar b about its mass center O is:

h = Jb × ḃ = m b 12 b × ḃ, ( 11
)
where J is moment of inertia of the rod J = 1 12 m b l 2 , m b is bar mass. Take the first time derivative of Eq. ( 11), we obtain:

ḣ = J ḃ × ḃ + Jb × b = Jb × b. ( 12
)
The time derivative of the angular momentum vector h equals to the sum of torques τ = τ e +τ 0 acting on the bar member about its center of mass, where τ 0 is the torque generated by the fluid force,

τ e = 1 2 b × (f 2 -f 1
) is torque given by the force at the two nodes of the bar. Then we have the following:

ḣ = Jb × b = τ , (13)
τ = τ e + τ 0 = 1 2 b × (f 2 -f 1 ) + τ 0 . ( 14
)
Eq. ( 13) can be written in any coordinates, but we choose inertial coordinates for simpler forms of final equations. To simplify notation hereafter, we define b = b E where b = Bb B = Eb E [26]. Using inertial coordinates, Eq. ( 13) can be written as:

J bb = 1 2 b(f 2 -f 1 ) + τ 0 . (15)
Since bar length is a constant, we have ||b|| = l or b T b = l 2 . Take the first and second derivative of the bar length constraint, we get:

ḃT b + b T ḃ = 2b T ḃ = 0, ( 16
) ḃT ḃ + b T b = 0, b T b = -ḃT ḃ = -|| ḃ|| 2 . ( 17
)
Write Eq. ( 15) and Eq. ( 17) together in matrix form as:

b b T b = 1 2J b(f 2 -f 1 ) -ḃT ḃ + τ0 J 0 . (18)
= b b T † 1 2J b(f 2 -f 1 ) -ḃT ḃ + b b T † τ0 J 0 (19) = 1 l 2 -b b 1 2J b(f 2 -f 1 ) -ḃT ḃ + 1 l 2 -b b τ0 J 0 (20) = 1 l 2 (- 1 2J bb (f 2 -f 1 ) -b ḃT ḃ) - 1 l 2 b τ 0 J (21) = - 1 2Jl 2 (-l 2 I + bb T )(f 2 -f 1 ) - 1 l 2 b ḃT ḃ - 1 Jl 2 bτ 0 , (22)
where superscript " †"denotes a matrix pseudoinverse.

Rearranging equations gives the final vector form of the single rod rotational dynamics:

J b = 1 2 (f 2 -f 1 ) - 1 2l 2 bb T (f 2 -f 1 ) - J l 2 b ḃT ḃ - 1 l 2 bτ 0 (23) = 1 2 (I - bb T ||b|| 2 )(f 2 -f 1 ) - J l 2 b ḃT ḃ - 1 l 2 bτ 0 . (24)

Translational dynamics

The position of center of mass of the bar is given by vector r. The fluid force working at the center of mass of the bar is f 0 and the sum of the internal forces from the strings and the external forces, acting on the two ends of the bar, is described by f 1 and f 2 , where r = Er E ,

f i = Ef E i . From Newton's second law, we have r = r E and f i = f E i . m b r = f 1 + f 2 + f 0 ,
which can be written in the inertial coordinates as:

m b r = f 0 + f 1 + f 2 . ( 25
)
3 Modeling of fluid forces on rod

The fluid force contributions are only considered on bar members due to their relatively bigger diameter than the strings [START_REF] Wroldsen | Proceedings of SPIE, the International Society for Optical Engineering[END_REF]. The continuously distributed fluid forces along the bar can be represented by a single force f 0 at the geometrical center O and a torque τ 0 , which are functions of velocities and angular velocities of the bar. Both vicious and fluid inertia forces are considered.

The fluid force on the geometrical center of the bar f 0 (a sum of tangent force f t and normal force f n ) and torque (τ 0 ) are given by [START_REF] Chen | Mechanisms underlying undulatory swimming: From neuromuscular activation to body-fluid interactions[END_REF]:

f 0 = f t + f n , (26)
-f t = 2.7c t l ρµd||v n ||v t , (27)
-f n = sgn(v n )c p (ρv 2 n /2)ld + c a ρπ(d/2) 2 l(a n -2v t θ), (28)
-τ 0 = 1 12 sgn(v n )c p ρdl 3 v n θ + 1 12 c a ρπ(d/2) 2 l 3 θ, ( 29
)
where θ = bḃ

||b|| 2 , θ = bb ||b|| 2 , ( 30
)
v t = bb T ||b|| 2 ṙ, v n = (I - bb T ||b|| 2 ) ṙ, ( 31
)
a t = bb T ||b|| 2 r, a n = (I - bb T ||b|| 2 )r, ( 32
)
where ρ and µ are the fluid density and viscosity, d is the diameter of the rod, c t = 1 and c p = 0.9 ∼ 1.1 are the drag coefficients on a smooth circular cylinder in the range of Reynolds number 20 < R e < 10 5 from experiments [37], c a is the dimensionless coefficient to account for different object shape (for cylinder, c a = 1). v t , v n are velocities and a t , a n are accelerations of the center of mass of the bar in tangential and normal direction, θ and θ are the angular velocity and acceleration of the rod about e 3 axis (the axis along the bar direction), sgn(•) takes the sign of the variable in the parentheses.

In order to write the dynamic equations in a standard second order form, we separate fluid force and torque into two parts, acceleration dependent one (with subscript A) and the rest (with subscript R), and move the acceleration dependent parts to the left hand side of Eq. ( 24). Let us first define:

f 0 = f R + f A , τ 0 = τ R + τ A , (33)
-f R := sgn(v n )c p (ρv 2 n /2)ld -2m A v t θ + 2.7c t l ρµd||v n ||v t , (34)
-f A := m A a n , (35)
-τ R := sgn(v n ) 1 12 c p ρdl 3 v n θ, -τ A := 1 12 m A l 2 θ, (36)
where m A = c a ρπ(d/2) 2 l is the mass of fluid being accelerated by the bar, usually called the added mass of fluid. The moment of inertia of a rod is J = 1 12 m b ||b|| 2 , replace τ 0 in Eq. ( 24) by Eq. ( 33), we have:

m b 12 + m A 12 I - bb T ||b|| 2 b + m b 12 || ḃ|| 2 ||b|| 2 b = 1 2 I - bb T ||b|| 2 (f 2 -f 1 ) - b ||b|| 2 τ R , (37)
m b + m A I - bb T ||b|| 2 r = f 2 + f 1 + f R . ( 38
)
4 Vector form of tensegrity system dynamics with fluid forces incorporated

H i := I 3 - b i b T i ||b i || 2 , m i := m bi I 3 + m Ai H i , (39)
g i := m bi 12 || ḃi || 2 ||b i || 2 , κ i := - bi ||b i || 2 . (40)
Then, Eq. ( 37) and Eq. ( 38) can be written into a compact form:

m 1 12 m r b + 0 0 0 ĝ ⊗ I 3 r b = 2C r ⊗ I 3 1 2 Ĥ(C b ⊗ I 3 ) f + f R κτ R , ( 41
)
where • is a diagonal operator that converts a vector into a diagonal matrix,

f = f T 1 f T 2 • • • f T 2β T is the vec- tor of nodal forces, n = n T 1 n T 2 • • • n T 2β
T is the vector of node coordinates.

The nodal force f is then given by:

f = w + f s + f d + (P T c ⊗ I 3 )t c , ( 42
)
where w is the external forces applied at each node, f s is the string forces, f d is the string damping force, t c is the internal force between joints. If the tensegrity system is subject to m constraints, we have:

(P c ⊗ I 3 )n = d c , (43)
where P c ∈ R m×2β is the constraint matrix, and d c ∈ R 3m is a constant vector. Eq. ( 43) describes two types of constraint: P cij = 1, P c ik = -1 and d ci = 0 when nodes j and k are connected by a universal joint, where the subscript i is the i th constraint, and j, k are the nodes index; P ci,j = 1 and d ci ∈ R 3 when the j th node is fixed to the ground. The forces in the strings f s can be written as:

f s = -(C T s γC s ⊗ I 3 )n, ( 44
)
γ i := k i (1 -l 0i /l i ), l i > l 0i 0, else , (45)
where k i , l i and l 0i are the string stiffness, length and rest length. Since we have l 2 i = s T i s i ⇒ li = s T i ṡi /l i for each string, the damping force in each string is calculated by:

f di := -µ i li s i l i = -µ i s T i ṡi l 2 i s i , (46)
written in a vector form, we have:

f d = -(C T s ζC s ⊗ I 3 )n, ( 47
)
ζ i : =    µ i s T i ṡi l 2 i , l i > l 0i 0, else , (48)
where µ i (N•s/m) is the damping coefficient. The angular displacement between the contact surface, as described by the constraint Eq. ( 43), introduces resistive torques to bars. The resistive torque on bar i by bar k due to angular velocities ω i and ω k is given by:

τ d ik = -η ik (ω i -ω k ), ω i = bi ḃi ||b i || 2 , ( 49
)
where η ik (N•m•s) is the damping coefficient between bar i and bar k. The vector form of resistive torque is given by:

τ d = -2C r P T c ηP c 2C T r ω = 4C r P T c ηP c C T r κḃ . ( 50
)
τ d is added to Eq. ( 41) by replacing τ R by τ R + τ d . Replace r, b, f and τ R in Eq. ( 41) by r = (C r ⊗I 3 )n, b = (C b ⊗ I 3 )n, Eq. ( 42) and τ R + τ d , where τ d is given by Eq. ( 50), and also left multiply C T r C T b ⊗ I 3 to Eq. ( 41) to make the mass matrix symmetric, we have:

M n + C ṅ + Kn = D(w + P T t c ) + z, (51)
where:

M := (C T r ⊗ I 3 ) m(C r ⊗ I 3 ) + 1 12 (C T b ⊗ I 3 ) m(C b ⊗ I 3 ), (52)
C := 4(C T b ⊗ I 3 )κ T (C r P T c ηP c C T r ⊗ I 3 )κ(C b ⊗ I 3 ), (53)
D := I - 1 2 (C T b ⊗ I 3 )∆(C b ⊗ I 3 ), (54)
∆ := block diagonal b i b T i ||b i || 2 , (55)
K := D(C T s (γ + ζ)C s ⊗ I 3 ) + (C T b ⊗ I 3 )ĝ(C b ⊗ I 3 ), (56)
z := (C T r ⊗ I 3 )f R + (C T b ⊗ I 3 )κτ R , (57)
P := P c ⊗ I 3 . ( 58
)
The vector z is how we add the fluid force, which is composed of a force vector f R at location r and a torque τ R both of which are acceleration independent, and m i is how we add the added mass effect of fluid (the fluid being accelerated by the bar) into the equation of motion. Note that Eq. ( 52) is the equation of motion for class-1 tensegrity with P = 0. Suppose the constraints are linearly independent, we can write:

(P c ⊗ I 3 )n = d c ⇒ n = P † d c + P ⊥ n c , (59)
P † = P † c ⊗ I 3 , P ⊥ = P ⊥ c ⊗ I 3 , (60)
where † means moore-penrose inverse, P ⊥ is the right null space of P . Replace n in Eq. ( 51) by above expression, we get:

M P ⊥ -DP T nc t c = -CP ⊥ ṅc -K(P † d c + P ⊥ n c ) + Dw + z. ( 61
)
Rewrite the above equation in a compact form, we have the equation of motion for class-k tensegrity:

nc t c = M P ⊥ T -(DP T ) T Π -1 ϕ, (62)
where:

Π := M P ⊥ (M P ⊥ ) T + DP T (DP T ) T , (63)
ϕ := -CP ⊥ ṅc -K(P † d c + P ⊥ n c ) + Dw + z. (64)
4.2 Normalization on P ⊥ to reduce numerical error

The condition number of Π is in order of 10 6 when the unit mass is in order of 10 -1 [START_REF] Chen | Mechanisms underlying undulatory swimming: From neuromuscular activation to body-fluid interactions[END_REF]. If 99% of unit mass is on bar 1, the condition number of Π will be in order of 10 10 . Note that the term DP T in Π has entries close to 1 when entries of P is so (the case of nodal constraint). To reduce numerical error, P ⊥ term in Π can be normalized as below:

M (P ⊥ /δ) -DP T nc t c = -C(P ⊥ /δ) ṅc - K(P † d c + (P ⊥ /δ)n c ) + Dw + z, nc := n c δ, ( 65
)
where δ is a scalar such that the entries of M/δ is close to 1. After such normalization, the condition number of Π can be reduced from 10 10 to 10 8 .

Computational effort and numerical error reduction

From the structure of the matrices, we know that there are many repeated block entries in Eq. ( 62) resulted from the structure of C r := 1 2 I I and C b := -I I . We expand Eq. (62) to identify those repeated blocks to allow the computer to compute them only once. This also reduces the order of matrices by half in arithmetic operations. Furthermore, most matrices in Eq. (62) are block diagonal, and the matrices P , P † , and P ⊥ are sparse, the functions for sparse matrix in matlab help to reduce the computation effort and possible numerical errors. For example, the command sparse creates a sparse matrix by saving only the nonzero entries and their indices, and the arithmetic operators operate, for example, when we multiply two matrices, only on nonzero entries, which avoid unnecessary computation and numerical error from those zero entries.

By expanding, matrices in Eq. ( 62) can be rewritten as:

M = 1 3 m 1 6 m 1 6 m 1 3 m , C = C 11 -C 11 -C 11 C 11 , (66)
C 11 := κT ((P T 1 + P T 2 )η(P 1 + P 2 ) ⊗ I 3 ) κ, (67)

P =: P 1 P 2 , D = I - 1 2 ∆ -∆ -∆ ∆ , (68)
K = K 11 K 12 K T 12 K 22 ⊗ I 3 + ĝ -ĝ -ĝ ĝ ⊗ I 3 - 1 2 ∆ ((K 11 -K T 12 ) ⊗ I 3 ) ∆ ((K 12 -K T 22 ) ⊗ I 3 ) , ( 69
)
K 11 := C T s1 (γ + ζ)C s1 , K 12 := C T s1 (γ + ζ)C s2 , (70)
K 22 := C T s2 (γ + ζ)C s2 , C s =: C s1 C s2 , (71)
where means the negative of the first row of the 2 by 1 block matrix. With above definition, the terms in ϕ in Eq. ( 62) were written as:

CP ⊥ ṅc = C 11 (P ⊥1 -P ⊥2 ) ṅc , P ⊥ =:

P ⊥1 P ⊥2 , (72)
P † d c + P ⊥ n c = p c1 + P ⊥1 n c p c2 + P ⊥2 n c , P † d c =: p c1 p c2 , (73)
Dw = w 1 w 2 - 1 2 ∆(w 1 -w 2 ) , w =: w 1 w 2 , (74)
The terms in Π matrix in Eq. ( 62) are written as:

M P ⊥ = m( 1 3 P ⊥1 + 1 6 P ⊥2 ) m( 1 6 P ⊥1 + 1 3 P ⊥2 ) , ( 75
)
DP T = P T 1 P T 2 - 1 2
∆(P T 1 -P T 2 ) , (76)
Π := M P ⊥ (P ⊥ ) T M + DP T P D =:

Π 11 Π 12 Π T 12 Π 22 . ( 77
)
Using block matrix inverse formula:

A B C D -1 = Θ 11 Θ 12 Θ 21 Θ 22 , (78)
where:

Θ 11 = (A -BD -1 C) -1 , (79)
Θ 12 = -A -1 B(D -CA -1 B) -1 , (80)
Θ 21 = -D -1 C(A -BD -1 C) -1 , (81)
Θ 22 = (D -CA -1 B) -1 . ( 82
)
The inverse of Π is then given by:

Π -1 = Λ 11 Λ 12 Λ 21 Λ 22 , (83)
where:

Λ 11 = (Π 11 -Π 12 Π -1 22 Π T 12 ) -1 , (84)
Λ 12 = -Π -1 11 Π 12 (Π 22 -Π T 12 Π -1 11 Π 12 ) -1 , (85)
Λ 21 = -Π -1 22 Π T 12 (Π 11 -Π 12 Π -1 22 Π T 12 ) -1 , (86)
Λ 22 = (Π 22 -Π T 12 Π -1 11 Π 12 ) -1 . ( 87
)
For class-1 tensegrity, we still need to inverse M directly, which is given by:

M -1 = 4 m-1 -2 m-1 -2 m-1 4 m-1 . ( 88
)

Calculate initial values of n c from n

In simulation, we set the initial conditions of the structure, and then check the time response. n is chosen to satisfy the constraints, the initial values of n c for ode solvers (ode45 is usually used) can be calculated from:

P ⊥ n c = (n -P † d c ), (89)
we have:

n c = (P ⊥ ) † (n -P † d c ). ( 90
)

Matrix formulation of tensegrity dynamics with fluid force incorporated

The dynamics of any given rod member in a tensegrity structure in the presence of fluid force are given by Eq. ( 37) and Eq. (38). To describe a full tensegrity structure, it basically stacks 2β vector equations for a system containing β rod members. In order to simplify the equation structure, one can assemble them in a matrix form.

We name nodes at the base of bar vectors as

N := n 1 • • • n 2β , for a network of β bars, define the 3 × β matrix B := b 1 • • • b β .
The mass centers of the bars is defined as R := r 1 • • • r β , where r i is the ith column of matrix R. The internal forces acting on nodes caused by string tensions is N C T s γC s , the full force matrix expression can then be written as F := W -N C T s γC s . Fluid forces can be assembled as

F R := f n1 1 + f t1 • • • f n1 β + f t β , ∆ R := τ R1 • • • τ R β . Let us also define B := b1 • • • bβ and e := 1 • • • 1 T ,
we have the following expression:

B m 1 12 + B ḂT Ḃ ml -2 1 12 + B T F C T b l-2 1 2 = F C T b 1 2 + e T ⊗ I 3 BT ∆ R l-2 , (91)
R m = 2F C T r + F R , (92)
where F R , τR and m are how we add the accelerationindependent and dependent fluid forces and torques into the equation, and l i := ||b i ||. Define [26]:

-λ := ḂT Ḃ ml -2 1 12 + B T F C T b l-2 1 2 , (93)
M := 1 12 C T b mC b + C T r mC r , (94)
K := C T s γC s -C T b λC b , (95)
with the following formula:

1 2 C T b C b + 2C T r C r = I, C b C r = 1 2 C T b 2C T r -1 , (96)
Eq. ( 91) and Eq. ( 92) can be written into:

N M + N K = W + e T ⊗ I 3 BT ∆ R l-2 C b + F R C r . ( 97
)
When the system is subject to constraint of form [26]:

N P = D, P ∈ R β×m . (98)
The constraint force T is given by:

T = ΩP T , (99)
and W in Eq. ( 97) is replaced by W +ΩP T . A analytical solution for Ω can be found in [26].

We can also add an additional constraint that the e 2 coordinate of n 1 remains constant d y , i.e.: e T 2 N e 1 = d y , e 2 = 0 1 0

T ∈ R 3 , (100)
e 1 = 1 0 T ∈ R 3β . (101)
The total constraint is N P e T 2 N e 1 = D d y which can not be written in the form of Eq. (98). So Eq. ( 98) is only applicable to constraints on vectors n i , i = 1, • • • , 2β, but not on individual entries n ix , n iy , n iz . The vector form equation of motion in section 4 is more suitable for compact expression:

e T 2 (e T 1 ⊗ I 3 )n = d y , Pn = d, (102)
where:

P := P ((e 1 ⊗ I 3 )e 2 )
T , d :=

d d y . (103)

Numerical study

It has been observed that DNA has a similar structure as tensegrity prism, where rigid bundles of DNA double helices resist compressive forces exerted by segments of single-stranded DNA that act as tension-bearing cables [START_REF] Liedl | [END_REF]. Many research on tensegrity prism has been conducted to study its properties in lightweight [START_REF] Skelton | Tensegrity systems[END_REF], configuration method for connecting [41], deployability by stimulus-responsive polymers [42], stiffness of the prism towers [43], etc. Thus, 3D three-bar prism example, shown in Fig. and6. The coordinate and velocity time history of nodes 4, 5, and 6 of the prism unit without fluid are given in Figs. 7 and8. From these two figures, we can see that the prism is experiencing a free periodical oscillation. The coordinate and velocity time history of nodes 4, 5, and 6 of the prism unit in the presence of fluid are given in Figs. 9 and 10. It is shown that the structure quickly damped to its equilibrium. Comparing the results with no fluid results, we can see that the water provides large damping to the dynamics response of the structure, which agrees with the physics. The bar length and fixed node errors of the prism are shown in Figs. 3 and time history of nodes 4, 5, and 6 of the prism unit in the presence of fluid with 3 m/s velocity in the x-direction are given in Figs. 12 and13. The result showed that the structure quickly damped to its equilibrium, and the whole structure moves towards the x-direction, shown in Fig. 11. Comparing the results for fluid with and without velocity, we can conclude that the fluid provides large damping to the dynamics response of the structure and a pushing force to the whole structure in the fluid velocity direction, which matches well with the physics.

Example 3: Landing in the presence of fluid

Landing problem has attracted a lot attentions from various researchers [START_REF] Kim | Soft robotics[END_REF]44,[START_REF] Rimoli | 57th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference[END_REF]. All these simulations are assuming a non-fluid environment. This is true for lunar or asteroid landing, but for other scenarios, the influence of fluid should be considered. For example, landing on the ocean floor, Earth, Mars, and Titan would require the consideration of fluid. The example is to demonstrate the capability of the formulation to perform the dynamic simulation with inputs from the external environment (gravity and fluid). A dynamic simulation result was shown when the prism lander was dropped from a height of 1 m (distance between the bottom nodes of the prism and the ground). For this simulation, the frictionless ground was modeled as a spring-damper system of stiffness k g = 3,000 N/m and damping c g = 3 N•s/m. An initial prestress value of all the bars and vertical strings λ = λ v = 15 N/m, force densities in the top and bottom strings are

γ t = γ b = λ/ √ 3 = 5 √ 3 N/m,
which result in selfequilibrium for the structure without gravity force [START_REF] Skelton | Tensegrity systems[END_REF]. We should point out that adding gravity in landing to the structure will slightly change the equilibrium configuration, and one can observe this phenomenon before the structure hitting the ground, shown in the x and y coordinates and velocities at 0 s ∼ 0.52 s in Fig. 16 and 0 s ∼ 0.84 s in Fig. 18.

For this case, we compare the motion of the flexible structure landing simulation with and without fluid. can see that it takes about 0.52 s for the prism to hit the ground. After hitting the ground, the prism is experiencing a big bouncing up and down with small damping from the ground. The coordinate and velocity time history of nodes 4, 5, and 6 of the prism unit in the presence of fluid are given in Figs. 18 and19. It is shown that it takes a little longer, 0.84 s, for the prism to hit the ground because of the fluid damping influence. After hitting the ground, the structure quickly damped to its equilibrium. Comparing the results with no fluid results, we can see that the water provides large damping to the landing process, which takes a longer time to land on the floor but damped quicker to its equilibrium.

Conclusion

This paper first derived the dynamics of one single rod. The continuously distributed fluid forces along the bar are represented by a single force at the geometrical center and a torque. Then, the fluid forces are added to the derived dynamics of equations. By stacking all the equations of motion for each rod, compact vector and matrix forms of class-1 and class-k tensegrity dynamics are formulated. Finally, based on a three-dimensional prism model, we simulate and compared results without considering fluid forces, fluid velocity, and gravity forces. Results show that fluid forces and velocity can give the structure damping and pushing force to the dynamics response of the structure, which agrees with the physics. This study gives an analytical dynamics formulation of fluid-structure interaction of any class-1 and class-k tensegrity structures. The integrated dynamic model of fluid-structure interaction developed in this paper is suitable for the design of feedback control or gives a quick and qualitative estimation of the behaviors of tensegrity systems in fluid in the initial design.

Fig. 1

 1 Fig. 1 Tensegrity bar member vector with fluid force nomenclature.

  Since the bar vector b described in body frame B is b = Bb B , where b B represents the components of the vector b as viewed in coordinate frame B. With b 3 along the rod, we have b B = 0 0 l T and ḃB = 0. The time derivative of vector b is:

  One can solve this linear algebra problem for b. Notice that b b T has full column rank, the unique solution for b is: b

4. 1

 1 Class-1 tensegrity with fluid incorporated Let the connectivity matrices of class-1 structure be C b = -I β I β for bars, where β is the number of bars in the system, and C s for strings, and C r = 1 2 I β I β [3]. Let us define:

Fig. 2 A

 2 Fig. 2 A 3D three-bar prism, thick black lines are bars and thin red lines are strings, prism height is 0.5 m, and bar length is 0.6952 m.

2 ,Fig. 3

 23 Fig. 3 Bar length errors of the prism.

4 . 6 . 2

 462 Example 2: Free oscillation in the fluid with an inlet velocity For this case, we compare the motion of the flexible structure dynamic with a fluid inlet velocity of 0 m/s and 3 m/s in the x-direction. The three-figure timelapse plots are given in Figs. 6 and 11. The coordinate and velocity time history of nodes 4, 5, and 6 of the prism unit in the presence of static fluid are given in Figs. 9 and 10. It is shown that the structure quickly damped to its equilibrium. The coordinate and velocity

Fig. 4

 4 Fig.4x, y, and z coordinates error of the fixed node.

Fig. 5

 5 Fig. 5 Free oscillation of the prism (no fluid) at T = 0s, 1s, and 2s.

Fig. 6

 6 Fig. 6 Free oscillation of the prism (immersed in the fluid) at T = 0s, 1s, and 2s.

Fig. 7

 7 Fig. 7 Coordinate time history of node 4, 5 and 6 from the free oscillation of the prism (no fluid).

Fig. 8

 8 Fig. 8 Velocity time history of node 4, 5 and 6 from the free oscillation of the prism (no fluid).

Fig. 9

 9 Fig. 9 Coordinate time history of node 4, 5 and 6 from the free oscillation of the prism (immersed in the fluid).

Fig. 10

 10 Fig. 10 Velocity time history of node 4, 5 and 6 from the free oscillation of the prism unit (immersed in the fluid).

Fig. 11

 11 Fig. 11 Dynamics simulation of the prism (fluid inlet velocity 3 m/s in x-direction) at T = 0s, 1s, and 2s.

Fig. 12

 12 Fig. 12 Coordinate time history of node 4, 5 and 6 from the dynamics simulation of the prism unit (fluid inlet velocity 3 m/s in x-direction).

Fig. 13

 13 Fig. 13 Velocity time history of node 4, 5 and 6 from the dynamics simulation of the prism unit (fluid inlet velocity 3 m/s in x-direction).

  The three-figure time-lapse plots are given in Figs. 14 and 15. The coordinate and velocity time history of nodes 4, 5, and 6 of the prism unit without fluid are given in Figs. 16 and 17. From these two figures, we

Fig. 14

 14 Fig. 14 Landing simulation of the prism (no fluid) at T = 0s, 1s, and 2s.

Fig. 15

 15 Fig. 15 Landing simulation of the prism unit (immersed in the fluid) at T = 0s, 1s, and 2s.

Fig. 16

 16 Fig. 16 Coordinate time history of node 4, 5 and 6 from the landing simulation of the prism unit (no fluid).

Fig. 17

 17 Fig. 17 Velocity time history of node 4, 5 and 6 from the landing simulation of the prism unit (no fluid).

Fig. 18

 18 Fig. 18 Coordinate time history of node 4, 5 and 6 from the landing simulation of the prism unit (with fluid).

Fig. 19

 19 Fig. 19 Velocity time history of node 4, 5 and 6 from the landing simulation of the prism unit (with fluid).