Application of discrete shape function in absolute nodal coordinate formulation

Zhicheng Song; Jinbao Chen; Chuanzhi Chen; Zhicheng Song; Jinbao Chen; Chuanzhi Chen

doi:10.3934/mbe.2021234

Mathematical Biosciences and Engineering

2021, Volume 18, Issue 4: 4603-4627. doi: 10.3934/mbe.2021234

Previous Article Next Article

Research article Special Issues

Application of discrete shape function in absolute nodal coordinate formulation

Key Laboratory of Exploration Mechanism of the Deep Space Planet Surface, Ministry of Industry and Information Technology, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China

Received: 09 April 2021 Accepted: 25 May 2021 Published: 27 May 2021

To solve integrals in the absolute nodal coordinate method and address the difficulty in applying it to an arbitrary-section beam, this paper focuses on two methods involving single integrals:the invariant matrix method and the Gerstmayr method, with cross-section characteristics by applying the interpolation of a discrete function. Such single integrals demonstrate that the nodal coordinate method can be applied to an arbitrary-section beam. The Euler–Bernoulli beam used in engineering structures is characterised by a symmetrical cross-section, small section size, zero odd integrals and negligible high-order even integrals, which simplifies the single integrals of the two methods. Finally, the Gaussian integration is adopted to improve the solving efficiency of elastic force and force Jacobian.

Keywords:

absolute nodal coordinate formulation (ANCF),
tensor of stress,
invariant matrices,
discrete shape function

Citation: Zhicheng Song, Jinbao Chen, Chuanzhi Chen. Application of discrete shape function in absolute nodal coordinate formulation[J]. Mathematical Biosciences and Engineering, 2021, 18(4): 4603-4627. doi: 10.3934/mbe.2021234

Related Papers:

[1]	Qinghua Liu, Jing Wang, Keying Song . Displacement identification of oil storage tank and calibration of tank capacity table model based on indicator function integral. Mathematical Biosciences and Engineering, 2024, 21(2): 1872-1883. doi: 10.3934/mbe.2024082
[2]	Ping Zhou, Yuting Zhang, Yunlei Yu, Weijia Cai, Guangquan Zhou . 3D shape measurement based on structured light field imaging. Mathematical Biosciences and Engineering, 2020, 17(1): 654-668. doi: 10.3934/mbe.2020034
[3]	Alessia Andò, Simone De Reggi, Davide Liessi, Francesca Scarabel . A pseudospectral method for investigating the stability of linear population models with two physiological structures. Mathematical Biosciences and Engineering, 2023, 20(3): 4493-4515. doi: 10.3934/mbe.2023208
[4]	Mohamed S. Eliwa, Buthaynah T. Alhumaidan, Raghad N. Alqefari . A discrete mixed distribution: Statistical and reliability properties with applications to model COVID-19 data in various countries. Mathematical Biosciences and Engineering, 2023, 20(5): 7859-7881. doi: 10.3934/mbe.2023340
[5]	Pratik Purohit, Prasun K. Roy . Interaction between spatial perception and temporal perception enables preservation of cause-effect relationship: Visual psychophysics and neuronal dynamics. Mathematical Biosciences and Engineering, 2023, 20(5): 9101-9134. doi: 10.3934/mbe.2023400
[6]	Jin Li . Linear barycentric rational collocation method to solve plane elasticity problems. Mathematical Biosciences and Engineering, 2023, 20(5): 8337-8357. doi: 10.3934/mbe.2023365
[7]	Kelum Gajamannage, Erik M. Bollt . Detecting phase transitions in collective behavior using manifold's curvature. Mathematical Biosciences and Engineering, 2017, 14(2): 437-453. doi: 10.3934/mbe.2017027
[8]	Maysaa Al Qurashi, Saima Rashid, Sobia Sultana, Hijaz Ahmad, Khaled A. Gepreel . New formulation for discrete dynamical type inequalities via $ h $-discrete fractional operator pertaining to nonsingular kernel. Mathematical Biosciences and Engineering, 2021, 18(2): 1794-1812. doi: 10.3934/mbe.2021093
[9]	Muhammad Bilal Khan, Pshtiwan Othman Mohammed, Muhammad Aslam Noor, Khadijah M. Abualnaja . Fuzzy integral inequalities on coordinates of convex fuzzy interval-valued functions. Mathematical Biosciences and Engineering, 2021, 18(5): 6552-6580. doi: 10.3934/mbe.2021325
[10]	Ke Bi, Yue Tan, Ke Cheng, Qingfang Chen, Yuanquan Wang . Sequential shape similarity for active contour based left ventricle segmentation in cardiac cine MR image. Mathematical Biosciences and Engineering, 2022, 19(2): 1591-1608. doi: 10.3934/mbe.2022074

Abstract

1. Introduction

Large size, high precision, lightweight and multi-shape have become the inevitable demand of future spacecraft, and flexible multi-body systems rich in large deformation components have been widely used in the aerospace field ^[1]. For example, spacecraft solar panels ^[2], satellite expandable antennas ^[3], space telescopes, etc. ^[4]. In order to obtain the mechanical mechanism of such mechanisms with large deformation and large rotation ^[5,6,7], much work has been done in this field.

The dynamics of the flexible multi-body system focused on the dynamic behavior of a system composed of a deformable object and rigid body when it experiences a wide range of spatial motions ^[8]. Turcic et al. ^[9] proposed Kineto Elastio Dynamic Analysis (KED) method firstly to deal with dynamics of flexible multi-body systems. This method decomposes flexible multi-body dynamics into multi-rigid body dynamics and structural dynamics provides high analytical accuracy for multi-body systems with low flexibility due to low-speed movement ^[10]. With the emergence of lightweight high-speed mechanism, the dynamic equation constructed by the KED method could not explain the coupling effect on the elastic deformation and large deformation caused by high-speed movements, so the ﬂoating frame of reference formulation (FFRF) has been proposed ^[11,12]. This method constructs the floating coordinate system on the flexible body and establishes the dynamic model by using the rigid body coordinates of a floating coordinate system and the node coordinates (or modal coordinates) with the flexible body. The geometrically exact formulation (GEF) can reasonably describe finite rotation, but it has higher requirements for interpolation strategy and numerical integration method. Shabana ^[13] proposed an absolute nodal coordinate formulation (ANCF) based on the finite element and continuum mechanics theory. In this method, the generalized coordinates of element nodes were defined in the global coordinate system, and slope vectors were used to replace the angular coordinates in the traditional finite element method. The dynamic equation has the characteristics of a constant inertia matrix, the absence of Coriolis force, and centrifugal force ^[14,15]. Compared with the traditional method of FFRF ^[16], the GEF ^[17], ANCF can describe the flexible multi-body system more accurately and plays an important role in dealing with the dynamics of flexible multi-body systems ^[18,19].

The ANCF is widely used in ﬂexible multibody simulations. Yakoub ^[20] applied the three-dimensional beam element to simulate the wheel-rail contact dynamics problem, which applied the ANCF to practical engineering firstly. Li et al. ^[21] simulated the antenna reﬂector mesh based on the slope deﬁcient ANCF beam element and applied it to the form-ﬁnding analysis of the reﬂector with the antenna. Tian et al. ^[22] introduced the ANCF in several types of deployable space structure and proposed a parallel computational framework to improve the computational efficiency. Lan et al. ^[23] used ANCF to model the multilayer beam with a circular cross-section. Wang et al. ^[24] studied the dynamics of flexible multi-body systems with hybrid uncertain parameters by using the ANCF, to predict they are truly dynamic behaviors.

This paper mainly studied the flexible beam element, and briefly reviews the beam element based on the absolute nodal coordinate formulation. Shabana et al. ^[25] considered a two-dimensional beam element based on ANCF, which could efficiently solve the elastic force and force Jacobi. However, it is impossible to describe the complex beam elements subjected to torsion and shear. For this reason, Omar and Shabana and Omar ^[26] further proposed a two-dimensional shear deformation beam element that could describe the shear effect. García-Vallejo et al. ^[27] pointed out that the shear lock problem existed in the two-dimensional shear deformation beam element, and a new locking-free shear deformable finite element has been put forward. To accurately describe the force characteristics of the flexible spatial beam, Shabana and Yakoub ^[28] proposed a three-dimensional beam element model based on the absolute nodal coordinate formulation. García-Vallejo et al. ^[29] proposed an invariant matrix method to efficiently formulate the elastic force and the elastic force Jacobi, pointed out the nonlinear stiffness matrix can be simplified using the invariant matrix. Based on the principle of virtual work, Gerstmayr and Shabana ^[30] obtained the expression of the elastic force expressed by the first Piola-Kirchhoff stress tensor and position gradient tensor, the calculation efficiency was improved. Based on Gerstmayr's research, Liu et al. ^[31] derived the analytical formula that could accurately describe the elastic force Jacobian matrix, further improved the solving efficiency of the elastic force Jacobian.

However, it is difficult to apply a three-dimensional ANCF beam element, on the one hand, the elastic force is derived from the volume integral of the elastic energy using Green-Lagrange strain and Lame's constants ^[32]. On the other hand, the ANCF fully parametritis beam element uses linear interpolation at the transverse directions, therefore, the cross-section of the beam element must keep rectangular. Therefore, what is the research emphases of this paper are to extend the applicability of the ANCF on arbitrary section beams and its succinct solution form.

The contents of the paper are as follows. Section 2 introduces the absolute nodal coordinate theory and the discrete form of its shape function and proposes a position gradient construction method that takes initial configurations into account. Section 3 describes the single-integral form of the invariant matrix obtained by applying the discrete form of the shape function in the García-Vallejo method. The other single-integral form is described in Section 4 by applying the discrete form of the shape function in the Gerstmayr method. Section 5 proposes a high-order interpolation function for the Poisson locking problem. Some numerical simulations are given in Section 6, and the correctness of the model is verified. Section 7 provides concluding remarks.

2. Three-dimensional beam elements based on the ANCF

2.1. Absolute nodal coordinate formulation

In this subsection, the three-dimensional beam element for simplicity is shown in . Three-dimensional beam elements have 2 nodes, $\boldsymbol{e} = {\left[{{\boldsymbol{e}}_{\boldsymbol{i}}}^{\boldsymbol{T}}\quad{{\boldsymbol{e}}_{\boldsymbol{j}}}^{\boldsymbol{T}}\right]}^{T}$ . There are 24 degrees of freedom, each node consists of three position coordinates and nine slope coordinates.

Figure 1. Position vector of any point in the beam elements.

DownLoad: Full-Size Img PowerPoint

The co-ordinates of node $i$ are represented as,

$\boldsymbol{e}_{\boldsymbol{i}}^{\mathrm{T}} = \left[\begin{array}{llll} \boldsymbol{r}^{\mathrm{T}} & \boldsymbol{r}_{\boldsymbol{x}}^{\mathrm{T}} & \boldsymbol{r}_{\boldsymbol{y}}^{\mathrm{T}} & \boldsymbol{r}_{\boldsymbol{z}}^{\mathrm{T}} \end{array}\right]^{\mathrm{T}}$

(1)

where, the vector $\boldsymbol{r}$ denotes node $i$ global position vector in the beam elements, described as,

$\boldsymbol{r} = \left[\begin{array}{lll} X_{i} & Y_{i} & Z_{i} \end{array}\right]^{\mathrm{T}}$

(2)

And ${\left[{{\boldsymbol{r}}_{\boldsymbol{x}}}^{\mathrm{T}}\quad{{\boldsymbol{r}}_{\boldsymbol{y}}}^{\mathrm{T}}\quad{{\boldsymbol{r}}_{\boldsymbol{z}}}^{\mathrm{T}}\right]}^{\mathrm{T}}$ is position vector gradient, derivatives of the position vector $\boldsymbol{r}$ with respect to co-ordinates $x, y$ and $z$ ,

${\boldsymbol{r}}_{\boldsymbol{x}} = \frac{\partial \boldsymbol{r}}{\partial x} ;\quad {\boldsymbol{r}}_{\boldsymbol{y}} = \frac{\partial \boldsymbol{r}}{\partial y} ;\quad {\boldsymbol{r}}_{\boldsymbol{z}} = \frac{\partial \boldsymbol{r}}{\partial z}$

(3)

Because the finite element is a continuous system, the global position vector $\boldsymbol{r}$ of an arbitrary point $\boldsymbol{p}$ in beam elements is the function of the local parameters $x, y$ and $z$ . Using the shape function as follows,

${\boldsymbol{r}}_{\boldsymbol{p}} = \boldsymbol{S}\left({\bf{x}}\right)\boldsymbol{e} = \left[{\mathrm{S}}_{1}{\bf{I}}\quad{\mathrm{S}}_{2}{\bf{I}}\quad{\mathrm{S}}_{3}{\bf{I}}\quad{\mathrm{S}}_{4}{\bf{I}}\quad{\mathrm{S}}_{5}{\bf{I}}\quad{\mathrm{S}}_{6}{\bf{I}}\quad{\mathrm{S}}_{7}{\bf{I}}\quad{\mathrm{S}}_{8}{\bf{I}}\right]\boldsymbol{e}$

(4)

And ${\bf{I}}$ is a 3 $\times$ 3 identity matrix and the shape functions are defined as ^[30],

$\left\{\begin{array}{c}\begin{array}{c}{\mathrm{S}}_{1} = 1-3{\xi }^{2}+2{\xi }^{3}, \quad\quad\quad\quad\quad\mathrm{ }{\mathrm{S}}_{2} = L\left(\xi -2{\xi }^{2}+{\xi }^{3}\right)\\ {\mathrm{S}}_{3} = L\left(\eta -\xi \eta \right), \quad\quad\quad{\mathrm{S}}_{4} = L\left(\zeta -\xi \zeta \right) \end{array}\\ \begin{array}{c}{\mathrm{S}}_{5} = 3{\xi }^{2}-2{\xi }^{3}, \quad\quad\quad{\mathrm{S}}_{6} = L\left(-{\xi }^{2}+{\xi }^{3}\right) \\ {\mathrm{S}}_{7} = L\xi \eta , \quad\quad\quad {\mathrm{S}}_{8} = L\xi \zeta \quad\quad\quad \end{array}\end{array}\right.$

(5)

where,

$\xi = \frac{x}{L};\quad\eta = \frac{y}{L};\quad\zeta = \frac{z}{L}$

(6)

$L$ represents the beam element length. The shape function of the displacement field is composed of the cubic interpolation of $x$ and linear interpolation of $y$ and $z$ in the local coordinate system.

According to the theory of continuum mechanics, the position gradient tensor $\boldsymbol{G}$ is used to describe the position relationship between the initial configuration ${\boldsymbol{r}}_{{\bf{0}}}$ of the elastomer and the current configuration $\boldsymbol{r}$ after movement. The local coordinate system fixed on the material is denoted as $\alpha \left(\alpha = x, y, z\right)$ .

$\boldsymbol{G} = \frac{\partial \boldsymbol{r}}{\partial {\boldsymbol{r}}_{{\bf{0}}}} = \frac{\partial \boldsymbol{r}}{\partial \pmb{\mathsf{α}}}\frac{\partial \pmb{\mathsf{α}}}{\partial {\boldsymbol{r}}_{{\bf{0}}}}$

(7)

where,

${\boldsymbol{r}}_{{\bf{0}}} = \boldsymbol{S}\left(x, y, z\right){\boldsymbol{e}}_{{\bf{0}}}, \boldsymbol{r} = \boldsymbol{S}\left(x, y, z\right)\boldsymbol{e}$

(8)

And,

$\frac{\partial \boldsymbol{r}}{\partial \pmb{\mathsf{α}}} = \frac{\partial \boldsymbol{S}\left(x, y, z\right)\boldsymbol{e}}{\partial \pmb{\mathsf{α}}} = \left[{\boldsymbol{S}}_{, x}\boldsymbol{e}\quad{\boldsymbol{S}}_{, y}\boldsymbol{e}\quad{\boldsymbol{S}}_{, z}\boldsymbol{e}\right]$

(9)

${\boldsymbol{e}}_{{\bf{0}}}$ represents the generalized coordinates of the beam element nodes in the reference configuration,

$\frac{\partial \pmb{\mathsf{α}}}{\partial {\boldsymbol{r}}_{{\bf{0}}}} = {\left(\frac{\partial {\boldsymbol{r}}_{{\bf{0}}}}{\partial \pmb{\mathsf{α}}}\right)}^{-1} = {\boldsymbol{G}}_{{\bf{0}}}^{-1}$

(10)

Since ${\boldsymbol{e}}_{{\bf{0}}}$ and ${\boldsymbol{G}}_{{\bf{0}}}^{-1}$ are constant, and the position gradient is written as,

$\boldsymbol{G} = \frac{\partial \boldsymbol{r}}{\partial \pmb{\mathsf{α}}}{\boldsymbol{G}}_{{\bf{0}}}^{-1}$

(11)

In order to construct the subsequent invariant matrix, the position gradient tensor should be rewritten here. García-Vallejo et al. ^[29] proposed a position gradient construction method to introduce the initial configuration, here, a new construction method was proposed, the position vector ${\boldsymbol{r}}_{\boldsymbol{p}}$ , can denote as,

${\boldsymbol{r}}_{\boldsymbol{p}} = \overline{\boldsymbol{e}}\overline{\boldsymbol{S}}$

(12)

$\overline{\boldsymbol{S}} \; \in {\mathcal{R}}^{8\times 1}$ , $\overline{\boldsymbol{e}}\in {\mathcal{R}}^{3\times 8}$ , the specific form is as follows.

$\left\{\begin{array}{c}\overline{\boldsymbol{S}} = {\left[{S}_{1}\quad{S}_{2}\quad{S}_{3}\quad{S}_{4}\quad{S}_{5}\quad{S}_{6} {S}_{7}\quad{S}_{8}\right]}^{\mathrm{T}} \\ \overline{\boldsymbol{e}} = \left[\begin{array}{ccc}\begin{array}{ccc}{e}_{1}& {e}_{4}& {e}_{7}\\ {e}_{2}& {e}_{5}& {e}_{8}\\ {e}_{3}& {e}_{6}& {e}_{9}\end{array}& \begin{array}{ccc}{e}_{10}& {e}_{13}& {e}_{16}\\ {e}_{11}& {e}_{14}& {e}_{17}\\ {e}_{12}& {e}_{15}& {e}_{18}\end{array}& \begin{array}{cc}{e}_{19}& {e}_{22}\\ {e}_{20}& {e}_{23}\\ {e}_{21}& {e}_{24}\end{array}\end{array}\right]\end{array}\right.$

(13)

Therefore, the position gradient is written as follows,

$\boldsymbol{G} = \overline{\boldsymbol{e}}\frac{\partial \overline{\boldsymbol{S}}}{\partial \pmb{\mathsf{α}}}{\left(\frac{\partial {\boldsymbol{r}}_{{\bf{0}}}}{\partial \pmb{\mathsf{α}}}\right)}^{-1} = \overline{\boldsymbol{e}}\boldsymbol{H}$

(14)

and

$\boldsymbol{H} = \frac{\partial \overline{\boldsymbol{S}}}{\partial \pmb{\mathsf{α}}}{\left(\frac{\partial {\boldsymbol{r}}_{{\bf{0}}}}{\partial \pmb{\mathsf{α}}}\right)}^{-1}$

(15)

$\boldsymbol{H}\in {\mathcal{R}}^{8\times 3}$ contains the reference configuration information of the beam element. According to the expression, $\boldsymbol{H}$ is a constant matrix, and the Green strain tensor can be expressed as,

$\left[\pmb{\mathsf{ε}}\right] = \frac{1}{2}({\boldsymbol{G}}^{\boldsymbol{T}}\boldsymbol{G}-{\bf{I}})$

(16)

The components of symmetric strain tensor are,

${\varepsilon }_{ij} = \left\{\begin{array}{c}\frac{1}{2}\left({{G}_{, i}}^{T}{G}_{, j}-1\right)\quad\quad(i = j, i, j = \mathrm{1, 2}, 3)\\ \frac{1}{2}{{G}_{, i}}^{T}{G}_{, j}\quad\quad\quad\quad(i\ne j, i, j = \mathrm{1, 2}, 3)\end{array}\right.$

(17)

${G}_{, i}$ represents the $i$ th ( $i = \mathrm{1, 2}, 3$ ) column of the position gradient.

${G}_{, i} = \overline{\boldsymbol{e}}{H}_{, i}$

(18)

${H}_{, i}$ is the $i$ th ( $i = \mathrm{1, 2}, 3$ ) column of the matrix H, and,

${H}_{, i} = {\left[{h}_{1i}\quad{h}_{2i}\quad{h}_{3i}\quad{h}_{4i}\quad{h}_{5i}\quad{h}_{6i}\quad{h}_{7i}\quad{h}_{8i}\right]}^{\mathrm{T}}$

(19)

Introduce ${\overline{\boldsymbol{H}}}_{i}\in {\mathcal{R}}^{3\times 24}$ ,

${\overline{\boldsymbol{H}}}_{i} = \left[{h}_{1i}{\bf{I}}\quad{h}_{2i}{\bf{I}}\quad{h}_{3i}{\bf{I}}\quad{h}_{4i}{\bf{I}}\quad{h}_{5i}{\bf{I}}\quad{h}_{6i}{\bf{I}}\quad{h}_{7i}{\bf{I}}\quad{h}_{8i}{\bf{I}}\right]$

(20)

Substitute into Eq (18), a new form of position gradient column vector is constructed.

${G}_{, i} = {\overline{\boldsymbol{H}}}_{i}\boldsymbol{e}$

(21)

Hence, the six different components of the symmetric strain tensor are,

$\left\{\begin{array}{c}{\mathtt{ε }}_{11} = \frac{1}{2}\left({\boldsymbol{e}}^{\mathrm{T}}{{\overline{\boldsymbol{H}}}_{1}}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{1}\boldsymbol{e}-1\right), {\mathtt{ε }}_{22} = \frac{1}{2}\left({\boldsymbol{e}}^{\mathrm{T}}{{\overline{\boldsymbol{H}}}_{2}}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{2}\boldsymbol{e}-1\right), {\mathtt{ε }}_{33} = \frac{1}{2}\left({\boldsymbol{e}}^{\mathrm{T}}{{\overline{\boldsymbol{H}}}_{3}}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{3}\boldsymbol{e}-1\right)\\ {\mathtt{ε }}_{12} = \frac{1}{2}{\boldsymbol{e}}^{\mathrm{T}}{{\overline{\boldsymbol{H}}}_{1}}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{2}\boldsymbol{e}, {\mathtt{ε }}_{13} = \frac{1}{2}{\boldsymbol{e}}^{\mathrm{T}}{{\overline{\boldsymbol{H}}}_{1}}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{3}\boldsymbol{e}\boldsymbol{ }, {\mathtt{ε }}_{23} = \frac{1}{2}{\boldsymbol{e}}^{T}{{\overline{\boldsymbol{H}}}_{2}}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{3}\boldsymbol{e}\end{array}\right.$

(22)

While, based on the principle of virtual work, the expression of the element elastic force can be obtained.

$\boldsymbol{Q}\left(\boldsymbol{e}\right) = -\frac{\partial U}{\partial \boldsymbol{e}} = -\boldsymbol{K}\left(\boldsymbol{e}\right)\boldsymbol{e}$

(23)

2.2. Discrete shape function

Obviously, the shape function $\overline{\boldsymbol{S}}$ , can be rewritten as follows,

$\overline{\boldsymbol{S}} = \boldsymbol{f}\left(\xi \right)+\eta \boldsymbol{g}\left(\xi \right)+\zeta \boldsymbol{h}\left(\xi \right)$

(24)

here,

$\left\{\begin{array}{c}\boldsymbol{f}\left(\xi \right) = {\left[{\mathrm{S}}_{1} \quad {\mathrm{S}}_{2} \quad 0 \quad 0 \quad {\mathrm{S}}_{5} \quad {\mathrm{S}}_{6} \quad 0 \quad 0\right]}^{\mathrm{T}} \\ \boldsymbol{g}\left(\xi \right) = {\left[0 \quad 0 \quad L \left(1-\xi \right) 0 \quad 0 \quad 0 \quad L \xi \quad 0\right]}^{\mathrm{T}} \\ \boldsymbol{h}\left(\xi \right) = {\left[0 \quad 0 \quad 0 \quad L \left(1-\xi \right) \quad 0 \quad 0 \quad 0 \quad L \xi \right]}^{\mathrm{T}} \end{array}\right.$

(25)

${\bf{I}}$ is a 3 $\times$ 3 identity matrix, ${\bf{O}}$ is a 3 $\times$ 3 zero matrix. $\boldsymbol{f}\left(\xi \right), \boldsymbol{g}\left(\xi \right), \boldsymbol{h}\left(\xi \right)$ just the function of $\xi$ , have nothing to do with $\eta$ and $\zeta$ . Based on Eq (24), the partial derivative of form function $\overline{\boldsymbol{S}}$ to local coordinates is,

$\frac{\partial \overline{\boldsymbol{S}}}{\partial \pmb{\mathsf{α}}} = \left[{\boldsymbol{f}}_{1}\left(\xi \right) \quad {\boldsymbol{f}}_{2}\left(\xi \right) \quad {\boldsymbol{f}}_{3}\left(\xi \right)\right]+\eta \left[{\boldsymbol{g}}_{1} \quad {{\bf{O}}}_{1} \quad {{\bf{O}}}_{1}\right]+\zeta \left[{\boldsymbol{h}}_{1} \quad {{\bf{O}}}_{1} \quad {{\bf{O}}}_{1}\right]$

(26)

The detailed expressions of ${\boldsymbol{f}}_{1}\left(\xi \right)$ , ${\boldsymbol{f}}_{2}\left(\xi \right)$ , ${\boldsymbol{f}}_{3}\left(\xi \right)$ and ${\boldsymbol{g}}_{1}$ , ${\boldsymbol{h}}_{1}$ are given in Appendix 1 Eq (A3), ${{\bf{O}}}_{1}$ is a $8\times 1$ zero vector. While, Eq (15) can be expressed as,

$\boldsymbol{H} = {\boldsymbol{w}}_{\bf 1}+\eta {\boldsymbol{w}}_{\bf 2}+\zeta {\boldsymbol{w}}_{\bf 3}$

(27)

Here, ${\boldsymbol{w}}_{\bf 1}, {\boldsymbol{w}}_{\bf 2}, {\boldsymbol{w}}_{\bf 3}\in {\mathcal{R}}^{8\times 3}$ ,

$\left\{\begin{array}{c}{\boldsymbol{w}}_{\bf 1} = \left[{\boldsymbol{f}}_{1}\left(\xi \right) \quad {\boldsymbol{f}}_{2}\left(\xi \right) \quad {\boldsymbol{f}}_{3}\left(\xi \right)\right]{\left(\frac{\partial {\boldsymbol{r}}_{{\bf{0}}}}{\partial \pmb{\mathsf{α}}}\right)}^{-1}\\ {\boldsymbol{w}}_{\bf 2} = \left[{\boldsymbol{g}}_{1} \quad {{\bf{O}}}_{1} \quad {{\bf{O}}}_{1}\right]{\left(\frac{\partial {\boldsymbol{r}}_{{\bf{0}}}}{\partial \pmb{\mathsf{α}}}\right)}^{-1} \quad \quad \quad \quad \quad \quad \quad \quad \\ {\boldsymbol{w}}_{\bf 3} = \left[{\boldsymbol{h}}_{1} \quad {{\bf{O}}}_{1} \quad {{\bf{O}}}_{1}\right]{\left(\frac{\partial {\boldsymbol{r}}_{{\bf{0}}}}{\partial \pmb{\mathsf{α}}}\right)}^{-1} \quad \quad \quad \quad \quad \quad \quad \quad \end{array}\right.$

(28)

And, ${H}_{, i}(i = \mathrm{1, 2}, 3)$

${H}_{, i} = {w}_{1, i}+\eta {w}_{2, i}+\zeta {w}_{3, i}$

(29)

From Eq (20),

${\overline{\boldsymbol{H}}}_{i} = {\overline{\boldsymbol{w}}}_{1i}+\eta {\overline{\boldsymbol{w}}}_{2i}+\zeta {\overline{\boldsymbol{w}}}_{3i}$

(30)

3. Application of discrete shape function to invariant matrix method

García-Vallejo pointed out that the nonlinear stiffness matrix $\boldsymbol{K}\left(\boldsymbol{e}\right)$ in Eq (23) can be divided into two parts,

$\boldsymbol{K}\left(\boldsymbol{e}\right) = {\boldsymbol{K}}_{1}+{\boldsymbol{K}}_{2}\left(\boldsymbol{e}\right) \\ {\boldsymbol{K}}_{1} = -\frac{3\mathtt{λ }+2\mathtt{μ }}{2}\sum\limits_{\alpha = 1}^{3}\int {{\overline{\boldsymbol{H}}}_{\alpha }}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{\alpha }\mathrm{d}V \\ {\boldsymbol{K}}_{2}\left(\boldsymbol{e}\right) = \frac{\mathtt{λ }+2\mathtt{μ }}{2}\sum\limits_{\alpha = 1}^{3}\int {{\overline{\boldsymbol{H}}}_{\alpha }}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{\alpha }\boldsymbol{e}{\boldsymbol{e}}^{\mathrm{T}}{{\overline{\boldsymbol{H}}}_{\alpha }}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{\alpha }\mathrm{d}V \\ +\sum\limits_{\alpha = 1}^{3}\sum\limits_{\beta = 1, \alpha \ne \beta }^{3}\int (\frac{\mathtt{λ }}{2}{{\overline{\boldsymbol{H}}}_{\alpha }}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{\alpha }\boldsymbol{e}{\boldsymbol{e}}^{\mathrm{T}}{{\overline{\boldsymbol{H}}}_{\beta }}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{\beta }+\mathtt{μ }{{\overline{\boldsymbol{H}}}_{\alpha }}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{\beta }\boldsymbol{e}{\boldsymbol{e}}^{\mathrm{T}}{{\overline{\boldsymbol{H}}}_{\alpha }}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{\beta })\mathrm{d}V$

(31)

The constant stiffness matrix ${\boldsymbol{K}}_{1}$ can be calculated and stored in the pre-processing process. The nonlinear stiffness matrix ${\boldsymbol{K}}_{2}\left(\boldsymbol{e}\right)$ , it can be solved by the invariant matrix, so the elements in the ${\boldsymbol{K}}_{2}\left(\boldsymbol{e}\right)$ matrix is expressed by the invariant matrix,

${\left({\boldsymbol{K}}_{2}\right(\boldsymbol{e}\left)\right)}_{ij} = {\boldsymbol{e}}^{\mathrm{T}}{\boldsymbol{C}}_{{\boldsymbol{K}}_{2}}^{ij}\boldsymbol{e} \\ {\boldsymbol{C}}_{{\boldsymbol{K}}_{2}}^{ij} = \frac{\mathtt{λ }+2\mathtt{μ }}{2}\sum\limits_{\alpha = 1}^{3}\int ({{{\overline{\boldsymbol{H}}}_{\alpha }}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{\alpha })}_{i, }^{\mathrm{T}}{\left({{\overline{\boldsymbol{H}}}_{\alpha }}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{\alpha }\right)}_{j, }\mathrm{d}V \\ +\sum\limits_{\alpha = 1}^{3}\sum\limits_{\beta = 1, \alpha \ne \beta }^{3}\int \left(\frac{\mathtt{λ }}{2}\right({{{\overline{\boldsymbol{H}}}_{\alpha }}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{\alpha })}_{i, }^{\mathrm{T}}{\left({{\overline{\boldsymbol{H}}}_{\beta }}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{\beta }\right)}_{j, }+\mathtt{μ }\left({{{\overline{\boldsymbol{H}}}_{\alpha }}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{\beta })}_{i, }^{\mathrm{T}}{\left({{\overline{\boldsymbol{H}}}_{\beta }}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{\alpha }\right)}_{j, }\right)\mathrm{d}V$

(32)

By substituting the discrete form of the shape function of Eq (30), ${\boldsymbol{K}}_{1}$ and ${\boldsymbol{K}}_{2}\left(\boldsymbol{e}\right)$ can be obtained. Please refer to Appendix 2 for the specific derivation process.

${\boldsymbol{K}}_{1} = -\frac{3\mathtt{λ }+2\mathtt{μ }}{2}\sum\limits_{\alpha = 1}^{3}\int \left({\boldsymbol{A}}_{\alpha }+\eta {\boldsymbol{B}}_{\alpha }+\zeta {\boldsymbol{C}}_{\alpha }+\eta \zeta {\boldsymbol{D}}_{\alpha }+{\eta }^{2}{\boldsymbol{E}}_{\alpha }+{\zeta }^{2}{\boldsymbol{F}}_{\alpha }\right)\mathrm{d}V \\ {\left({\boldsymbol{K}}_{2}\right(\boldsymbol{e}\left)\right)}_{ij} = {\boldsymbol{e}}^{\mathrm{T}}\left\{\sum\limits_{k = 1}^{15}\int {{\mathit{ \Pi}} }_{k}\left[\begin{array}{c}{\boldsymbol{N}}_{1}^{k}+\eta {\boldsymbol{N}}_{2}^{k}+\zeta {\boldsymbol{N}}_{3}^{k}+\eta \zeta {\boldsymbol{N}}_{4}^{k}+{\eta }^{2}{\boldsymbol{N}}_{5}^{k}\\ +{\zeta }^{2}{\boldsymbol{N}}_{6}^{k}+\eta {\zeta }^{2}{\boldsymbol{N}}_{7}^{k}+{\eta }^{2}\zeta {\boldsymbol{N}}_{8}^{k}+{\eta }^{3}{\boldsymbol{N}}_{9}^{k}+{\zeta }^{3}{\boldsymbol{N}}_{10}^{k}\\ +{\eta }^{2}{\zeta }^{2}{\boldsymbol{N}}_{11}^{k}+\eta {\zeta }^{3}{\boldsymbol{N}}_{12}^{k}+{\eta }^{3}\zeta {\boldsymbol{N}}_{13}^{k}+{\eta }^{4}{\boldsymbol{N}}_{14}^{k}+{\zeta }^{4}{\boldsymbol{N}}_{15}^{k}\end{array}\right]\mathrm{d}V\right\}\boldsymbol{e}$

(33)

Orzechowski et al. ^[33] pointed out the ANCF is valid theoretically for the arbitrarily shaped beam elements. Change the integral form of ${\boldsymbol{K}}_{1}$ and ${\boldsymbol{K}}_{2}\left(\boldsymbol{e}\right)$ , and write it as,

${\boldsymbol{K}}_{1} = -\frac{3\mathtt{λ }+2\mathtt{μ }}{2}\sum\limits_{\alpha = 1}^{3}{\int }_{{\bf{0}}}^{L}dx\int \left({\boldsymbol{A}}_{\alpha }+\eta {\boldsymbol{B}}_{\alpha }+\zeta {\boldsymbol{C}}_{\alpha }+\eta \zeta {\boldsymbol{D}}_{\alpha }+{\eta }^{2}{\boldsymbol{E}}_{\alpha }+{\zeta }^{2}{\boldsymbol{F}}_{\alpha }\right)\mathrm{d}\mathrm{A} \\ {\left({\boldsymbol{K}}_{2}\right(\boldsymbol{e}\left)\right)}_{ij}\\ = {\boldsymbol{e}}^{\mathrm{T}}\left\{\sum\limits_{k = 1}^{15}{\int }_{{\bf{0}}}^{L}dx\int {{\mathit{ \Pi}} }_{k}\left[\begin{array}{c}{\boldsymbol{N}}_{1}^{k}+\eta {\boldsymbol{N}}_{2}^{k}+\zeta {\boldsymbol{N}}_{3}^{k}+\eta \zeta {\boldsymbol{N}}_{4}^{k}+{\eta }^{2}{\boldsymbol{N}}_{5}^{k}\\ +{\zeta }^{2}{\boldsymbol{N}}_{6}^{k}+\eta {\zeta }^{2}{\boldsymbol{N}}_{7}^{k}+{\eta }^{2}\zeta {\boldsymbol{N}}_{8}^{k}+{\eta }^{3}{\boldsymbol{N}}_{9}^{k}+{\zeta }^{3}{\boldsymbol{N}}_{10}^{k}\\ +{\eta }^{2}{\zeta }^{2}{\boldsymbol{N}}_{11}^{k}+\eta {\zeta }^{3}{\boldsymbol{N}}_{12}^{k}+{\eta }^{3}\zeta {\boldsymbol{N}}_{13}^{k}+{\eta }^{4}{\boldsymbol{N}}_{14}^{k}+{\zeta }^{4}{\boldsymbol{N}}_{15}^{k}\end{array}\right]\mathrm{d}\mathrm{A}\right\}\boldsymbol{e}$

(34)

In practice, the beam element commonly used in engineering is the Euler-Bernoulli beam with a symmetrical section and a small section size. Considering this situation, the cross-section integral can be simplified, that is, the odd term integral is zero, even higher-order integral tends to zero. The greatest advantage of this is that the stiffness matrix ${\boldsymbol{K}}_{1}$ , ${\boldsymbol{K}}_{2}\left(\boldsymbol{e}\right)$ is simplified to the single integral form which is only related to the area and moment of inertia of the section.

${\boldsymbol{K}}_{1} = -\frac{3\lambda +2\mu }{2}\sum\limits_{\alpha = 1}^{3}{\int }_{{\bf{0}}}^{L}\left[\mathrm{A}{\boldsymbol{A}}_{\alpha }+{\mathrm{I}}_{z}{\boldsymbol{E}}_{\alpha }+{\mathrm{I}}_{y}{\boldsymbol{F}}_{\alpha }\right]\mathrm{d}x \\ {\left({\boldsymbol{K}}_{2}\left(\boldsymbol{e}\right)\right)}_{ij} = {\boldsymbol{e}}^{\mathrm{T}}\left\{\sum\limits_{k = 1}^{15}{\int }_{{\bf{0}}}^{L}{{\mathit{ \Pi}} }_{k}[\mathrm{A}{\boldsymbol{N}}_{1}^{k}+{\mathrm{I}}_{z}{\boldsymbol{N}}_{5}^{k}+{\mathrm{I}}_{y}{\boldsymbol{N}}_{6}^{k}]dx\right\}\boldsymbol{e}$

(35)

$\mathrm{A}$ is the section area of the beam element with arbitrary section, ${\mathrm{I}}_{z}$ is the moment of inertia around the $z$ -axis of the beam section, and ${\mathrm{I}}_{y}$ is the moment of inertia around the $y$ -axis of the beam section.

The Jacobian of the elastic force can also be derived by the invariant matrix method,

$\boldsymbol{J}\left(\boldsymbol{e}\right) = \frac{\partial \boldsymbol{Q}\left(\boldsymbol{e}\right)}{\partial \boldsymbol{e}}$

(36)

And the element of the Jacobian can express as,

${J}_{ik} = -\sum\limits_{j = 1}^{n}({{\boldsymbol{K}}_{1})}_{ij}{\delta }_{jk}-\sum\limits_{j = 1}^{n}\left({\boldsymbol{e}}^{T}{\boldsymbol{C}}_{{\boldsymbol{K}}_{2}}^{ij}\boldsymbol{e}\right){\delta }_{jk}-\sum\limits_{j = 1}^{n}\frac{\partial \left({\boldsymbol{e}}^{T}{\boldsymbol{C}}_{{\boldsymbol{K}}_{2}}^{ij}\boldsymbol{e}\right)}{\partial {e}_{k}}{e}_{j}$

(37)

4. Application of discrete shape function to Gerstmayr method

The Gerstmayr method is based on the first Piola-Kirchhoff stress tensor and position gradient tensor ^[31], the virtual work of the elastic force is written as,

$\delta W = \int \mathcal{S}:\delta \mathcal{E}dV = \int \mathcal{P}:\delta \boldsymbol{G}dV$

(38)

here, $\mathcal{S}$ is the second Piola-Kirchhoff stress tensor, $\mathcal{E}$ is the Green strain tensor, $\mathcal{P}$ is the first Piola-Kirchhoff stress tensor, and the component of the elastic force can be expressed as,

${Q}_{i}^{\boldsymbol{e}} = -\int {\mathcal{P}}_{kl}\frac{\partial {G}_{kl}}{\partial {e}_{i}}dV$

(39)

Introduce,

${e}_{i} = {\overline{e}}_{rs}, i = 3\left(s-1\right)+r$

(40)

Equation (39) can be expressed as,

${Q}_{i}^{\boldsymbol{e}} = -\int {\mathcal{P}}_{kl}\frac{\partial {G}_{kl}}{\partial {\overline{e}}_{rs}}dV = -\int {\mathcal{P}}_{kl}\frac{\partial \left({\overline{e}}_{km}{H}_{ml}\right)}{\partial {\overline{e}}_{rs}}dV \\ = -\int {\mathcal{P}}_{kl}{\delta }_{kr}{\delta }_{ms}{H}_{ml}dV = -\int {\mathcal{P}}_{rl}{H}_{sl}dV$

(41)

${\delta }_{kr}$ is Kronecker $\delta$ , $\mathcal{P}$ can be represented by $\boldsymbol{G}$ and $\mathcal{S}$ ,

${\mathcal{P}}_{rl} = {G}_{rn}{\mathcal{S}}_{nl} = \left({\overline{e}}_{rt}{H}_{tn}\right)\left[\frac{1}{2}{\mathcal{D}}_{nlab}({{G}_{ac}}^{\mathrm{T}}{G}_{cb}-{\mathcal{J}}_{ab})\right]$

(42)

here, ${\mathcal{D}}_{nlab}$ represents the component of the fourth-order elastic tensor, ${\mathcal{J}}_{ab}$ represents the component of the second-order unit tensor,

$\mathcal{D} = \lambda {\delta }_{ij}{\delta }_{kl}{\overrightarrow{\boldsymbol{e}}}_{i} \otimes {\overrightarrow{\boldsymbol{e}}}_{j} \otimes {\overrightarrow{\boldsymbol{e}}}_{k} \otimes {\overrightarrow{\boldsymbol{e}}}_{l}+2\mu {\delta }_{ik}{\delta }_{jl}{\overrightarrow{\boldsymbol{e}}}_{i} \otimes {\overrightarrow{\boldsymbol{e}}}_{k} \otimes {\overrightarrow{\boldsymbol{e}}}_{j} \otimes {\overrightarrow{\boldsymbol{e}}}_{l} \\ \mathcal{J} = {\delta }_{ij}{\overrightarrow{\boldsymbol{e}}}_{i} \otimes {\overrightarrow{\boldsymbol{e}}}_{j}$

(43)

Substituting Eq (42) into Eq (41),

${Q}_{i}^{\boldsymbol{e}} = -\int \left({\overline{e}}_{rt}{H}_{tn}\right)\left[\frac{1}{2}{\mathcal{D}}_{nlab}({{G}_{ac}}^{\mathrm{T}}{G}_{cb}-{\mathcal{J}}_{ab})\right]{H}_{sl}dV \\ = -\frac{1}{2}\int {\mathcal{D}}_{nlab}\left({\overline{e}}_{rt}{H}_{tn}\right)({\overline{e}}_{cf}{H}_{fa}{\overline{e}}_{cv}{H}_{vb}-{\mathcal{J}}_{ab}){H}_{sl}dV \\ = -\frac{1}{2}\int {\mathcal{D}}_{nlab}{H}_{fa}{H}_{vb}{H}_{tn}{H}_{sl}dV{\overline{e}}_{rt}{\overline{e}}_{cf}{\overline{e}}_{cv}+\frac{1}{2}\int {\mathcal{D}}_{nlab}{\mathcal{J}}_{ab}{H}_{tn}{H}_{sl}dV{\overline{e}}_{rt} \\ = \frac{1}{2}{K}_{ts}^{\mathrm{I}}{\overline{e}}_{rt}-\frac{1}{2}{K}_{fvts}^{\mathrm{I}\mathrm{I}}{\overline{e}}_{rt}{\overline{e}}_{cf}{\overline{e}}_{cv}$

(44)

Introduce the discrete form of the shape function to Eq (27), and write ${w}_{i}(i = \mathrm{1, 2}, 3)$ as ${w}^{i}(i = \mathrm{1, 2}, 3)$ for convenience,

$\left\{\begin{array}{c}\begin{array}{c}{H}_{fa} = {{w}^{1}}_{fa}+\eta {{w}^{2}}_{fa}+\zeta {{w}^{3}}_{fa}\\ {H}_{vb} = {{w}^{1}}_{vb}+\eta {{w}^{2}}_{vb}+\zeta {{w}^{3}}_{vb}\end{array}\\ \begin{array}{c}{H}_{tn} = {{w}^{1}}_{tn}+\eta {{w}^{2}}_{tn}+\zeta {{w}^{3}}_{tn} \\ {H}_{sl} = {{w}^{1}}_{sl}+\eta {{w}^{2}}_{sl}+\zeta {{w}^{3}}_{sl} \end{array}\end{array}\right.$

(45)

So,

${H}_{tn}{H}_{sl} = \mathbb{A}+\eta \mathbb{B}+\zeta \mathbb{C}+\eta \zeta \mathbb{D}+{\eta }^{2}\mathbb{E}+{\zeta }^{2}\mathbb{F} \\ {H}_{fa}{H}_{vb}{H}_{tn}{H}_{sl} = {\mathbb{N}}_{1}+\eta {\mathbb{N}}_{2}+\zeta {\mathbb{N}}_{3}+\eta \zeta {\mathbb{N}}_{4}+{\eta }^{2}{\mathbb{N}}_{5} \\ +{\zeta }^{2}{\mathbb{N}}_{6}+\eta {\zeta }^{2}{\mathbb{N}}_{7}+{\eta }^{2}\zeta {\mathbb{N}}_{8}+{\eta }^{3}{\mathbb{N}}_{9}+{\zeta }^{3}{\mathbb{N}}_{10} \\ +{\eta }^{2}{\zeta }^{2}{\mathbb{N}}_{11}+\eta {\zeta }^{3}{\mathbb{N}}_{12}+{\eta }^{3}\zeta {\mathbb{N}}_{13}+{\eta }^{4}{\mathbb{N}}_{14}+{\zeta }^{4}{\mathbb{N}}_{15}$

(46)

The specific derivation process can refer to Appendix 3. When the beam element is a symmetrical cross-section and small section size, the Eq (46) can be simplified to,

${H}_{tn}{H}_{sl} = \mathbb{A}+{\eta }^{2}\mathbb{E}+{\zeta }^{2}\mathbb{F}, {H}_{fa}{H}_{vb}{H}_{tn}{H}_{sl} = {\mathbb{N}}_{1}+{\eta }^{2}{\mathbb{N}}_{5}+{\zeta }^{2}{\mathbb{N}}_{6}$

(47)

and

${K}_{ts}^{\mathrm{I}} = \int {\mathcal{D}}_{nlab}{\mathcal{J}}_{ab}\left(\mathbb{A}+{\eta }^{2}\mathbb{E}+{\zeta }^{2}\mathbb{F}\right)dV = {\int }_{{{0}}}^{l}{\mathcal{D}}_{nlab}{\mathcal{J}}_{ab}[\mathrm{A}\mathbb{A}+{\mathrm{I}}_{\mathrm{z}}\mathbb{E}+{\mathrm{I}}_{y}\mathbb{F}]\mathrm{d}x \\ {K}_{fvts}^{\mathrm{I}\mathrm{I}} = \int {\mathcal{D}}_{nlab}\left({\mathbb{N}}_{1}+{\eta }^{2}{\mathbb{N}}_{5}+{\zeta }^{2}{\mathbb{N}}_{6}\right)dV = {\int }_{{{0}}}^{l}{\mathcal{D}}_{nlab}[\mathrm{A}{\mathbb{N}}_{1}+{\mathrm{I}}_{\mathrm{z}}{\mathbb{N}}_{5}+{\mathrm{I}}_{y}{\mathbb{N}}_{6}]\mathrm{d}x$

(48)

${K}_{ts}^{\mathrm{I}}$ and ${K}_{fvts}^{\mathrm{I}\mathrm{I}}$ are constant stiffness matrices that can be stored in the computer during the pre-processing, compared with the García-Vallejo method, the storage scale of the Gerstmayr method is significantly reduced. Jacobi of elastic force is often used in numerical analysis, Liu et al. ^[31] further derived the Jacobi analytic formula for elastic force on the basis of Gerstmayr.

${J}_{ij} = \frac{\partial {Q}_{i}^{\boldsymbol{e}}}{\partial {e}_{j}}$

(49)

Introduce,

${e}_{j} = {\overline{e}}_{\alpha \beta }, j = 3\left(\beta -1\right)+\alpha$

(50)

The analytic formula is,

$\frac{\partial {Q}_{i}^{\boldsymbol{e}}}{\partial {e}_{j}} = \frac{\partial \left(\frac{1}{2}{K}_{ts}^{\mathrm{I}}{\overline{e}}_{rt}-\frac{1}{2}{K}_{fvts}^{\mathrm{I}\mathrm{I}}{\overline{e}}_{rt}{\overline{e}}_{cf}{\overline{e}}_{cv}\right)}{\partial {\overline{e}}_{\alpha \beta }} \\ = \frac{1}{2}{K}_{ts}^{\mathrm{I}}{\delta }_{r\alpha }{\delta }_{t\beta }-\frac{1}{2}{K}_{fvts}^{\mathrm{I}\mathrm{I}}{\delta }_{r\alpha }{\delta }_{t\beta }{\overline{e}}_{cf}{\overline{e}}_{cv}-\frac{1}{2}{K}_{fvts}^{\mathrm{I}\mathrm{I}}{\overline{e}}_{rt}({\delta }_{c\alpha }{\delta }_{f\beta }{\overline{e}}_{cv}+{\overline{e}}_{cf}{\delta }_{c\alpha }{\delta }_{v\beta })$

(51)

5. Poisson locking

Since most of the interpolating polynomials into ANCF beam elements in the transverse direction only use the first-order polynomial, the cross-section is still planar, thus prone to various locking effects. What is worth mentioning is that using higher-order interpolation polynomial can avoid locking effects ^[34].

$\mathbb{r} = \left[\begin{array}{c}{a}_{{\bf{0}}}+{a}_{1}x+{a}_{2}y+{a}_{3}z+{a}_{4}xy+{a}_{5}xz+{a}_{6}{x}^{2}+{a}_{7}{x}^{3}+{a}_{8}{x}^{2}y+{a}_{9}{x}^{2}z\\ {b}_{{\bf{0}}}+{b}_{1}x+{b}_{2}y+{b}_{3}z+{b}_{4}xy+{b}_{5}xz+{b}_{6}{x}^{2}+{b}_{7}{x}^{3}+{b}_{8}{x}^{2}y+{b}_{9}{x}^{2}z\\ {c}_{{\bf{0}}}+{c}_{1}x+{c}_{2}y+{c}_{3}z+{c}_{4}xy+{c}_{5}xz+{c}_{6}{x}^{2}+{c}_{7}{x}^{3}+{c}_{8}{x}^{2}y+{c}_{9}{x}^{2}z\end{array}\right]$

(52)

The position vector expressed in terms of the nodal coordinates is given by,

$\mathbb{r} = \overline{\mathbb{e}}\overline{\mathbb{S}}$

(53)

$\overline{\mathbb{S}} \; \in {\mathcal{R}}^{10\times 1}$ , $\overline{\mathbb{e}}\in {\mathcal{R}}^{3\times 10}$ , the specific form is as follows.

$\left\{\begin{array}{c}\overline{\mathbb{S}} = {\left[{S}_{1}\quad{S}_{2}\quad{S}_{3}\quad{S}_{4}\quad{S}_{5}\quad{S}_{6}\quad{S}_{7}\quad{S}_{8}\quad{S}_{9}\quad{S}_{10}\right]}^{\mathrm{T}} \\ \overline{\mathbb{e}} = \left[\begin{array}{ccc}\begin{array}{ccc}{e}_{1}& {e}_{4}& {e}_{7}\\ {e}_{2}& {e}_{5}& {e}_{8}\\ {e}_{3}& {e}_{6}& {e}_{9}\end{array}& \begin{array}{ccc}{e}_{10}& {e}_{13}& {e}_{16}\\ {e}_{11}& {e}_{14}& {e}_{17}\\ {e}_{12}& {e}_{15}& {e}_{18}\end{array}& \begin{array}{cc}{e}_{19}& {e}_{22} \\ {e}_{20}& {e}_{23} \\ {e}_{21}& {e}_{24} \end{array}\end{array}\begin{array}{cc}{e}_{25}& {e}_{28}\\ {e}_{26}& {e}_{29}\\ {e}_{27}& {e}_{30}\end{array}\right]\end{array}\right.$

(54)

here,

$\left\{\begin{array}{c}\begin{array}{c}\begin{array}{c}{\mathrm{S}}_{1} = 1-3{\xi }^{2}+2{\xi }^{3}, \quad\quad\quad\quad{\mathrm{S}}_{2} = L\left(\xi -2{\xi }^{2}+{\xi }^{3}\right)\\ {\mathrm{S}}_{3} = L\left(1-{\xi }^{2}\right)\eta , \quad\quad\quad\quad{\mathrm{S}}_{4} = L\left(1-{\xi }^{2}\right) \zeta \end{array}\\ \begin{array}{c}{\mathrm{S}}_{5} = 3{\xi }^{2}-2{\xi }^{3}, \quad\quad\quad\quad {\mathrm{S}}_{6} = L\left(-{\xi }^{2}+{\xi }^{3}\right) \\ {\mathrm{S}}_{7} = L{\xi }^{2}\eta , \quad\quad\quad\quad {\mathrm{S}}_{8} = L{\xi }^{2}\zeta \quad\quad\quad\quad \end{array}\end{array}\\ {\mathrm{S}}_{9} = \left(\xi -{\xi }^{2}\right){L}^{2}\eta , \quad\quad\quad\quad{\mathrm{S}}_{10} = \left(\xi -{\xi }^{2}\right){L}^{2}\zeta \end{array}\right.$

(55)

The discrete form of shape function is,

$\overline{\mathbb{S}} = \mathbb{f}\left(\xi \right)+\eta \mathbb{g}\left(\xi \right)+\zeta \mathbb{h}\left(\xi \right)$

(56)

And,

$\left\{\begin{array}{c}\mathbb{f}\left(\xi \right) = {\left[{\mathrm{S}}_{1}\quad{\mathrm{S}}_{2}\quad0\quad0\quad{\mathrm{S}}_{5}\quad{\mathrm{S}}_{6}\quad0\quad0\quad0\quad0\right]}^{\mathrm{T}} \\ \mathbb{g}\left(\xi \right) = {\left[0\quad0 \quad L\left(1-{\xi }^{2}\right) 0\quad0\quad0 \quad L{\xi }^{2} \quad 0\quad{L}^{2}\left(\xi -{\xi }^{2}\right) 0\right]}^{\mathrm{T}} \\ \mathbb{h}\left(\xi \right) = {\left[0\quad0\quad0 \quad L\left(1-{\xi }^{2}\right) 0\quad0\quad0 \quad L{\xi }^{2} \quad 0\quad{L}^{2}\left(\xi -{\xi }^{2}\right)\right]}^{\mathrm{T}}\end{array}\right.$

(57)

According to the research content of Sections 3 and 4, the solution forms of elastic force and force Jacobi based on the García-Vallejo method and the Gerstmayr method can be obtained by reconstructing the $\mathbb{H}$ (as given in Appendix 4).

6. Numerical example

In this section, some numerical examples are used to verify the correctness of the proposed model.

6.1. Statics problem

A cantilever beam with a circular cross-section is shown in . Geometry properties of the beam are: the length is $L = 1m$ , the radius of the circular section is $r = 0.01 \; m$ . And its material properties are: Young's modulus E = 72 GPa, the Poisson's ratio is $\nu = 0.33$ .

Figure 2. Flexible pendulum.

DownLoad: Full-Size Img PowerPoint

Under the externally applied force $F = 10N$ , the displacement of the cantilever beam endpoint can be expressed as,

$\Delta y = \frac{F{l}^{3}}{3EI}$

(58)

The theoretical value of $\Delta y$ is 5.8946 mm, the beam element is simulated by Method Ⅰ (Modified García-Vallejo Method), Method Ⅱ (Modified Gerstmayr Method), Method Ⅲ (Modified Higher-Order Element), and the calculation values are shown in Table 1.

Table 1. The calculation values of

$F = 10N$ .

The theoretical value (mm)	Method Ⅲ (mm)	Method Ⅱ (mm)	Method Ⅰ (mm)
5.8946	5.7921	4.1328	4.1331

| Show Table

DownLoad: CSV

According to the data in , when $\nu = 0.33$ , the Poisson locking leads to the wrong results from Methods Ⅰ and Ⅱ, while Method Ⅲ is consistent with the theoretical results.

6.2. Flexible pendulum

Considering the motion trajectory of the right endpoint B of the flexible pendulum with circular cross-section under the action of gravity, the length of the pendulum is $2 \;m$ , as shown in . The beam element section is a circle with a radius of $0.02\; m$ , the elastic modulus of the material E = 720 MPa, the density is 7200 kg/m³, Poisson's ratio $\nu = 0$ and the gravitational acceleration is 9.81 m/s². In order to verify the correctness of the results, the same flexible pendulum model is established in MSC.ADAMS and the iteration steps length are set to 0.001 seconds. The simulation results are shown in Figures 3–6.

Figure 3. Flexible pendulum.

DownLoad: Full-Size Img PowerPoint

Figure 4. The trajectory of point B in the X direction.

DownLoad: Full-Size Img PowerPoint

Figure 5. The trajectory of point B in the Y direction.

DownLoad: Full-Size Img PowerPoint

Figure 6. The velocity trajectory of point B in the X direction.

DownLoad: Full-Size Img PowerPoint

Figures 4 and 5 show the trajectory of the flexible pendulum right endpoint B in the X direction and Y direction, it can be seen that the simulation results of MSC.ADAMS is almost consistent with the absolute nodal coordinate formulation, Method Ⅰ (Modified García-Vallejo Method), Method Ⅱ (Modified Gerstmayr Method), Method Ⅲ (Modified Higher-Order Element). Figures 6 and 7 show the velocity trajectory of the flexible pendulum right endpoint B in the X direction and Y direction. Compared with the velocity trajectory of MSC.ADAMS, the velocity curve obtained by the absolute nodal coordinate formulation shows that point B has a certain degree of jitter, indicating that the ANCF can describe the flexible characteristics of three-dimensional beams perfectly. The results of MSC. ADAMS demonstrated that Methods Ⅰ, Ⅱ and Ⅲ proposed in this paper is reliable.

Figure 7. The velocity trajectory of point B in the Y direction.

DownLoad: Full-Size Img PowerPoint

6.3. Spatial flexible double pendulum

shows the spatial flexible double pendulum with a circular cross-section and under the effect of gravity, O and A are ball hinges. The length of $\mathrm{O}\mathrm{A}$ and $\mathrm{A}\mathrm{B}$ are 0.5 m and 0.8 m, respectively. The beam element section is circular with a diameter of 0.02 m. The gravitational acceleration is 9.81 m/s². The material elastic modulus E = 300 MPa, Poisson's ratio $\nu = 0$ , the density is 7200 kg/m³, and the simulation time is 3 seconds. Figure 9 presents the displacement of the endpoint B of the spatial flexible double pendulum in the X-direction. Compared with the displacement results of Method Ⅰ and Method Ⅱ, Figure 8, Method Ⅲ shows the difference. As shown in Figure 10, the velocity curve of the Modified Higher-Order Element has obvious jitter, especially in the later stage, the results show that Method Ⅲ makes the beam more flexible.

Figure 8. The spatial flexible double pendulum.

DownLoad: Full-Size Img PowerPoint

Figure 9. The trajectory of point B in the X direction.

DownLoad: Full-Size Img PowerPoint

Figure 10. The velocity of point B in the X direction.

DownLoad: Full-Size Img PowerPoint

Run MATLAB on a laptop with an Intel Core 2.6GHz processor and 6GB of RAM, different numbers of elements were used to divide the beam, and the CPU time necessary to solve the problem as shown in Table 2.

Table 2. Computation time(s) spent by different methods.

Number of elements		OA (2); AB (5)	OA (4); AB (6)
Method Ⅰ	Step size	$1\times {10}^{-4}$	$1\times {10}^{-4}$
	Pre-processing	14.82	19.76
	Simulation	1371.2	1922.7
	Total	1386.0	1942.46
Method Ⅱ	Step size	$1\times {10}^{-4}$	$1\times {10}^{-4}$
	Pre-processing	9	11.98
	Simulation	573.74	906.59
	Total	582.74	918.57
Method Ⅲ	Step size	$1\times {10}^{-4}$	$1\times {10}^{-4}$
	Pre-processing	34.03	49.79
	Simulation	1210.5	1877.6
	Total	1244.53	1927.39

| Show Table

DownLoad: CSV

Table 1 shows the calculation time of the three methods. Method Ⅱ has better advantages, both the pre-processing speed and the simulation speed are greatly improved. While the pre-processing speed and the simulation speed of Method Ⅲ are significantly reduced due to higher-order interpolation polynomials and more nodal coordinates. What is worth mentioning is that, the invariant of the modified García-Vallejo method can be shared for elements with equal dimensions, it decreases significantly the amount of data to compute in the pre-processing stage. Figure 11 shows the spatial configuration of the flexible double pendulum.

Figure 11. The spatial falling flexible double pendulum at different time steps.

DownLoad: Full-Size Img PowerPoint

7. Summary and conclusions

In this paper, the shape function describing the absolute nodal coordinate formulation is divided into three parts, which are independent of the local coordinates in the floating coordinate system, the advantage of this is that the absolute nodal coordinate method can be applied to an arbitrary-section beam via the discrete form of the shape function. Moreover, considering the Euler-Bernoulli beam used in engineering structures is characterised by a symmetrical cross-section, small section size, zero odd integrals, and negligible high-order even integrals, hence the single-integral form reduces the difficulty of applying the absolute nodal coordinate method. Furthermore, single integrals are more efficient and improve pre-processing speed compared to volume integrals. Finally, the correctness of the proposed methods is verified through simulation examples. By comparing the computational costs of different methods, modified Gerstmayr Method obvious advantages in pre-processing speed and computational efficiency.

Conflict of interest

The authors declare that there are no conflicts of interest related to this paper.

Appendix 1

Discrete shape function,

$\overline{\boldsymbol{S}} = \boldsymbol{f}\left(\xi \right)+\eta \boldsymbol{g}\left(\xi \right)+\zeta \boldsymbol{h}\left(\xi \right)$

(A1)

The partial derivative with respect to local coordinates,

$\frac{\partial \overline{\boldsymbol{S}}}{\partial \pmb{\mathsf{α}}} = \left[{\overline{\boldsymbol{S}}, }_{x}{\overline{\boldsymbol{S}}, }_{y}{\overline{\boldsymbol{S}}, }_{z}\right]$

$\left\{\begin{array}{c}{\overline{\boldsymbol{S}}, }_{x} = \frac{\partial \overline{\boldsymbol{S}}}{\partial \xi }\frac{\partial \xi }{\partial x} = {\boldsymbol{f}}_{1}\left(\xi \right)+\eta {\boldsymbol{g}}_{1}+\zeta {\boldsymbol{h}}_{1}\\ {\overline{\boldsymbol{S}}, }_{y} = \frac{\partial \overline{\boldsymbol{S}}}{\partial \eta }\frac{\partial \eta }{\partial y} = {\boldsymbol{f}}_{2}\left(\xi \right) \\ {\overline{\boldsymbol{S}}, }_{y} = \frac{\partial \overline{\boldsymbol{S}}}{\partial \eta }\frac{\partial \eta }{\partial y} = {\boldsymbol{f}}_{2}\left(\xi \right) \end{array}\right.$

(A2)

here,

$\left\{\begin{array}{c}\begin{array}{c}{\boldsymbol{f}}_{1}\left(\xi \right) = {\left[-\frac{6}{l}\xi \left(1-\xi \right) \quad \left(1-4\xi +3{\xi }^{2}\right) \quad 0 \quad \quad 0 \quad \quad \frac{6}{l}\xi \left(1-\xi \right) \quad \left(-2\xi +3{\xi }^{2}\right) \quad 0 \quad 0\right]}^{\mathrm{T}}\\ {\boldsymbol{g}}_{1} = {\left[0 \quad 0 -1 \quad 0 \quad 0 \quad 0 \quad 1 \quad 0\right]}^{\mathrm{T}}\\ {\boldsymbol{h}}_{1} = {[0 \quad 0 \quad 0 -1 \quad 0 \quad 0 \quad 0 \quad 1]}^{\mathrm{T}} \end{array}\\ {\boldsymbol{f}}_{2}\left(\xi \right) = {\left[0 \quad 0 \quad \right(1-\xi \left) \quad 0 \quad 0 \quad 0 \quad \xi 0\right]}^{\mathrm{T}} \\ {\boldsymbol{f}}_{3}\left(\xi \right) = {\left[0 \quad 0 \quad 0 \quad \right(1-\xi \left) \quad 0 \quad 0 \quad 0 \quad \xi \right]}^{\mathrm{T}} \end{array}\right.$

(A3)

The discrete form of $\boldsymbol{H}$ ,

${\bf{H}} = {{\bf{w}}}_{\bf 1}+\mathrm{\eta }{{\bf{w}}}_{\bf 2}+\mathrm{\zeta }{{\bf{w}}}_{\bf 3}$

(A4)

The specific expression form of ${\boldsymbol{w}}_{\bf 1}, {\boldsymbol{w}}_{\bf 2}, {\boldsymbol{w}}_{\bf 3}$ are,

$\left\{\begin{array}{c}{\boldsymbol{w}}_{\bf 1} = \left[{\boldsymbol{f}}_{1}\left(\xi \right) \quad {\boldsymbol{f}}_{2}\left(\xi \right) \quad {\boldsymbol{f}}_{3}\left(\xi \right)\right]{\left(\frac{\partial {\boldsymbol{r}}_{{\bf{0}}}}{\partial \pmb{\mathsf{α}}}\right)}^{-1}\\ {\boldsymbol{w}}_{\bf 2} = \left[{\boldsymbol{g}}_{1} \quad {{\bf{O}}}_{1} \quad {{\bf{O}}}_{1}\right]{\left(\frac{\partial {\boldsymbol{r}}_{{\bf{0}}}}{\partial \pmb{\mathsf{α}}}\right)}^{-1} \quad \\ {\boldsymbol{w}}_{\bf 3} = \left[{\boldsymbol{h}}_{1} \quad {{\bf{O}}}_{1} \quad {{\bf{O}}}_{1}\right]{\left(\frac{\partial {\boldsymbol{r}}_{{\bf{0}}}}{\partial \pmb{\mathsf{α}}}\right)}^{-1} \end{array}\right.$

(A5)

and,

$\left\{\begin{array}{c}{\overline{\boldsymbol{H}}}_{1} = {\overline{w}}_{11}+\eta {\overline{w}}_{21}+\zeta {\overline{w}}_{31}\\ {\overline{\boldsymbol{H}}}_{2} = {\overline{w}}_{12}+\eta {\overline{w}}_{22}+\zeta {\overline{w}}_{32}\\ {\overline{\boldsymbol{H}}}_{3} = {\overline{w}}_{13}+\eta {\overline{w}}_{23}+\zeta {\overline{w}}_{33}\end{array}\right.$

(A6)

Appendix 2

$\mathrm{E}$ is the material modulus of elasticity, $\nu$ is the material Poisson's ratio, then, the lamer constant can be expressed as,

$\lambda = \frac{E\nu }{(1+\nu )(1-2\nu )}, \mu = \frac{E}{2(1+\nu )}$

1) Calculate ${K}_{1}$

${{\overline{\boldsymbol{H}}}_{1}}^{\boldsymbol{T}}{\overline{\boldsymbol{H}}}_{1} \\ = {{\overline{w}}_{11}}^{T}{\overline{w}}_{11}+\eta \left[{{\overline{w}}_{11}}^{T}{\overline{w}}_{21}+{{\overline{w}}_{21}}^{T}{\overline{w}}_{11}\right]+\zeta \left[{{\overline{w}}_{11}}^{T}{\overline{w}}_{31}+{{\overline{w}}_{31}}^{T}{\overline{w}}_{11}\right]\\+\eta \zeta \left[{{\overline{w}}_{21}}^{T}{\overline{w}}_{31}+{{\overline{w}}_{31}}^{T}{\overline{w}}_{21}\right]+{\eta }^{2}{{\overline{w}}_{21}}^{T}{\overline{w}}_{21}+{\zeta }^{2}{{\overline{w}}_{31}}^{T}{\overline{w}}_{31} \\ = {\boldsymbol{A}}_{1}+\eta {\boldsymbol{B}}_{1}+\zeta {\boldsymbol{C}}_{1}+\eta \zeta {\boldsymbol{D}}_{1}+{\eta }^{2}{\boldsymbol{E}}_{1}+{\zeta }^{2}{\boldsymbol{F}}_{1} \\ {{\overline{\boldsymbol{H}}}_{2}}^{\boldsymbol{T}}{\overline{\boldsymbol{H}}}_{2} \\ = {{\overline{w}}_{12}}^{T}{\overline{w}}_{12}+\eta \left[{{\overline{w}}_{12}}^{T}{\overline{w}}_{22}+{{\overline{w}}_{22}}^{T}{\overline{w}}_{12}\right]+\zeta \left[{{\overline{w}}_{12}}^{T}{\overline{w}}_{32}+{{\overline{w}}_{32}}^{T}{\overline{w}}_{12}\right]\\+\eta \zeta \left[{{\overline{w}}_{22}}^{T}{\overline{w}}_{32}+{{\overline{w}}_{32}}^{T}{\overline{w}}_{22}\right]+{\eta }^{2}{{\overline{w}}_{22}}^{T}{\overline{w}}_{22}+{\zeta }^{2}{{\overline{w}}_{32}}^{T}{\overline{w}}_{32} \\ = {\boldsymbol{A}}_{2}+\eta {\boldsymbol{B}}_{2}+\zeta {\boldsymbol{C}}_{2}+\eta \zeta {\boldsymbol{D}}_{2}+{\eta }^{2}{\boldsymbol{E}}_{2}+{\zeta }^{2}{\boldsymbol{F}}_{2} \\ {{\overline{\boldsymbol{H}}}_{3}}^{\boldsymbol{T}}{\overline{\boldsymbol{H}}}_{3} \\ = {{\overline{w}}_{13}}^{T}{\overline{w}}_{13}+\eta \left[{{\overline{w}}_{13}}^{T}{\overline{w}}_{23}+{{\overline{w}}_{23}}^{T}{\overline{w}}_{13}\right]+\zeta \left[{{\overline{w}}_{13}}^{T}{\overline{w}}_{33}+{{\overline{w}}_{33}}^{T}{\overline{w}}_{13}\right]\\+\eta \zeta \left[{{\overline{w}}_{23}}^{T}{\overline{w}}_{33}+{{\overline{w}}_{33}}^{T}{\overline{w}}_{23}\right]+{\eta }^{2}{{\overline{w}}_{23}}^{T}{\overline{w}}_{23}+{\zeta }^{2}{{\overline{w}}_{33}}^{T}{\overline{w}}_{33} \\ = {\boldsymbol{A}}_{3}+\eta {\boldsymbol{B}}_{3}+\zeta {\boldsymbol{C}}_{3}+\eta \zeta {\boldsymbol{D}}_{3}+{\eta }^{2}{\boldsymbol{E}}_{3}+{\zeta }^{2}{\boldsymbol{F}}_{3}$

(B1)

Substitute ${{\overline{\boldsymbol{H}}}_{\alpha }}^{\mathrm{T}}{\overline{\boldsymbol{H}}}_{\alpha }(\alpha = \mathrm{1, 2}, 3)$ into ${\boldsymbol{K}}_{1}$

${\boldsymbol{K}}_{1} = -\frac{3\mathtt{λ }+2\mathtt{μ }}{2}\sum\limits_{\alpha = 1}^{3}\int ({\boldsymbol{A}}_{\alpha }+\eta {\boldsymbol{B}}_{\alpha }+\zeta {\boldsymbol{C}}_{\alpha }+\eta \zeta {\boldsymbol{D}}_{\alpha }+{\eta }^{2}{\boldsymbol{E}}_{\alpha }+{\zeta }^{2}{\boldsymbol{F}}_{\alpha })\mathrm{d}V$

(B2)

2) Calculate ${K}_{2}\left(e\right)$

Let,

$\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_{k}^{i j} = \left(\overline{\boldsymbol{H}}_{l}^{\mathrm{T}} \overline{\boldsymbol{H}}_{m}\right)_{i, }^{\mathrm{T}}\left(\overline{\boldsymbol{H}}_{p}{ }^{\mathrm{T}} \overline{\boldsymbol{H}}_{q}\right)_{j}, \\ = \left(\boldsymbol{A}_{\mathrm{I}}+\eta \boldsymbol{B}_{\mathrm{I}}+\zeta \boldsymbol{C}_{\mathrm{I}}+\eta \zeta \boldsymbol{D}_{\mathrm{I}}+\eta^{2} \boldsymbol{E}_{\mathrm{I}}+\zeta^{2} \boldsymbol{F}_{\mathrm{I}}\right)_{i, }^{T}\left(\boldsymbol{A}_{\mathrm{II}}+\eta \boldsymbol{B}_{\mathrm{II}}+\zeta \boldsymbol{C}_{\mathrm{II}}+\eta \zeta \boldsymbol{D}_{\mathrm{II}}\right.\\ \left.+\eta^{2} \boldsymbol{E}_{\mathrm{II}}+\zeta^{2} \boldsymbol{F}_{\mathrm{II}}\right)_{j, }\\ = \boldsymbol{N}_{1}+\eta \boldsymbol{N}_{2}+\zeta \boldsymbol{N}_{3}+\eta \zeta \boldsymbol{N}_{4}+\eta^{2} \boldsymbol{N}_{5}+\zeta^{2} \boldsymbol{N}_{6}+\eta \zeta^{2} \boldsymbol{N}_{7}+\eta^{2} \zeta \boldsymbol{N}_{8}+\eta^{3} \boldsymbol{N}_{9}+\zeta^{3} \boldsymbol{N}_{10}\\ +\eta^{2} \zeta^{2} \boldsymbol{N}_{11}+\eta \zeta^{3} \boldsymbol{N}_{12}+\eta^{3} \zeta \boldsymbol{N}_{13}+\eta^{4} \boldsymbol{N}_{14}+\zeta^{4} \boldsymbol{N}_{15}$

(B3)

Here, $l, m, p, q = 1, 2, 3;\mathrm{N}\mathrm{o}\mathrm{t}\mathrm{e}:1.\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}\;l\ne m, p = m, q = l; 2.\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}\;l = m, p = q$ .

$\left\{\begin{array}{c}\begin{array}{c}{\boldsymbol{A}}_{\mathrm{{\rm I}}} = {{\overline{w}}_{1l}}^{\mathrm{T}}{\overline{w}}_{1m} \\ {\boldsymbol{B}}_{\mathrm{{\rm I}}} = \left[{{\overline{w}}_{1l}}^{\mathrm{T}}{\overline{w}}_{2m}+{{\overline{w}}_{2m}}^{\mathrm{T}}{\overline{w}}_{1l}\right]\\ {\boldsymbol{C}}_{\mathrm{{\rm I}}} = \left[{{\overline{w}}_{1l}}^{\mathrm{T}}{\overline{w}}_{3m}+{{\overline{w}}_{3m}}^{\mathrm{T}}{\overline{w}}_{1l}\right]\end{array}\\ \begin{array}{c}{\boldsymbol{D}}_{\mathrm{{\rm I}}} = \left[{{\overline{w}}_{2l}}^{\mathrm{T}}{\overline{w}}_{3m}+{{\overline{w}}_{3m}}^{\mathrm{T}}{\overline{w}}_{2l}\right]\\ {\boldsymbol{E}}_{\mathrm{{\rm I}}} = {{\overline{w}}_{2l}}^{\mathrm{T}}{\overline{w}}_{2m} \\ {\boldsymbol{F}}_{\mathrm{{\rm I}}} = {{\overline{w}}_{3l}}^{\mathrm{T}}{\overline{w}}_{3m} \end{array}\end{array}\right.$

$\left\{\begin{array}{c}\begin{array}{c}{\boldsymbol{A}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}} = {{\overline{w}}_{1p}}^{\mathrm{T}}{\overline{w}}_{1q} \\ {\boldsymbol{B}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}} = \left[{{\overline{w}}_{1p}}^{\mathrm{T}}{\overline{w}}_{2q}+{{\overline{w}}_{2q}}^{\mathrm{T}}{\overline{w}}_{1p}\right]\\ {\boldsymbol{C}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}} = \left[{{\overline{w}}_{1p}}^{\mathrm{T}}{\overline{w}}_{3q}+{{\overline{w}}_{3q}}^{\mathrm{T}}{\overline{w}}_{1p}\right]\end{array}\\ \begin{array}{c}{\boldsymbol{D}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}} = \left[{{\overline{w}}_{2p}}^{\mathrm{T}}{\overline{w}}_{3q}+{{\overline{w}}_{3q}}^{\mathrm{T}}{\overline{w}}_{2p}\right]\\ {\boldsymbol{E}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}} = {{\overline{w}}_{2p}}^{\mathrm{T}}{\overline{w}}_{2q} \\ {\boldsymbol{F}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}} = {{\overline{w}}_{3p}}^{\mathrm{T}}{\overline{w}}_{3q} \end{array}\end{array}\right.$

(B4)

and,

${\boldsymbol{N}}_{1} = {{A}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}{A}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, } \\ {\boldsymbol{N}}_{2} = \left[{{A}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{B}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\right]{\left[{B}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{A}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\right]}^{\mathrm{T}} \\ {\boldsymbol{N}}_{3} = \left[{{A}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{C}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\right]{\left[{C}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{A}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\right]}^{\mathrm{T}} \\ {\boldsymbol{N}}_{4} = \left[{{A}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{B}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{C}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{D}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\right]{\left[{D}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{C}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{B}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{A}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\right]}^{\mathrm{T}} \\ {\boldsymbol{N}}_{5} = \left[{{A}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{B}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{E}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\right]{\left[{E}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{B}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{A}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\right]}^{\mathrm{T}} \\ {\boldsymbol{N}}_{6} = \left[{{A}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{C}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{F}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\right]{\left[{F}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{C}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{A}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\right]}^{\mathrm{T}} \\ {\boldsymbol{N}}_{7} = \left[{{B}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{C}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{D}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{F}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\right]{\left[{F}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{D}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{C}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{B}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\right]}^{\mathrm{T}} \\ {\boldsymbol{N}}_{8} = \left[{{B}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{C}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{D}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{E}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\right]{\left[{E}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{D}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{C}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{B}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\right]}^{\mathrm{T}} \\ {\boldsymbol{N}}_{9} = \left[{{B}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{E}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\right]{\left[{E}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{B}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\right]}^{\mathrm{T}} \\ {\boldsymbol{N}}_{10} = \left[{{C}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{F}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\right]{\left[{F}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{C}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\right]}^{\mathrm{T}} \\ {\boldsymbol{N}}_{11} = \left[{{E}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{D}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{F}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\right]{\left[{F}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{D}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{E}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\right]}^{\mathrm{T}} \\ {\boldsymbol{N}}_{12} = \left[{{D}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{F}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\right]{\left[{F}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{D}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\right]}^{\mathrm{T}} \\ {\boldsymbol{N}}_{13} = \left[{{D}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\quad{{E}_{\mathrm{{\rm I}}i, }}^{\mathrm{T}}\right]{\left[{E}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\quad{D}_{\mathrm{{\rm I}}\mathrm{{\rm I}}j, }\right]}^{\mathrm{T}} \\ {\boldsymbol{N}}_{14} = {{E}_{1i, }}^{T}{E}_{1j, } \\ {\boldsymbol{N}}_{15} = {{F}_{1i, }}^{T}{F}_{1j, }$

(B5)

The elasticity coefficient of the corresponding ${\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_{k}}^{ij}$ may be denoted as ${{\mathit{ \Pi}} }_{k}$ , and the invariant matrix can be expressed as,

${\boldsymbol{C}}_{{\boldsymbol{K}}_{2}}^{ij} = \sum\limits_{k = 1}^{15}\int {{\mathit{ \Pi}} }_{k}{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_{k}}^{ij}\mathrm{d}V \\ = \sum\limits_{k = 1}^{15}\int {{\mathit{ \Pi}} }_{k}\left[\begin{array}{c}{\boldsymbol{N}}_{1}^{k}+\eta {\boldsymbol{N}}_{2}^{k}+\zeta {\boldsymbol{N}}_{3}^{k}+\eta \zeta {\boldsymbol{N}}_{4}^{k}+{\eta }^{2}{\boldsymbol{N}}_{5}^{k}\\ +{\zeta }^{2}{\boldsymbol{N}}_{6}^{k}+\eta {\zeta }^{2}{\boldsymbol{N}}_{7}^{k}+{\eta }^{2}\zeta {\boldsymbol{N}}_{8}^{k}+{\eta }^{3}{\boldsymbol{N}}_{9}^{k}+{\zeta }^{3}{\boldsymbol{N}}_{10}^{k}\\ +{\eta }^{2}{\zeta }^{2}{\boldsymbol{N}}_{11}^{k}+\eta {\zeta }^{3}{\boldsymbol{N}}_{12}^{k}+{\eta }^{3}\zeta {\boldsymbol{N}}_{13}^{k}+{\eta }^{4}{\boldsymbol{N}}_{14}^{k}+{\zeta }^{4}{\boldsymbol{N}}_{15}^{k}\end{array}\right]\mathrm{d}V$

(B6)

and,

${\left({\boldsymbol{K}}_{2}\right(\boldsymbol{e}\left)\right)}_{ij} = {\boldsymbol{e}}^{\mathrm{T}}\left\{\sum\limits_{k = 1}^{15}\int {{\mathit{ \Pi}} }_{k}\left[\begin{array}{c}{\boldsymbol{N}}_{1}^{k}+\eta {\boldsymbol{N}}_{2}^{k}+\zeta {\boldsymbol{N}}_{3}^{k}+\eta \zeta {\boldsymbol{N}}_{4}^{k}+{\eta }^{2}{\boldsymbol{N}}_{5}^{k}\\ +{\zeta }^{2}{\boldsymbol{N}}_{6}^{k}+\eta {\zeta }^{2}{\boldsymbol{N}}_{7}^{k}+{\eta }^{2}\zeta {\boldsymbol{N}}_{8}^{k}+{\eta }^{3}{\boldsymbol{N}}_{9}^{k}+{\zeta }^{3}{\boldsymbol{N}}_{10}^{k}\\ +{\eta }^{2}{\zeta }^{2}{\boldsymbol{N}}_{11}^{k}+\eta {\zeta }^{3}{\boldsymbol{N}}_{12}^{k}+{\eta }^{3}\zeta {\boldsymbol{N}}_{13}^{k}+{\eta }^{4}{\boldsymbol{N}}_{14}^{k}+{\zeta }^{4}{\boldsymbol{N}}_{15}^{k}\end{array}\right]\mathrm{d}V\right\}\boldsymbol{e}$

(B7)

Appendix 3

1) Calculate ${H}_{tn}{H}_{sl}$

${H}_{tn}{H}_{sl} = ({{w}^{1}}_{tn}+\eta {{w}^{2}}_{tn}+\zeta {{w}^{3}}_{tn})$ (

${{w}^{1}}_{sl}+\eta {{w}^{2}}_{sl}+\zeta {{w}^{3}}_{sl}$ )

$= {{w}^{1}}_{tn}{{w}^{1}}_{sl}+\eta \left[{{w}^{1}}_{tn}{{w}^{2}}_{sl}+{{w}^{2}}_{tn}{{w}^{1}}_{sl}\right]+\zeta \left[{{w}^{1}}_{tn}{{w}^{3}}_{sl}+{{w}^{3}}_{tn}{{w}^{1}}_{sl}\right] \\ +\eta \zeta \left[{{w}^{2}}_{tn}{{w}^{3}}_{sl}+{{w}^{3}}_{tn}{{w}^{2}}_{sl}\right]+{\eta }^{2}{{w}^{2}}_{tn}{{w}^{2}}_{sl}+{\zeta }^{2}{{w}^{3}}_{tn}{{w}^{3}}_{sl} \\ = \mathbb{A}+\eta \mathbb{B}+\zeta \mathbb{C}+\eta \zeta \mathbb{D}+{\eta }^{2}\mathbb{E}+{\zeta }^{2}\mathbb{F}$

(C1)

2) Calculate ${H}_{fa}{H}_{vb}{H}_{tn}{H}_{sl}$

${H}_{fa}{H}_{vb}{H}_{tn}{H}_{sl} \\ = ({\mathbb{A}}_{\mathrm{I}}+\eta {\mathbb{B}}_{\mathrm{I}}+\zeta {\mathbb{C}}_{\mathrm{I}}+\eta \zeta {\mathbb{D}}_{\mathrm{I}}+{\eta }^{2}{\mathbb{E}}_{\mathrm{I}}+{\zeta }^{2}{\mathbb{F}}_{\mathrm{I}})({\mathbb{A}}_{\mathrm{I}\mathrm{I}}+\eta {\mathbb{B}}_{\mathrm{I}\mathrm{I}}+\zeta {\mathbb{C}}_{\mathrm{I}\mathrm{I}}+\eta \zeta {\mathbb{D}}_{\mathrm{I}\mathrm{I}}+{\eta }^{2}{\mathbb{E}}_{\mathrm{I}\mathrm{I}}+{\zeta }^{2}{\mathbb{F}}_{\mathrm{I}\mathrm{I}}) \\ = {\mathbb{N}}_{1}+\eta {\mathbb{N}}_{2}+\zeta {\mathbb{N}}_{3}+\eta \zeta {\mathbb{N}}_{4}+{\eta }^{2}{\mathbb{N}}_{5}+{\zeta }^{2}{\mathbb{N}}_{6}+\eta {\zeta }^{2}{\mathbb{N}}_{7}+{\eta }^{2}\zeta {\mathbb{N}}_{8}+{\eta }^{3}{\mathbb{N}}_{9}+{\zeta }^{3}{\mathbb{N}}_{10} \\ +{\eta }^{2}{\zeta }^{2}{\mathbb{N}}_{11}+\eta {\zeta }^{3}{\mathbb{N}}_{12}+{\eta }^{3}\zeta {\mathbb{N}}_{13}+{\eta }^{4}{\mathbb{N}}_{14}+{\zeta }^{4}{\mathbb{N}}_{15}$

(C2)

Here, $f, v, t, s = 1, \cdots, 8;a, b, n, l = \mathrm{1, 2}, 3$ .

$\left\{\begin{array}{c}\begin{array}{c}{\mathbb{A}}_{\mathrm{{\rm I}}} = {{w}^{1}}_{fa}{{w}^{1}}_{vb} \\ {\mathbb{B}}_{\mathrm{{\rm I}}} = \left[{{w}^{1}}_{fa}{{w}^{2}}_{vb}+{{w}^{2}}_{fa}{{w}^{1}}_{vb}\right] \\ {\mathbb{C}}_{\mathrm{{\rm I}}} = \left[{{w}^{1}}_{fa}{{w}^{3}}_{vb}+{{w}^{3}}_{fa}{{w}^{1}}_{vb}\right] \end{array}\\ \begin{array}{c}{\mathbb{D}}_{\mathrm{{\rm I}}} = \left[{{w}^{2}}_{fa}{{w}^{3}}_{vb}+{{w}^{3}}_{fa}{{w}^{2}}_{vb}\right] \\ {\mathbb{E}}_{\mathrm{{\rm I}}} = {{w}^{2}}_{fa}{{w}^{2}}_{vb} \\ {\mathbb{F}}_{\mathrm{{\rm I}}} = {{w}^{3}}_{fa}{{w}^{3}}_{vb} \end{array}\end{array}\right. \\ \left\{\begin{array}{c}\begin{array}{c}{\mathbb{A}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}} = {{w}^{1}}_{tn}{{w}^{1}}_{sl} \\ {\mathbb{B}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}} = \left[{{w}^{1}}_{tn}{{w}^{2}}_{sl}+{{w}^{2}}_{tn}{{w}^{1}}_{sl}\right] \\ {\mathbb{C}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}} = \left[{{w}^{1}}_{tn}{{w}^{3}}_{sl}+{{w}^{3}}_{tn}{{w}^{1}}_{sl}\right] \end{array}\\ \begin{array}{c}{\mathbb{D}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}} = \left[{{w}^{2}}_{tn}{{w}^{3}}_{sl}+{{w}^{3}}_{tn}{{w}^{2}}_{sl}\right] \\ {\mathbb{E}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}} = {{w}^{2}}_{tn}{{w}^{2}}_{sl} \\ {\mathbb{F}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}} = {{w}^{3}}_{tn}{{w}^{3}}_{sl} \end{array}\end{array}\right.$

(C3)

and,

${\mathbb{N}}_{1} = {\mathbb{A}}_{\mathrm{{\rm I}}}{\mathbb{A}}_{\mathrm{I}\mathrm{{\rm I}}} \\ {\mathbb{N}}_{2} = \left[{\mathbb{A}}_{\mathrm{{\rm I}}}\quad{\mathbb{B}}_{\mathrm{{\rm I}}}\right]{\left[{\mathbb{B}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{A}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\right]}^{\mathrm{T}} \\ {\mathbb{N}}_{3} = \left[{\mathbb{A}}_{\mathrm{{\rm I}}}\quad{\mathbb{C}}_{\mathrm{{\rm I}}}\right]{\left[{\mathbb{C}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{A}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\right]}^{\mathrm{T}} \\ {\mathbb{N}}_{4} = \left[{\mathbb{A}}_{\mathrm{{\rm I}}}\quad{\mathbb{B}}_{\mathrm{{\rm I}}}\quad{\mathbb{C}}_{\mathrm{{\rm I}}}\quad{\mathbb{D}}_{\mathrm{{\rm I}}}\right]{\left[{\mathbb{D}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{C}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{B}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{A}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\right]}^{\mathrm{T}} \\ {\mathbb{N}}_{5} = \left[{\mathbb{A}}_{\mathrm{{\rm I}}}\quad{\mathbb{B}}_{\mathrm{{\rm I}}}\quad{\mathbb{E}}_{\mathrm{{\rm I}}}\right]{\left[{\mathbb{E}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{B}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{A}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\right]}^{\mathrm{T}} \\ {\mathbb{N}}_{6} = \left[{\mathbb{A}}_{\mathrm{{\rm I}}}\quad{\mathbb{C}}_{\mathrm{{\rm I}}}\quad{\mathbb{F}}_{\mathrm{{\rm I}}}\right]{\left[{\mathbb{F}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{C}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{A}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\right]}^{\mathrm{T}} \\ {\mathbb{N}}_{7} = \left[{\mathbb{B}}_{\mathrm{{\rm I}}}\quad{\mathbb{C}}_{\mathrm{{\rm I}}}\quad{\mathbb{D}}_{\mathrm{{\rm I}}}\quad{\mathbb{F}}_{\mathrm{{\rm I}}}\right]{\left[{\mathbb{F}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{D}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{C}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{B}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\right]}^{\mathrm{T}} \\ {\mathbb{N}}_{8} = \left[{\mathbb{B}}_{\mathrm{{\rm I}}}\quad{\mathbb{C}}_{\mathrm{{\rm I}}}\quad{\mathbb{D}}_{\mathrm{{\rm I}}}\quad{\mathbb{E}}_{\mathrm{{\rm I}}}\right]{\left[{\mathbb{E}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{D}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{C}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{B}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\right]}^{\mathrm{T}} \\ {\mathbb{N}}_{9} = \left[{\mathbb{B}}_{\mathrm{{\rm I}}}\quad{\mathbb{E}}_{\mathrm{{\rm I}}}\right]{\left[{\mathbb{E}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{B}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\right]}^{\mathrm{T}} \\ {\mathbb{N}}_{10} = \left[{\mathbb{C}}_{\mathrm{{\rm I}}}\quad{\mathbb{F}}_{\mathrm{{\rm I}}}\right]{\left[{\mathbb{F}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{C}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\right]}^{\mathrm{T}} \\ {\mathbb{N}}_{11} = \left[{\mathbb{E}}_{\mathrm{{\rm I}}}\quad{\mathbb{D}}_{\mathrm{{\rm I}}}\quad{\mathbb{F}}_{\mathrm{{\rm I}}}\right]{\left[{\mathbb{F}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{D}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{E}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\right]}^{\mathrm{T}} \\ {\mathbb{N}}_{12} = \left[{\mathbb{D}}_{\mathrm{{\rm I}}}\quad{\mathbb{F}}_{\mathrm{{\rm I}}}\right]{\left[{\mathbb{F}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{D}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\right]}^{\mathrm{T}} \\ {\mathbb{N}}_{13} = \left[{\mathbb{D}}_{\mathrm{{\rm I}}}\quad{\mathbb{E}}_{\mathrm{{\rm I}}}\right]{\left[{\mathbb{E}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\quad{\mathbb{D}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}\right]}^{\mathrm{T}} \\ {\mathbb{N}}_{14} = {\mathbb{E}}_{\mathrm{{\rm I}}}{\mathbb{E}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}} \\ {\mathbb{N}}_{15} = {\mathbb{F}}_{\mathrm{{\rm I}}}{\mathbb{F}}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}$

(C4)

Appendix 4

$\overline{\mathbb{S}} = \mathbb{f}\left(\xi \right)+\eta \mathbb{g}\left(\xi \right)+\zeta \mathbb{h}\left(\xi \right)$

(D1)

The partial derivative with respect to local coordinates,

$\frac{\partial \overline{\mathbb{S}}}{\partial \pmb{\mathsf{α}}} = \left[{\overline{\mathbb{S}}, }_{x}{\overline{\mathbb{S}}, }_{y}{\overline{\mathbb{S}}, }_{z}\right] \\ \left\{\begin{array}{c}{\overline{\mathbb{S}}, }_{x} = \frac{\partial \overline{\mathbb{S}}}{\partial \xi }\frac{\partial \xi }{\partial x} = {\mathbb{f}}_{1}\left(\xi \right)+\eta {\mathbb{g}}_{1}\left(\xi \right)+\zeta {\mathbb{h}}_{1}\left(\xi \right)\\ {\overline{\mathbb{S}}, }_{y} = \frac{\partial \overline{\mathbb{S}}}{\partial \eta }\frac{\partial \eta }{\partial y} = {\mathbb{f}}_{2}\left(\xi \right) \\ {\overline{\mathbb{S}}, }_{z} = \frac{\partial \overline{\mathbb{S}}}{\partial \zeta }\frac{\partial \zeta }{\partial z} = {\mathbb{f}}_{3}\left(\xi \right)\end{array}\right.$

(D2)

here,

$\left\{\begin{array}{c}{\mathbb{f}}_{1}\left(\xi \right) = {\left[-\frac{6}{l}\xi \left(1-\xi \right)\left(1-4\xi +3{\xi }^{2}\right) \quad 0 \quad 0 \quad \frac{6}{l}\xi \left(1-\xi \right)\left(-2\xi +3{\xi }^{2}\right) \quad 0 \quad 0 \quad 0 \quad 0\right]}^{\mathrm{T}}\\ {\mathbb{g}}_{1}\left(\xi \right) = {\left[0 \quad 0 -2\xi \quad 0 \quad 0 \quad 0 \quad 2\xi \quad 0 \quad l\left(1-2\xi \right) \quad 0\right]}^{\mathrm{T}}\\ {\mathbb{h}}_{1}\left(\xi \right) = {\left[0 \quad 0 \quad 0 \quad -2\xi \quad 0 \quad 0 \quad 0 \quad 2\xi \quad 0 \quad l\left(1-2\xi \right)\right]}^{\mathrm{T}}\\ {\mathbb{f}}_{2}\left(\xi \right) = {\left[0 \quad 0 \quad \left(1-{\xi }^{2}\right) \quad 0 \quad 0 \quad 0 \quad {\xi }^{2} \quad 0 \quad l\left(\xi -{\xi }^{2}\right)0\right]}^{\mathrm{T}}\\ {\mathbb{f}}_{3}\left(\xi \right) = {\left[0 \quad 0 \quad 0 \quad \left(1-{\xi }^{2}\right) \quad 0 \quad 0 \quad 0 \quad {\xi }^{2} \quad 0 \quad l\left(\xi -{\xi }^{2}\right)\right]}^{\mathrm{T}}\end{array}\right.$

(D3)

The discrete form of $\mathbb{H}$ ,

$\mathbb{H} = {\mathbb{w}}_{{{\bf{1}}}}+\eta {\mathbb{w}}_{{{\bf{2}}}}+\zeta {\mathbb{w}}_{{{\bf{3}}}}$

(D4)

The specific expression form of ${\mathbb{w}}_{1}, {\mathbb{w}}_{2}, {\mathbb{w}}_{3}$ are

$\left\{\begin{array}{c}{\mathbb{w}}_{\bf 1} = \left[{\mathbb{f}}_{1}\left(\xi \right) \quad {\mathbb{f}}_{2}\left(\xi \right) \quad {\mathbb{f}}_{3}\left(\xi \right)\right]{\left(\frac{\partial {\mathbb{r}}_{{\bf{0}}}}{\partial \pmb{\mathsf{α}}}\right)}^{-1}\\ {\mathbb{w}}_{\bf 2} = \left[{\mathbb{g}}_{1} \quad {\mathbb{O}}_{1} \quad {\mathbb{O}}_{1}\right]{\left(\frac{\partial {\mathbb{r}}_{{\bf{0}}}}{\partial \pmb{\mathsf{α}}}\right)}^{-1} \\ {\mathbb{w}}_{\bf 3} = \left[{\mathbb{h}}_{1} \quad {\mathbb{O}}_{1} \quad {\mathbb{O}}_{1}\right]{\left(\frac{\partial {\mathbb{r}}_{{\bf{0}}}}{\partial \pmb{\mathsf{α}}}\right)}^{-1} \end{array}\right.$

(D5)

and,

$\left\{\begin{array}{c}{\overline{\mathbb{H}}}_{1} = {\overline{w}}_{11}+\eta {\overline{w}}_{21}+\zeta {\overline{w}}_{31}\\ {\overline{\mathbb{H}}}_{2} = {\overline{w}}_{12}+\eta {\overline{w}}_{22}+\zeta {\overline{w}}_{32}\\ {\overline{\mathbb{H}}}_{3} = {\overline{w}}_{13}+\eta {\overline{w}}_{23}+\zeta {\overline{w}}_{33}\end{array}\right.$

(D6)

References

[1]	B. Y. Duan, The state-of-the-art and development trend of large space-borne deployable antenna, Electro-Mech. Eng., 33 (2017), 1-14.
[2]	J. Wu, S. Yan, L. Xie, G. Peng, Reliability apportionment approach for spacecraft solar array using fuzzy reasoning petri net and fuzzy comprehensive evaluation, Acta Astronaut., 76 (2012), 136-144. doi: 10.1016/j.actaastro.2012.02.023
[3]	X. M. Gao, D. P. Jin, T. Chen, Nonlinear analysis and experimental investigation of a rigid-flexible antenna system, Meccanica, 53 (2018), 33-48. doi: 10.1007/s11012-017-0708-z
[4]	Y. H. Zheng, P. Ruan, S. Cao, Deployable structure design and analysis for space membrane diffractive telescope, Hongwai Yu Jiguang Gongcheng/Infrared Laser Eng., 45 (2016), 133-137.
[5]	P. Li, Q. W. Ma, Y. P. Song, C. Liu, Q. Tian, S. P. Ma, et al., Deployment dynamics simulation and ground test of a large space hoop truss antenna reflector, Sci. Sin. Phys., Mech. Astron., 47 (2017), 104602.
[6]	X. L. Du, J. L. Du, H. Bao, G. H. Sun, Deployment analysis of deployable antennas considering cable net and truss flexibility, Aerosp. Sci. Technol., 82-83 (2018), 557-565.
[7]	K. N. Urata, J. Sumantyo, N. Imura, K. Ito, S. Gao, Development of a circularly polarized L-band SAR deployable mesh reflector antenna for microsatellite earth observation, in 2017 IEEE International Symposium on Antennas and Propagation & USNC/URSI National Radio Science Meeting, IEEE, (2017).
[8]	A. A. Shabana, Flexible multibody dynamics: Review of past and recent developments, Multibody Syst. Dyn., 1 (1997), 189-222. doi: 10.1023/A:1009773505418
[9]	D. A. Turcic, A. Midha, J. R. Bosnik, Dynamic analysis of elastic mechanism systems. Part Ⅱ: Experimental results, J. Dyn. Syst. Meas. Control, 106 (1984), 255.
[10]	M. Manish, Kineto-elastodynamic analysis of Watt's mechanism using ANSYS, Procedia Technol., 23 (2016), 51-59. doi: 10.1016/j.protcy.2016.03.073
[11]	P. W. Likins, Finite element appendage equations for hybrid coordinate dynamic analysis, Int. J. Solids Struct., 8 (1972), 709-731. doi: 10.1016/0020-7683(72)90038-8
[12]	J. L. Escalona, H. A. Hussien, A. A. Shabana, Application of the absolute nodal co-ordinate formulation to multibody system dynamics, J. Sound Vib., 214 (1998), 833-851. doi: 10.1006/jsvi.1998.1563
[13]	A. A. Shabana, An absolute nodal coordinates formulation for the large rotation and deformation analysis of flexible bodies, No. MBS96-1-UIC, University of Illinois at Chicago, 1996.
[14]	A. A. Shabana, C. J. Desai, E. Grossi, M. Patel, Generalization of the strain-split method and evaluation of the nonlinear ANCF finite elements, Acta Mech., 231 (2020), 1365-1376. doi: 10.1007/s00707-019-02558-w
[15]	Y. Zhang, C. Wei, Y. Zhao, C. L. Tan, Y. J. Liu, Adaptive ANCF method and its application in planar flexible cables, Acta Mech. Sin., 34 (2018), 199-213. doi: 10.1007/s10409-017-0721-4
[16]	C. H. Zhao, Dynamic model of flexible beam with arbitrary rigid motion based on floating frame of reference formulation, J. Shanghai Univ. Eng., 27 (2013), 39-42.
[17]	A. A. Shabana, Geometrically accurate infinitesimal-rotation spatial finite elements, Proc. Inst. Mech. Eng., 233 (2019), 182-187. doi: 10.1177/0954405417712550
[18]	J. L. Sun, Q. Tian, H. Y. Hu, Advances in dynamic modeling and optimization of flexible multibody systems, Chin. J. Theor. Appl. Mech., 51 (2019), 1565-1586.
[19]	M. Dibold, J. Gerstmayr, H. Irschik, A detailed comparison of the absolute nodal coordinate and the floating frame of reference formulation in deformable multibody systems, J. Comput. Nonlinear Dyn., 4 (2009).
[20]	R. Y. Yakoub, A New Three-Dimensional Absolute Coordinate Based Beam Element with Application to Wheel/Rail Interaction, University of Illinois at Chicago, 2001.
[21]	P. Li, C. Liu, Q. Tian, H. Y. Hu, Dynamics of a deployable mesh reflector of satellite antenna: Form-finding and modal analysis, J. Comput. Nonlinear Dyn., 11 (2016), 041017.
[22]	Q. Tian, J. Zhao, C. L. Liu, C. Y. Zhou, Dynamics of space deployable structures, in ASME 2015 international design engineering technical conferences and computers and information in engineering conference, 2015.
[23]	P. Lan, Q. L. Tian, Z. Q. Yu, A new absolute nodal coordinate formulation beam element with multilayer circular cross section, Acta Mech. Sin., 36 (2020), 82-96. doi: 10.1007/s10409-019-00897-4
[24]	Z. Wang, Q. Tian, H. Y. Hu, Dynamics of flexible multibody systems with hybrid uncertain parameters, Mech. Mach. Theory, 121 (2018), 128-147. doi: 10.1016/j.mechmachtheory.2017.09.024
[25]	M. Campanelli, M. Berzeri, A. A. Shabana, Performance of the incremental and non-incremental finite element formulations in flexible multibody problems, J. Mech. Design, 122 (2000), 498-507. doi: 10.1115/1.1289636
[26]	M. A. Omar, A. A. Shabana, A two-dimensional shear deformation beam for large rotation and deformation, J. Sound Vib., 243 (2001), 565-576. doi: 10.1006/jsvi.2000.3416
[27]	D. García-Vallejo, A. M. Mikkola, J. L. Escalona, A new locking-free shear deformable finite element based on absolute nodal coordinates, Nonlinear Dyn., 50 (2007), 249-264. doi: 10.1007/s11071-006-9155-4
[28]	A. A. Shabana, R. Y. Yakoub, Three-dimensional absolute nodal coordinate formulation for beam elements: theory, J. Mech. Design, 123 (2001), 614-621. doi: 10.1115/1.1410099
[29]	D. García-Vallejo, J. Mayo, J. L. Escalona, J. Domínguez, Efficient evaluation of the elastic forces and the Jacobian in the absolute nodal coordinate formulation, Nonlinear Dyn., 35 (2004), 313-329. doi: 10.1023/B:NODY.0000027747.41604.20
[30]	J. Gerstmayr, A. A. Shabana, A. Strain, Efficient integration of the elastic forces and thin three-dimensional beam elements in the absolute nodal coordinate formulation, Multibody Dyn., 2005.
[31]	C. Liu, Q. Tian, H. Y. Hu, Efficient computational method for dynamics of flexible multibody systems based on absolute nodal coordinate, Lixue Xuebao/Chin. J. Theor. Appl. Mech., 42 (2010), 1197-1205.
[32]	K. Otsuka, K. Makihara, ANCF-ICE beam element for modeling highly flexible and deployable aerospace structures, in AIAA Scitech 2019 Forum, 2019.
[33]	G. Orzechowski, Analysis of beam elements of circular cross section using the absolute nodal coordinate formulation, Arch. Mech. Eng., 59 (2012), 283-296. doi: 10.2478/v10180-012-0014-1
[34]	J. Gerstmayr, A. A. Shabana, Analysis of thin beams and cables using the absolute nodal co-ordinate formulation, Nonlinear Dyn., 45 (2006), 109-130. doi: 10.1007/s11071-006-1856-1

This article has been cited by:

1.	Ahmed A. Shabana, An overview of the ANCF approach, justifications for its use, implementation issues, and future research directions, 2023, 1384-5640, 10.1007/s11044-023-09890-z
2.	Ahmed A. Shabana, Mechanically induced temperature oscillations in the coupled thermo-elasticity analysis of articulated dynamical systems, 2024, 238, 1464-4193, 38, 10.1177/14644193231197262

Reader Comments

Your name:*

Email:*
© 2021 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

Mathematical Biosciences and Engineering

3.9

Metrics

Article views(2810) PDF downloads(98) Cited by(2)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(11) / Tables(2)

Mathematical Biosciences and Engineering

Application of discrete shape function in absolute nodal coordinate formulation

Related Papers:

Abstract

1. Introduction

2. Three-dimensional beam elements based on the ANCF

2.1. Absolute nodal coordinate formulation

2.2. Discrete shape function

3. Application of discrete shape function to invariant matrix method

4. Application of discrete shape function to Gerstmayr method

5. Poisson locking

6. Numerical example

6.1. Statics problem

6.2. Flexible pendulum

6.3. Spatial flexible double pendulum

7. Summary and conclusions

Conflict of interest

Appendix 1

Appendix 2

Appendix 3

Appendix 4

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Catalog

Mathematical Biosciences and Engineering

Application of discrete shape function in absolute nodal coordinate formulation

Related Papers:

Abstract

1. Introduction

2. Three-dimensional beam elements based on the ANCF

2.1. Absolute nodal coordinate formulation

2.2. Discrete shape function

3. Application of discrete shape function to invariant matrix method

4. Application of discrete shape function to Gerstmayr method

5. Poisson locking

6. Numerical example

6.1. Statics problem

6.2. Flexible pendulum

6.3. Spatial flexible double pendulum

7. Summary and conclusions

Conflict of interest

Appendix 1

Appendix 2

Appendix 3

Appendix 4

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog