# Vibrations of complex shells of revolution using B-spline finite elements

Ayech Benjeddou

Structural Mechanics and Coupled Systems Laboratory, CNAM, 2 rue ConteÂ, F-75003 Paris, France

This paper presents and evaluates a local-global B-spline finite element method for the modal analysis of joined thin shells of revolution. Each shell component is modelled by a super element. For smooth meridian, local displacements of the shell middle surface and their first derivatives are matched at the junction. However, for discontinuous shells, the $C^1$ continuity of tangential displacements is relaxed and a local to global frame transformation is applied for normal and longitudinal displacements. The super element formulation in local-global frame reference and the minimum number of degrees of freedom (three nodal displacements augmented by the bending rotation at the joints only), provided by the uniform cubic B-spline interpolation, are the main originalities of the present work. Efficiency and accuracy of the proposed technique are evaluated through various well known two-branch shell combinations. Results are compared to analytical, experimental and numerical solutions found in the literature. A very satisfactory behaviour of the approach was observed.

Keywords: Modal analysis; Junctions; Shells of revolution; B-spline finite elements; Local-global frame formulation

## 1. Introduction

Joined shells of revolution find wide use in pressure vessels, hermetic systems, nuclear reactors, space vehicles, aircrafts and rockets. They consist of thin-walled structures composed of two or more simple components having one common axis of symmetry and presenting a sudden change in any relevant shell parameter across the joint. Here, we are interested in geometric discontinuities due to (i) a discontinuity of the radius of the meridian, as in an hemisphere-cylinder combination (Fig. 1b), (ii) a discontinuity of the colatitude angle (Fig. 1c) or, equivalently, the second principal radius of curvature, as in a cone-cylinder combination, and (iii) an abrupt change of the thick-

ness. Major difficulties are generally met with the second type (ii).

Early theoretical and experimental investigations indicate that during the modal vibrations of joined shells, localized high bending stresses generally develop within a narrow region enclosing the joint. The vibratory nature of this phenomenon may lead to fatigue and thus requires some attention in the analysis and design of joined, particularly branched, shells. Indeed, the proper matching of boundary conditions between the shell components offers a challenging problem when these structures are of concern. Hence, in the following, a review of different methods used for the vibration analysis of shell combinations is made with a special emphasis on the continuity conditions enforced at the joints.
<figure><img src="image_1.png" /><figcaption>Fig. 1. Axisymmetric element notations and meridian geometric discontinuities.</figcaption></figure>

Vibration investigations of joined axisymmetric shells were firstly based on analytical or semi-analytical methods, supported sometimes by experimental measurements. This approach limits naturally the field of investigation to particular configurations. Hence, by the end of the 50s, Saunders and Pasley [1] presented experimental and analytical results for a sphere-cone combination. The analytical procedure was based on an adaptation of the approximate Rayleigh inextensional theory. Only tangential and circumferential displacements were then made continuous across the junction. Nevertheless, they obtained an agreement between predicted and measured natural frequencies within 5% for the first five modes of vibration. This particular configuration was also considered by Kalnins [2] as an illustrative example of a bending theory which introduced a notion of transfer matrix between the state space vector of the eight fundamental shell variables at a current point and that at the boundary. The continuity of all relevant variables at the end nodes of the segments was ensured. The solution of the eigenvalue problem was obtained by a multisegment, direct numerical integration approach and for the particular sphere-cone combination considered, the lowest natural frequency for the circumferential number $n = 3$ was 9% from that of reference [1].

Hu and Raney [3] generalized this treatment of the joint by introducing a coupling matrix which is necessary only when the colatitude angle discontinuity exists at the joint. This matrix transfers across the joint the eight fundamental shell variables from local co-ordinate systems to a common global one. The vibrational characteristics of a cone-cylinder combination were then investigated analytically and experimentally. It is worthwhile to notice that the State Space Method (or Transfer Matrix Method), known nowadays in sub-structure synthesis methods, finds its origin in references [2,3].

Vibrations of combined cylindrical shells with or without circular plates were investigated by Takahashi

et al. [4,5]. The Lagrangian of the combined system was expressed by unknown boundary values alone and frequency equations and modal patterns of vibration were obtained by minimizing the Lagrangian. The shell continuity conditions were expressed by equating the axial, circumferential and radial displacements and the slope due to radial displacement.

Recently, Tavakoli and Singh [6,7] proposed a sub-structure synthesis method based on the State Space Method (SSM) for the eigensolution of axisymmetric joined/hermetic thin shell structures. They developed a library of simple shell transfer matrices, that they used later to analyse joined/hermetic shells, such as hermetic capsule and a refrigeration compressor shell. Later, they supported their method experimentally by the modal analysis of a hermetic can which is composed of a circular cylinder with two circular end plates [7]. Discrepancies were noticed, in Ref. [6,7], between finite element and SSM results for the first circumferential mode.

Let us recall that the analysis of shells of revolution have been carried out using finite difference energy method [8] as well as the finite element method. Usually, continuity conditions between adjacent elements can be ensured simply by letting the two elements have the same set of nodal Degrees of Freedom (DOF) at the common node. However, these must be expressed in a common and uniquely defined co-ordinate system. Almost exclusively, the global co-ordinate system is retained, as in the wide used conical axisymmetric thin shell element (see Ref. [9], for example). Recently, Özakça and Hinton [10] used global displacement-based Mindlin-Reissner formulation to develop a family of $n$-noded $C^0$ elements. For branched axisymmetric shells, they expressed all displacements in the global co-ordinate system. Therefore, unexplained discrepancies were also noticed between SSM [6,7] and FEM results for the first circumferential mode, when symmetry conditions were enforced at the shell apices.
General free vibration analysis of shells of revolution has also been made using two-node local frame curved B-spline finite elements [11,13]. Hence, Luah and Fan [11] presented an isoparametric C² B-spline element having 9 DOFs/node. These are the three displacements and their first and second derivatives. However, Benjeddou [12,13] proposed a simpler super parametric C¹ B-spline element having fewer DOFs. That is, only at the meridian end nodes that first derivatives are used as DOFs. For example, for a finite element model of N elements, Luah and Fan [11] should get 9(N + 1) DOFs, whereas, only 3(N + 3) will be obtained by the B-spline formulation [12,13]. 6N DOFs are then saved by the latter approach. Detailed numerical comparisons of both elements can be found in Ref. [13] for single branched shells only. In fact, neither of them can deal with geometric discontinuities due to their local formulations. Advantages of using B-spline functions for the finite element method were discussed elsewhere [12,14] and are beyond the scope of the present work which focuses on the local-global super element approach to deal with joined shells of revolution.

The objective of this paper is to present and evaluate a local-global displacement-based finite element approach for the vibrations of joined shells of revolution. This work could be seen as an extension of the recent B-spline Finite Element Method (BSFEM) [12], initially devoted to simple shells of revolution. Hence, the COREBS element [13] will be enhanced in order to handle physically discontinuous axisymmetrical shells. Section 2 recalls the B-spline finite element formulation, in a local frame reference, for single shells of revolution and details its generalization to joint shells. Applications of the present approach for the vibration of various well known shell combinations are presented in Section 3. Finally, some concluding remarks on modelling and computation of joined axisymmetric shells close the paper.

## 2. Local-global formulation for joined axisymmetric shells

The technique presented here to treat joined thin shells of revolution takes advantage of the local frame reference BSFEM formulation and the several features of the COREBS element presented and evaluated in [12,13] on single branched shells. A practical feature of our BSFEM, is the separation between geometry and behaviour. This is useful, for example, for the higher modes representation or to decrease the mesh size without any change of geometric data. An abrupt change in shell thickness or material constants is merely handled by assigning different values of these parameters for the super elements (or shell components) meeting at the point of discontinuity.

Tangents, for geometric C¹ continuity, are evaluated automatically and independently for each super element [14].

The joined thin shell is idealized by a number of super elements, each representing a simple shell component or a part of a simple shell, with abrupt change of any of its parameters. For smooth continuous meridians, these super elements are then assembled in the shell middle surface local co-ordinate system with full C¹ continuities at the joint. For physically discontinuous shells, the tangential displacements C¹ continuity is relaxed, and a local to global transformation is used for longitudinal and normal displacements. A modified version, COREBS_M, of the above element is then obtained. COREBS_M can also use special end elements without derivatives of the tangential displacements, thus reducing the number of DOF at super element end nodes from six to four.

In the following, the local frame reference formulation of the COREBS element is briefly described. It can be used for smooth continuous shell meridians. For discontinuous shells, an enhanced version, COREBS_M, of the above element is presented at the end of this section. Both elements are based on a classical thin shell theory of Kirchhoff-Love type.

### 2.1. Geometry and displacement interpolations

A super element representing a simple shell component of the joined shell, is geometrically defined as a set of spans constituting the geometric mesh (or master mesh). In the kth span, of reduced parameter $\mu \in [-1,1]$, the cylindrical co-ordinates of a current point are represented by a uniform cubic B-spline interpolation as follows [14],

$$r_k(\mu) = \langle F_{k4}(\mu) \rangle \{R_{k4}\}, \quad z_k(\mu) = \langle F_{k4}(\mu) \rangle \{Z_{k4}\} \quad (1)$$

where,

$$\langle F_{k4}(\mu) \rangle = \langle F_1(\mu)F_2(\mu)F_3(\mu)F_4(\mu) \rangle$$

$$\langle R_{k4} \rangle = \langle R_{k-1}R_kR_{k+1}R_{k+2} \rangle,$$

$$\langle Z_{k4} \rangle = \langle Z_{k-1}Z_kZ_{k+1}Z_{k+2} \rangle$$

and $F_i(\mu), i = 1, \dots, 4$, are the cubic B-spline functions given in Appendix A. ($R_k, Z_k$) represent the cylindrical co-ordinates of the kth vertex of the control polygon of a super element. They are computed separately as indicated in [13,14].

An optional secondary (or slave) mesh can be automatically created merely by dividing each geometric span into M sub-spans which are the actual finite elements. The geometry of the ith finite element (subspan) of reduced parameter $\xi \in [-1,1]$, is controlled by
that of the kth span through the following relation,

$$ \mu = \frac{1}{M}(\xi + 2i - M - 1), \quad i = 1, \dots, M \qquad (2) $$

Hence, no additional computations are needed to deduce the slave mesh from the master one. Moreover, this separation of geometry and behaviour can be used beneficially to better represent particular situations as high modes or point forces, since coarse master mesh is often sufficient to describe geometry but a refined slave mesh is usually needed for behaviour modelling. Note that for $M = 1$, the slave mesh coincides with the given master mesh.

In the ith element (Fig. 1a) local meridional u, circumferential v and normal w displacements are classically decomposed in Fourier series. Only θ-symmetric components are retained in the following, since the θ-anti symmetric amplitudes can be deduced just by changing to the opposite the sign of v-values [13]. Then, the same cubic B-spline representation is used for the θ-symmetric Fourier components of the nth order harmonic.

$$ \begin{aligned} u_n(\xi) &= \langle \mathbf{F}_{i4}(\xi) \rangle \{\mathbf{U}_{ni4}\}, & w_n(\xi) &= \langle \mathbf{F}_{i4}(\xi) \rangle \{\mathbf{W}_{ni4}\}, \\ v_n(\xi) &= \langle \mathbf{F}_{i4}(\xi) \rangle \{\mathbf{V}_{ni4}\} \end{aligned} \qquad (3) $$

Note that the variables $U_{ni}$, $W_{ni}$, $V_{ni}$, ..., etc, of the displacement control polygon have no physical significance. They are simply variables of convenience (non-nodal control DOF) that control physical displacements (nodal physical DOF).

As can be seen from the above representations Eqs. (1) and (3) and Eq. (2), the element, at its minimum capability, is isoparametric ($M = 1$) in addition to the optional number of B-spline functions of the slave mesh ($M > 1$).

## 2.2. B-spline finite element stiffness and mass matrices

Standard manipulation with respect to local non-nodal (control) co-ordinates, defined in Eq. (3), gives the element potential energy, for the nth order harmonic, in the form,

$$ U_{ne} = \frac{1}{2} \langle \mathbf{D}_{ne} \rangle [\mathbf{K}_{ne}] \{\mathbf{D}_{ne}\} \qquad (4) $$

where,

$$ \langle \mathbf{D}_{ne} \rangle = \langle U_{ni-1} W_{ni-1} V_{ni-1} U_{ni} W_{ni} V_{ni} U_{ni+1} W_{ni+1} \\ V_{ni+1} U_{ni+2} W_{ni+2} V_{ni+2} \rangle $$

is the element vector of control DOF and,

$$ [\mathbf{K}_{ne}] = k\pi \int_{-1}^{+1} [\mathbf{B}_n(\mu, \xi)]^T [\mathbf{H}] [\mathbf{B}_n(\mu, \xi)] r(\mu) J(\xi) \, d\xi $$

defines the (12 × 12) sized element stiffness matrix. In the above expression, k takes the values 2 for n = 0 and 1 otherwise. [$\mathbf{B}_n(\mu, \xi)$] and [H] indicate the strain and elasticity matrices, respectively. They are given in the appendix. $J = ds/d\xi$ is the determinant of the Jacobian.

Similarly, for a harmonic motion with circular frequency $\omega$, the element kinetic energy can be written as,

$$ V_e = \frac{1}{2} \omega^2 \langle \mathbf{D}_{ne} \rangle [\mathbf{M}_e] \{\mathbf{D}_{ne}\} \qquad (5) $$

where,

$$ [\mathbf{M}_e] = k\pi \int_{-1}^{+1} \rho_s [N(\xi)]^T [N(\xi)] r(\mu) J(\xi) \, d\xi $$

is the (12 × 12) sized element mass matrix. [$N(\xi)$], designates the (3 × 12) sized interpolation matrix deduced from Eq. (3). In the above expression, rotatory inertia were neglected, that is why the mass matrix is independent of the circumferential number n except the fact that k = 2 for n = 0 and 1 for n > 0. Stiffness and mass matrices are evaluated numerically.

## 2.3. Assembly of element matrices and eigenvalue problem

Element stiffness and mass matrices are first assembled up to the geometric span level i.e., slave mesh assembly. Then, the span matrices are assembled again up to the super element level i.e., master mesh assembly. At this stage, the non-nodal DOFs are replaced by nodal DOFs so that the potential and kinetic energies of a super element can be put in the form,

$$ \begin{aligned} U_{ns} &= \frac{1}{2} \langle \mathbf{d}_{ns} \rangle [\mathbf{K}_{ns}] \{\mathbf{d}_{ns}\}, \\ V_s &= \frac{1}{2} \omega^2 \langle \mathbf{d}_{ns} \rangle [\mathbf{M}_s] \{\mathbf{d}_{ns}\} \end{aligned} \qquad (6) $$

where,

$$ \langle \mathbf{d}_{ns} \rangle = \langle u_{n1} w_{n1} v_{n1} \frac{du}{ds_{n1}} \frac{dw}{ds_{n1}} \frac{dv}{ds_{n1}} \dots u_{ni} w_{ni} \frac{du}{ds_{nN}} \frac{dw}{ds_{nN}} \frac{dv}{ds_{nN}} \rangle_{v_{ni} \dots u_{nN} w_{nN} v_{nN}} $$

is the 3(N + 3) × 1 sized super element vector of nodal DOF defined from the non-nodal DOF vector {$\mathbf{D}_{ns}$} as,

$$ \{\mathbf{d}_{ns}\} = [\mathbf{T}_s]\{\mathbf{D}_{ns}\} \qquad (7) $$

where,
$$
\langle \boldsymbol{D}_{ns} \rangle = \langle U_{n0} W_{n0} V_{n0} U_{n1} W_{n1} V_{n1} \dots U_{ni} W_{ni} \\
\qquad V_{ni} \dots U_{nN} W_{nN} V_{nN} U_{nN+1} W_{nN+1} V_{nN+1} \rangle
$$

Note that the 6 DOF at fictive nodes 0 and N + 1 are eliminated due to the six derivative DOFs and that, there are only 3 DOF per internal node of the super element and 6 DOF at its end nodes. N is the number of nodes in a super element (master and slave nodes). In Eq. (6), the inverted version of the above relation is used. A special algorithm was developed to take advantage of the particular topology of the [T_s] matrix [12].

A third step assembly on super elements is made to
obtain the whole shell stiffness [K_n] and mass [M]
matrices. Hence, Eq. (6) becomes, for the n-th order
harmonic,

$$
U_n = \frac{1}{2} \langle \mathbf{d}_n \rangle [\mathbf{K}_n] \{\mathbf{d}_n\}, \quad V = \frac{1}{2} \omega^2 \langle \mathbf{d}_n \rangle [\mathbf{M}] \{\mathbf{d}_n\} \quad (8)
$$

where, {$\mathbf{d}_n$} is now the total vector of DOFs of the
whole structure for the harmonic order $n$.

It is worthwhile to precise that all assembly steps are
made in the shell middle surface local co-ordinate sys-
tem. First (slave) and second (master) steps respect C²
continuity between the elements and between spans;
whereas, in the third step, only C¹ continuity is
ensured. Thus, across the middle surface junction of a
smooth joined shell the local displacements and their
first derivatives are continuous.

From the potential and kinetic energies, the total
energy functional of the joined structure $H_n = U_n - V$
is made stationary,

$$
\delta H_n = 0, \quad \forall \{\delta \mathbf{d}_n\} \tag{9}
$$

leading to the following generalized symmetric eigen-
value problem,

$$
([\mathbf{K}_n] - \omega^2[\mathbf{M}])\{\mathbf{d}_n\} = \{\mathbf{0}\} \quad (10)
$$

which is solved, for each harmonic order *n*, by a sub-
spaces iterative method.

2.4. Junction modelling

For colatitude angle φ discontinuity, as in a cylin-
der-cone or cylinder-plate combinations, the local co-
ordinate system is not unique. Continuity conditions at
the joint should be written in a global co-ordinate sys-
tem or as linear constraints in local frame systems. The
former solution was preferred because it is easier to
implement without increasing the system size. Since v
and dw/ds are also global DOF, only u and w DOF
should be transformed to radial ū and axial w̄ displace-
ments for each super element end node,

$$
u = \cos \varphi \bar{u} + \sin \varphi \bar{w}, \quad w = -\sin \varphi \bar{u} + \cos \varphi \bar{w} \quad (11)
$$

Prior to this transformation, the derivative DOF du/
ds and dv/ds are eliminated by making vanish the third
parametric derivatives of tangential displacements at
each super element end node. This leads to new interp-
olation functions for u and v displacements there.
Hence, for the super element first node,

$$
u_{n1}(\xi) = \langle N_{13}(\xi) \rangle \{U_{n13}\}, \tag{12}
$$

$$
v_{n1}(\xi) = \langle N_{13}(\xi) \rangle \{V_{n13}\}
$$

where,

$$
\langle N_{n13}(\xi) \rangle = \langle N_1^1(\xi) N_2^1(\xi) N_3^1(\xi) \rangle
$$

$$
\langle U_{13} \rangle = \langle U_1 U_2 U_3 \rangle, \quad \langle V_{13} \rangle = \langle V_1 V_2 V_3 \rangle
$$

and for the super element end node,

$$
\begin{equation}
\begin{aligned}
u_{nN}(\xi) &= \langle \mathbf{N}_{N3}(\xi) \rangle \{ \mathbf{U}_{nN3} \}, \\
v_{nN}(\xi) &= \langle \mathbf{N}_{N3}(\xi) \rangle \{ \mathbf{V}_{nN3} \}
\end{aligned}
\tag{13}
\end{equation}
$$

$$
v_{nN}(\xi) = \langle \mathbf{N}_{N3}(\xi) \rangle \{ \mathbf{V}_{nN3} \}
$$

where,

$$
\langle \mathbf{N}_{N3}(\xi) \rangle = \langle N_1^N(\xi) N_2^N(\xi) N_3^N(\xi) \rangle
$$

$$
\langle U_{N3} \rangle = \langle U_{N-2} U_{N-1} U_N \rangle, \quad \langle V_{N3} \rangle = \langle V_{N-2} V_{N-1} V_N \rangle
$$

The modified B-spline functions $N_i^j(i = 1,2,3; j = 1,N)$, detailed in Appendix A, are linear combinations of the original $F_i$ functions. This formulation eliminates the need for fictive nodes 0 and $N+1$ leading to lesser DOF at the end nodes (four instead of six).

Using the above tangential displacement interpolations, super element end finite elements have the following size reduced stiffness and mass matrices (10 × 10 sized instead of 12 × 12),

$$
\begin{align}
[\mathbf{K}_{ne}^j] &= k\pi \int_{-1}^{+1} [\mathbf{B}_n^j(\mu, \xi)]^\mathrm{T} [\mathbf{H}] [\mathbf{B}_n^j(\mu, \xi)] r(\mu) J(\xi) \, d\xi; \tag{14} \\
j &= 1, N \nonumber
\end{align}
$$

and,

$$
[\mathbf{M}_e^j] = k\pi \int_{-1}^{+1} \rho_s [N^j(\xi)]^\mathrm{T} [N^j(\xi)] r(\mu) J(\xi) \, d\xi \quad (15)
$$

where the new matrices [B^j] and [N^j] are given in
Appendix A.

After the elementary matrices assembling till the
super element level, the Eq. (7) is applied in the form,

$$
\{\tilde{\mathbf{d}}_{ns}\} = [\tilde{\mathbf{T}}_s]\{\tilde{\mathbf{D}}_{ns}\} \quad (16)
$$

where,
<figure><img src="image_2.png" /><figcaption>Fig. 2. Geometry and material properties of a hemisphere-cylinder combination ($E = 2.07 \ 10^{11}$ N/m$^2$, $\rho = 7800$ kg/m$^3$, $\nu = 0.3$).</figcaption></figure>

$$ \langle \tilde{d}_{ns} \rangle = \langle u_{n1} w_{n1} v_{n1} \frac{dw}{ds_{n1}} \dots u_{ni} w_{ni} v_{ni} \dots u_{nN} w_{nN} v_{nN} \frac{dw}{ds_{nN}} \rangle $$

is the super element nodal DOF vector, and

$$ \langle \tilde{D}_{ns} \rangle = \langle W_{n0} U_{n1} W_{n1} V_{n1} \dots U_{ni} W_{ni} V_{ni} \dots U_{nN} W_{nN} V_{nN} W_{nN+1} \rangle $$

the corresponding non-nodal control DOF vector.

In comparison to the above COREBS formulation, there is a gain of 4 DOF per super element for the present modified formulation. Moreover, COREBS_M super elements become compatible with standard axisymmetric finite elements such as those of Hermite type. The remaining steps follow the above COREBS

formulation till the following eigenvalue problem,

$$ ([\tilde{K}_n] - \omega^2[\tilde{M}])\{\tilde{d}_n\} = \{0\} \quad (17) $$

which is solved by a subspaces iterative method.

The present COREBS_M formulation can also be used in complete local formulation (without Eq. (11)) for continuous joined shells.

### 3. Application to vibrations of joined shells

The objective of this section is to show the validity of the above BSFEM approach to treat joined thin shells of revolution. Our results are compared to known analytical [6], experimental [7] and numerical results [10]. Possible explanations are also given to the deviations, mentioned in [6,7,10], between SSM and FEM solutions for the first order circumferential number ($n = 1$).

#### 3.1. Hermetic capsule

The modal characteristics of a hermetic capsule, which consists of a circular cylinder closed at both ends with hemispherical end-caps, was investigated by Tavakoli and Singh using a semi-analytical State Space Method [6], and Özakça and Hinton [10] by a 305 DOF $C^0$ Mindlin-Reissner axisymmetric shell finite element model. Only a half of the capsule was analyzed due to the structure mid-plane symmetry (Fig. 2).

To obtain the corresponding symmetric modes, the

<figure><img src="image_3.png" /><figcaption>Fig. 3. Natural frequencies of a hemisphere-cylinder combination.</figcaption></figure>
<table><caption>Table 1
Natural frequencies (Hz) of a hemisphere-cylinder combination</caption><thead><tr><th>n</th><th>m</th><th>Ref. [6]<br>SSM (ref.)</th><th>Without<br>B.C.<br>(Eq.<br>(18)) at the apex<br>[10]<br>C<sup>0</sup> FEM (305 DOF)</th><th>Err<sup>a</sup> (%)</th><th>COREBS<br>C<sup>1</sup> BSFEM (72 DOF)</th><th>Err<sup>a</sup> (%)</th><th>With<br>B.C.<br>(Eq.<br>(18)) at the apex<br>Ref. [10]<br>C<sup>0</sup> FEM (305 DOF)</th><th>Err<sup>a</sup> (%)</th><th>COREBS<br>C<sup>1</sup> BSFEM (72 DOF)</th><th>Err<sup>a</sup> (%)</th></tr></thead><tbody><tr><td>0</td><td>1</td><td>4011</td><td>4004</td><td>-0.17</td><td>4015</td><td>0.10</td><td>4004</td><td>-0.17</td><td>4015</td><td>0.10</td></tr><tr><td>2</td><td>2</td><td>6292</td><td>6292</td><td>0.</td><td>6313</td><td>0.33</td><td>6292</td><td>0.</td><td>6313</td><td>0.33</td></tr><tr><td>3</td><td>3</td><td>6812</td><td>6802</td><td>-0.15</td><td>6838</td><td>0.38</td><td>6802</td><td>-0.15</td><td>6838</td><td>0.38</td></tr><tr><td>1</td><td>1</td><td>2740</td><td>2740</td><td>0.</td><td>2746</td><td>0.22</td><td><b>1016</b></td><td><b>-62.9</b></td><td>2746</td><td>0.22</td></tr><tr><td>2</td><td>2</td><td>5790</td><td>5793</td><td>0.05</td><td>5804</td><td>0.24</td><td><b>3548</b></td><td><b>-38.7</b></td><td>5804</td><td>0.24</td></tr><tr><td>3</td><td>3</td><td>5933</td><td>5933</td><td>0.</td><td>5960</td><td>0.45</td><td><b>5803</b></td><td><b>-2.19</b></td><td>5960</td><td>0.45</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 4. Geometry and material properties of a cone-cylinder combination ($E = 2.07 \ 10^{11}$ N/m$^2$, $\rho = 7800$ kg/m$^3$, $v = 0.3$).</figcaption></figure>

boundary conditions $u = dw/ds = 0$ were enforced at the middle of the cylinder. At the apex of the hemisphere the following symmetry conditions (in local coordinates) are theoretically required,

$$u = v = \frac{dw}{ds} \quad \text{for } n = 0;$$

$$u + v = 0, \quad w = 0, \quad \frac{dw}{ds} + \frac{u}{R} = 0 \quad \text{for } n = 1$$

$$u = w = v = \frac{dw}{ds} = 0 \quad \text{for } n \ge 2 \qquad (18)$$

Due to numerical integrations, these relations could be ignored, but in no way to be applied incompletely. This can lead to biased errors as was observed in Ref. [10], where only radial displacement was fixed.

The hermetic capsule was analyzed using two super elements with 10 COREBS elements for each component of the compound shell. This gives a total of 72 DOF (to compare to the 305 DOF used in [10]). Fig. 3 presents the first three frequencies for the first four circumferential mode numbers, in reference to those of [6]. Good accuracy was obtained with a few number of elements as well as DOFs, showing good correlation between analytical (SSM) and numerical (BSFEM) results. Similar frequencies were obtained with and without applying symmetry conditions (18) at the apex (cf. Table 1). However, a slight reduction of about 5% of some mode amplitudes was observed when Eqs. (18) are enforced.

The complete $C^1$ continuity formulation (COREBS) behaves very well in this test example because the joined shell is smoothly continuous.

## 3.2. Cone-cylinder combination

Experimental and analytical results for a truncated cone-cylinder combination were first given by Hu and

<footnote><sup>a</sup> Err = 100(FEM-SSM)/SSM.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 5. Natural frequencies of a cone-cylinder combination.</figcaption></figure>

Raney [3]. Both freely supported and clamped edges were studied for two shell models of different geometrical data. A complete cone-cylinder combination, as that shown in Fig. 4, was then investigated analytically by Tavakoli and Singh using SSM [6], and by the FEM by Özakça and Hinton [10] with a 305 DOF cubic four-noded C⁰ Mindlin-Reissner finite element model.

The cone-cylinder combination, shown in Fig. 4 was analyzed by the present BSFEM formulation using two super elements with 10 COREBS_M type finite element (without special elements) for each branch to model the joined shell. The required symmetry conditions at the apex could be written, in the cylindrical co-ordinate system, as

$$ \tilde{u} = v = \frac{dw}{ds} = 0, \quad \text{for } n = 0 $$

$$ \tilde{u} + v = 0, \quad \tilde{w} = \frac{dw}{ds} = 0, \quad \text{for } n = 1 $$

$$ \tilde{u} = \tilde{w} = v = \frac{dw}{ds} = 0, \quad \text{for } n \ge 2 \qquad (19) $$

The base of the cylinder was also clamped i.e., $\tilde{u} = \tilde{w} = v = dw/ds = 0$.

As can be seen in Fig. 5 representing the first three natural frequencies of the joined shell, for $n = 0$ to 3, COREBS_M provides good results with a small number of DOFs.

Similar frequencies were obtained with or without applying the symmetry conditions at the cone apex (cf.

Table 2). Hence, COREBS_M represents correctly and implicitly the symmetry conditions (19). A possible explanation of the deviations between SSM and FEM results, observed in Refs. [6,10] for $n = 1$, is that the symmetry conditions (19) were incompletely applied. In fact only the radial displacement was fixed in [10].

The local formulation (COREBS) does not give good results for this test case because of the colatitude angle discontinuity between the cylinder and the cone. Local to global transformation was then needed to obtain correct results (COREBS_M).

### 3.3. Hermetic can

A hermetic can structure, composed of a circular cylinder with two circular end plates is now investigated by the BSFEM approach. Taking advantage of the structure mid-plane symmetry, only half of the shell combination is analyzed (Fig. 6). Symmetry conditions enforced at the mid-section of the cylinder are

<figure><img src="image_6.png" /><figcaption>Fig. 6. Geometry and material properties of a plate-cylinder combination ($E = 2.07 \ 10^{11}$ N/m², $\rho = 7800$ kg/m³, $v = 0.3$).</figcaption></figure>
<table><caption>Table 2
Natural frequencies (Hz) of a cone-cylinder combination</caption><thead><tr><th>n</th><th>m</th><th>Ref. [6]<br>SSM (ref.)</th><th>Without<br>B.C.<br>(Eq.<br>(19)) at the apex<br>Ref. [10]<br>C<sup>0</sup> FEM (305 DOF)</th><th>Err<sup>a</sup> (%)</th><th>COREBS_M<br>C<sup>1</sup> BSFEM (72 DOF)</th><th>Err<sup>a</sup> (%)</th><th>With<br>B.C.<br>(Eq.<br>(19)) at the apex<br>Ref. [10]<br>C<sup>0</sup> FEM (305 DOF)</th><th>Err<sup>a</sup> (%)</th><th>COREBS_M<br>C<sup>1</sup> BSFEM (72 DOF)</th><th>Err<sup>a</sup> (%)</th></tr></thead><tbody><tr><td>0</td><td>1</td><td>2525</td><td>2516</td><td>-0.36</td><td>2346</td><td>-7.09</td><td>2516</td><td>-0.36</td><td>2346</td><td>-7.09</td></tr><tr><td>2</td><td>2</td><td>5585</td><td>5577</td><td>-0.14</td><td>5253</td><td>-5.95</td><td>5577</td><td>-0.14</td><td>5253</td><td>-5.95</td></tr><tr><td>3</td><td>3</td><td>6470</td><td>6470</td><td>0</td><td>6393</td><td>-1.19</td><td>6470</td><td>0</td><td>6393</td><td>-1.19</td></tr><tr><td>1</td><td>1</td><td>1962</td><td>1962</td><td>0</td><td>1985</td><td>1.17</td><td><strong>2446</strong></td><td><strong>24.7</strong></td><td>1986</td><td>1.22</td></tr><tr><td>2</td><td>2</td><td>3432</td><td>3423</td><td>-0.26</td><td>3204</td><td>-6.64</td><td><strong>3817</strong></td><td><strong>11.2</strong></td><td>3207</td><td>-6.55</td></tr><tr><td>3</td><td>3</td><td>5349</td><td>5350</td><td>0.02</td><td>5345</td><td>-0.07</td><td>5582</td><td>4.35</td><td>5345</td><td>-0.07</td></tr></tbody></table>

those used in the first example. Only symmetric modes are computed.

The joined shell was modelled with two super elements containing nine COREBS_M elements for each branch. The first three natural frequencies for circular mode numbers $n = 0, 1$ compare favourably to experimental, SSM [7], and FEM [10] results (see Table 3).

Here, the discontinuity of the colatitude angle is very sharp, that is why the local to global transformation was not sufficient to obtain good results. Moreover, the tangents are discontinuous and thus $du/ds$ and $dv/ds$ end nodal DOF should be eliminated using COREBS_M. With these enhancements to the original full local formulation (COREBS), good results were reached with few elements as well as DOFs, as can be seen from Table 3.

Note that, as in the previous examples, our results are similar with or without application of the symmetry conditions at the centre of the circular end plate. Thus, our BSFEM formulation represents implicitly these conditions.

## 4. Concluding remarks

A local $C^1$ B-spline Finite Element formulation was presented and evaluated on joined shells of revolution. Each branch of the combined shell was modelled by a super element. Some modifications were made to our previous element, COREBS [12], to treat shell discontinuities. Hence, a local to global transformation of normal and longitudinal displacements was applied at the joint and the $C^1$ continuity was relaxed for the tangential components. An enhanced version of the above element was obtained (COREBS_M) which is compatible with standard axisymmetric finite elements, but with fewer DOF since it has only 3 DOF inside each super element and 4 DOF at its ends.

Validation of the present approach, through several benchmarks shows good accuracy and compares favourably to other numerical, analytical and experimental results. Moreover, only a few elements were needed to model the different test examples accurately. Application of the symmetry conditions at apices has no significant perturbations on the results in the contrary to other formulations, for which deviations were met between analytical and numerical results. It has also been shown that full local $C^1$ continuity BSFEM can handle accurately joined shells which do not present colatitude angle discontinuity. However, super elements should be used for mesh refinement even for smooth meridians. Therefore, non-uniform or rational B-splines could be used to overcome this limitation if it is so considered.

Some features of the present element were reported on its sandwich version [15] to study dynamic response

<footnote><sup>a</sup> Err = 100(FEM-SSM)/SSM.</footnote>
of damped sandwich shells of revolution. The local for-
mulation and the minimum number of rotational DOF
make this element very attractive for the study of the
vibrations of structures containing liquids or a gas.
Hence, the enhanced version of the sandwich element
is being extended to be used in hybrid active/passive
control of fluid-structure coupling systems.

Acknowledgements

The author would like to thank Prof. R. Ohayon for
helpful suggestions and comments on the preparation
of this paper.

**Appendix A**

A.1. Uniform cubic B-spline interpolation

Inside a super element, the interpolation matrix of
the symmetric Fourier components of the local displa-
cements u (longitudinal), w (normal) and v (circumfer-
ential) can be written as,

$$
\begin{equation}
\begin{split}
[N(\xi)] = \begin{bmatrix}
F_1(\xi) & 0 & 0 \\
0 & F_1(\xi) & 0 \\
0 & 0 & F_1(\xi)
\end{bmatrix}
& \begin{bmatrix}
F_2(\xi) & 0 & 0 \\
0 & F_2(\xi) & 0 \\
0 & 0 & F_2(\xi)
\end{bmatrix} \\
& \begin{bmatrix}
F_3(\xi) & 0 & 0 \\
0 & F_3(\xi) & 0 \\
0 & 0 & F_3(\xi)
\end{bmatrix} \\
& \begin{bmatrix}
F_4(\xi) & 0 & 0 \\
0 & F_4(\xi) & 0 \\
0 & 0 & F_4(\xi)
\end{bmatrix}
\end{split}
\tag{A1}
\end{equation}
$$

where,

$$
F_1(\xi) = \frac{1}{48} (-\xi^3 + 3\xi^2 - 3\xi + 1),
$$

$$
F_2(\xi) = \frac{1}{48} (3\xi^3 - 3\xi^2 - 15\xi + 23),
$$

$$
F_3(\xi) = F_2(-\xi), \quad F_4(\xi) = F_1(-\xi) \tag{A2}
$$

For the first element of a super element, the modified
interpolation matrix is,

$$
\begin{equation}
\begin{split}
[N^1(\xi)] = \begin{bmatrix}
0 & N_1^1(\xi) & 0 & 0 \\
F_1(\xi) & |0 & F_2(\xi) & 0 \\
0 & 0 & 0 & N_1^1(\xi)
\end{bmatrix} & \\
& \begin{bmatrix}
N_2^1(\xi) & 0 & 0 \\
|0 & F_3(\xi) & 0 \\
0 & 0 & N_2^1(\xi)
\end{bmatrix} \\
& \begin{bmatrix}
N_3^1(\xi) & 0 & 0 \\
|0 & F_4(\xi) & 0 \\
0 & 0 & N_3^1(\xi)
\end{bmatrix}
\end{split}
\tag{A3}
\end{equation}
$$

where the $N_i^1$, are linear combinations of the original
basis functions $F_i$,

$$
\begin{align}
N_1^1(\xi) &= F_2(\xi) + 3F_1(\xi), \notag \\
N_2^1(\xi) &= F_3(\xi) - 3F_1(\xi), \tag{A4} \\
N_3^1(\xi) &= F_4(\xi) + F_1(\xi) \notag
\end{align}
$$

Similar interpolation matrix is obtained for the Nth ele-
ment of a super element,

$$
\begin{equation}
\begin{split}
[N^N(\xi)] = \begin{bmatrix}
N_1^N(\xi) & 0 & 0 \\
0 & F_1(\xi) & 0 \\
0 & 0 & N_1^N(\xi)
\end{bmatrix} & \\
& \begin{bmatrix}
N_2^N(\xi) & 0 & 0 \\
0 & F_2(\xi) & 0 \\
0 & 0 & N_2^N(\xi)
\end{bmatrix} \\
& \begin{bmatrix}
N_3^N(\xi) & 0 & 0 & 0 \\
0 & F_3(\xi) & 0 & F_4(\xi) \\
0 & 0 & N_3^N(\xi) & 0
\end{bmatrix}
\end{split}
\tag{A5}
\end{equation}
$$

here, the modified B-spline functions can be deduced
from those of the first element by the following re-
lation,

$$
N_{4-i}^{N}(\xi) = N_{i}^{1}(-\xi), \quad i = 1, 2, 3 \qquad (A6)
$$

Note that, in Eqs. (A3) and (A5) only the interpola-
tions of the tangential displacements are modified.

A.2. Elementary strain matrices

For any element inside a super element, the strain
matrix of the n-th harmonic order can be expressed as,

$$
[\mathbf{B}_n(\mu, \xi)] = [\mathbf{L}_n(\mu, \xi)][\mathbf{N}(\xi)] \quad (A7)
$$

whereas, for the special end elements, the following re-
lations are used,

$$
[\mathbf{B}_n^j(\mu, \xi)] = [\mathbf{L}_n(\mu, \xi)][\mathbf{N}^j(\xi)] \quad j = 1, N \quad (A8)
$$
<table><caption>Table 3
Natural frequencies (Hz) of a plate-cylinder combination</caption><thead><tr><th>n</th><th>m</th><th>Ref. [7] Exp.</th><th>Ref. [7] SSM (ref.)</th><th>Ref. [10] C<sup>0</sup> FEM (275 DOF)</th><th>Err<sup>a</sup> (%)</th><th>COREBS_M C<sup>1</sup> BSFEM (66 DOF)</th><th>Err<sup>a</sup> (%)</th></tr></thead><tbody><tr><td>0</td><td>1</td><td>360</td><td>359</td><td>353</td><td>-1.67</td><td>365</td><td>1.67</td></tr><tr><td>2</td><td>-</td><td>-</td><td>1412</td><td>1384</td><td>-1.98</td><td>1426</td><td>0.99</td></tr><tr><td>3</td><td>-</td><td>-</td><td>3173</td><td>3104</td><td>-2.17</td><td>3195</td><td>0.69</td></tr><tr><td>1</td><td>1</td><td>620</td><td>783</td><td>733</td><td>-6.38</td><td>754</td><td>-3.70</td></tr><tr><td>2</td><td>-</td><td>-</td><td>2218</td><td>2157</td><td>-2.75</td><td>2187</td><td>-1.40</td></tr><tr><td>3</td><td>-</td><td>-</td><td>3799</td><td>3716</td><td>-2.18</td><td>3918</td><td>3.13</td></tr></tbody></table>

<sup>a</sup> Err = 100(FEM-SSM)/SSM.

where,

$$ [L_n(\mu, \xi)] = \begin{bmatrix} \frac{d}{ds} & -\frac{1}{R_s} & 0 \\ \frac{\cos \varphi}{r} & -\frac{\sin \varphi}{r} & \frac{n}{r} \\ -\frac{n}{r} & 0 & \frac{d}{ds} - \frac{\cos \varphi}{r} \\ -\frac{1}{R_s} \frac{d}{ds} + \frac{1}{R_s^2} \frac{dR_s}{ds} & -\frac{d^2}{ds^2} & 0 \\ -\frac{\cos \varphi}{rR_s} & -\frac{\cos \varphi}{r} \frac{d}{ds} + \frac{n^2}{r^2} & -\frac{n}{r^2} \sin \varphi \\ \frac{n}{rR_s} & \frac{2n}{r} \frac{d}{ds} - \frac{2n}{r^2} \cos \varphi & -\frac{\sin \varphi}{r} \frac{d}{ds} - \frac{\cos \varphi}{r} \left( \frac{1}{R_s} - 2 \frac{\sin \varphi}{r} \right) \end{bmatrix} \quad (A9) $$

is the derivative operator of the strain tensor.

### A.3. Elasticity matrix

For an isotropic elastic material, the elasticity matrix is given by,

$$ [H] = \frac{Eh}{1 - v^2} $$

$$ \begin{bmatrix} 1 & v & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \\ & \frac{1-v}{2} & 0 & 0 & 0 & 0 \\ & & \frac{h^2}{12} & \frac{h^2}{12} & 0 & 0 \\ \text{sym.} & & & \frac{h^2}{12} & 0 & 0 \\ & & & & \frac{h^2}{12} & \frac{1-v}{2} \end{bmatrix} \quad (A10) $$

where, *E*, is the Young's modulus, *v*, is the Poisson's ratio and *h*, the shell thickness.

### References

[1] Saunders H, Pasley PR. Inextensional vibrations of a sphere-cone shell combination. Journal of the Acoustical Society of America 1959;31:579-83.

[2] Kalnins A. Free vibrations of rotationally symmetric shells. Journal of the Acoustical Society of America 1964;36:1355-65.

[3] Hu WCL, Raney JP. Experimental and analytical study of vibrations of joined shells. AIAA Journal 1967;5:976-80.

[4] Takahashi S, Hirano Y, Suzuki K. Vibrations of a system connecting two cylindrical shells. Bulletin of the JSME 1970;13:850-7.

[5] Takahashi S, Hirano Y. Vibration of a combination of circular plates and cylindrical shells. Bulletin of the JSME 1971;14:20-8.

[6] Tavakoli MS, Singh R. Eigensolutions of joined/hermetic shell structures using the state space method. Journal of Sound and Vibration 1989;100:97-123.

[7] Tavakoli MS, Singh R. Modal analysis of a hermetic can. Journal of Sound and Vibration 1990;136:141-5.

[8] Bushnell D. Stress, stability and vibration of complex,
branched shells of revolution. *Computers and Structures* 1974;4:399–435.

[9] Ohayon R, Nicolas-Vullierme B. An efficient shell finite element for the computation of the vibrations of fluid-structure systems of revolution. SMIRT 6, TP No. 1981-76, 1981.

[10] Özakça M, Hinton E. Free vibration analysis and optimization of axisymmetric plates and shells-I. finite element formulation. *Computers and Structures* 1994;52:1181–97.

[11] Luah MH, Fan SC. General free vibration analysis of shells of revolution using the spline finite element method. *Computers and Structures* 1989;33:1153–62.

[12] Benjeddou A. A B-spline finite element method for the analysis of homogeneous and sandwich shells with elastic

or viscoelastic core (in French). PhD thesis. Compiegne University of Technology, France, 1994.

[13] Benjeddou A, Hamdi MA. A new B-spline finite element for the dynamic analysis of shells of revolution. *European Journal of Finite Elements* 1994;3:101–26 in French.

[14] Guerich M, Benjeddou A, Hamdi MA. Uniform B₃-spline interpolation: A finite element approach. In: *Proceedings of the 4th Maghrebian Colloquium on Numurical Models for Engineers, Algiers, Algeria, vol. 2*. 1993. p. 724–9 (in French).

[15] Benjeddou A, Hamdi MA. A B-spline finite element for the dynamic analysis of sandwich shells of revolution. *Engineering Computations* 1996;13:240–62.