Beam Structural Modelling in Hydroelastic Analysis of Ultra Large Container Ships

1. Introduction

Ultra large container ships are very sensitive to the wave load of quartering seas due to considerably reduced torsional stiffness caused by large deck openings. As a result, their natural frequencies can fall into the range of encounter frequencies in an ordinary sea spectrum. Therefore, the wave induced hydroelastic response of large container ships becomes an important issue in structural design. Mathematical hydroelastic model incorporates structural, hydrostatic and hydrodynamic parts (Senjanović et al. 2007, 2008a, 2009b, 2010b). Beam structural model is preferable in the early design stage and for determining global response, while for more detailed analyses 3D FEM model has to be used. The hydroelastic analysis is performed by the modal superposition method, which requires dry natural vibrations of the structure to be determined. For each mode dynamic coefficients (added mass and damping) and wave load are calculated based on velocity potential. The governing equation of ship motion in rough sea specified for the impulsive (slamming) load as a transient problem is solved in time domain. The motion equation is also given for the case of harmonic wave excitation (springing), which is solved in the frequency domain.

In the chapter, methodology of the ship hydroelastic analysis is described, and position and role of the beam structural model is explained. Beam finite element for coupled horizontal and torsional vibrations, that includes warping of ship cross-section, is constructed. Shear influence on both bending and torsion is taken into account. The strip element method is used for determination of normal and shear stress flows, and stiffness moduli, i.e. shear area, torsional modulus, shear inertia modulus (as a novelty), and warping modulus.

In the modelling of large container ships it is important to appropriately account for the contribution of transverse bulkheads to hull stiffness and the behavior of relatively short engine room structure. In the former case, the equivalent torsional modulus is determined by increasing ordinary (St. Venant) value, depending on the ratio of the strain energy of a bulkhead and corresponding hull portion. Equivalent torsional modulus of the engine room structure is also determined utilizing the energy approach. It is assumed that a short closed structure behaves as an open one with the contribution of decks.

Application of the beam structural model for ship hydroelastic analysis is illustrated in case of a very large container ship. Correlation of dry natural vibrations analysis results for the beam model with those for 3D FEM model shows very good agreement. Hydroelastic analysis emphasizes peak values of transfer functions of displacements and sectional forces in resonances, i.e. in the case when the encounter frequency is equal to one of the natural frequencies.

2. Methodology of ship hydroelastic analysis

A structural model, ship and cargo mass distributions and geometrical model of ship surface have to be defined to perform ship hydroelastic analysis. At the beginning, dry natural vibrations have to be calculated, and after that modal hydrostatic stiffness, modal added mass, damping and modal wave load are determined. Finally, wet natural vibrations as well as the transfer functions (RAO) for determining ship structural response to wave excitation are obtained (Senjanović et al. 2008a, 2009b), Fig. 1.

Figure 1.

Methodology of ship hydroelastic analysis

3. General remarks on structural model

A ship hull, as an elastic non-prismatic thin-walled girder, performs longitudinal, vertical, horizontal and torsional vibrations. Since the cross-sectional centre of gravity and centroid, as well as the shear centre positions are not identical, coupled longitudinal and vertical, and horizontal and torsional vibrations occur, respectively. The shear centre in ships with large hatch openings is located below the keel and therefore the coupling of horizontal and torsional vibrations is extremely high. The above problem is rather complex due to geometrical discontinuity of the hull cross-section, Fig. 2.

The accuracy of the solution depends on the reliability of stiffness parameters determination, i.e. of bending, shear, torsional and warping moduli. The finite element method is a powerful tool to solve the above problem in a successful way. One of the first solutions for coupled horizontal and torsional hull vibrations, dealing with the finite element technique, is given in (Kawai, 1973, Senjanović & Grubišić, 1991). Generalised and improved solutions are presented in (Pedersen, 1985, Wu & Ho, 1987). In all these references, the determination of hull stiffness is based on the classical thin-walled girder theory, which doesn’t give a satisfactory value for the warping modulus of the open cross-section (Haslum & Tonnessen, 1972, Vlasov, 1961). Apart from that, the fixed values of stiffness moduli are determined, so that the application of the beam theory for hull vibration analysis is limited to a few lowest natural modes only. Otherwise, if the mode dependent stiffness parameters are used the application of the beam theory can be extended up to the tenth natural mode (Senjanović & Fan, 1989, 1992, 1997).

.

Figure 2.

Discontinuities of ship hullMethodology of ship hydroelastic analysis

4. Consistent differential equations of beam vibrations

Referring to the flexural beam theory (Timoshenko & Young, 1955, Senjanović, 1990), the total beam deflection, w, consists of the bending deflection, w _b , and the shear deflection, w _s , while the angle of cross-section rotation depends only on the former, Fig. 3

$$w=w_b+w_s,\phi =\frac{\partial w_b}{\partial x}$$ E1

The cross-sectional forces are the bending moment and the shear force

$$M=-E{I}_b\frac{\partial \phi }{\partial x},$$ E2

$$Q=G{A}_s\frac{\partial w_s}{\partial x},$$ E3

where E and G are the Young's and shear modulus, respectively, while I _b and A _s are the moment of inertia of cross-section and shear area, respectively.

The inertia load consists of the distributed transverse load, q _i , and the bending moment, μ _i , and in the case of coupled horizontal and torsional vibration is specified as

$$q_i=-m\left(\frac{{\partial }^{2}w}{\partial t^2}+c\frac{{\partial }^2\psi }{\partial t^2}\right),$$ E4

$${\mu }_i=-J_b\frac{{\partial }^2\phi }{\partial t^2},$$ E5

where m is the distributed mass, J _b is the mass moment of inertia about z-axis, and c is the distance between the centre of gravity and the shear centre, $$c=z_G-z_S$$ , Fig. 4.

Figure 3.

Beam bending and torsion

Figure 4.

Cross-section of a thin-walled girder

In a similar way the total twist angle, ψ, consists of the pure twist angle, ψ _t , and the shear contribution, ψ _s , while the second torsional displacement, which causes warping of cross-section, is variation of the pure twist angle, i.e. Fig. 3 (Pavazza, 2005)

$$\psi ={\psi }_t+{\psi }_s,\vartheta =\frac{\partial {\psi }_t}{\partial x}.$$ E6

The cross-sectional forces include the pure torsional torque, T _t , warping bimoment, B _w , and additional torque due to restrained warping, T _w

$$T_t=G{I}_t\vartheta ,$$ E7

$$B_w=-E{I}_w\frac{\partial \vartheta }{\partial x},$$ E8

$$T_w=G{I}_s\frac{\partial {\psi }_s}{\partial x},$$ E9

where I _t , I _w and I _s are the torsional modulus, warping modulus and shear inertia modulus, respectively.

The inertia load consists of the distributed torque, μ _ti , and the bimoment, b _i , presented in the following form:

$${\mu }_{ti}=-J_t\frac{{\partial }^2\psi }{\partial t^2}-mc\frac{{\partial }^{2}w}{\partial t^2},$$ E10

$$b_i=-J_w\frac{{\partial }^2\vartheta }{\partial t^2},$$ E11

where J _t is the mass polar moment of inertia about the shear centre, and J _w is the mass bimoment of inertia with respect to the warping centre, Fig. 4.

Considering the equilibrium of a beam differential element, one can write for flexural vibrations

$$\frac{\partial M}{\partial x}=Q+{\mu }_i,$$ E12

$$\frac{\partial Q}{\partial x}=-q_i-q,$$ E13

and for torsional vibrations (Pavazza, 1991)

$$\frac{\partial B_w}{\partial x}=T_w+b_i,$$ E14

$$\frac{\partial T_t}{\partial x}+\frac{\partial T_w}{\partial x}=-{\mu }_{ti}-\mu .$$ E15

The above equations can be reduced to two coupled partial differential equations as follows. Substituting Eqs. (2) and (3) into (12) yields

$$\frac{\partial w_s}{\partial x}=-\frac{E{I}_b}{G{A}_s}\frac{{\partial }^2\phi }{\partial x^2}+\frac{J_b}{G{A}_s}\frac{{\partial }^2\phi }{\partial t^2}.$$ E16

By inserting Eqs. (3) and (4) into (13) leads to

$$E{I}_b\frac{{\partial }^4\phi }{\partial x^4}+m\frac{{\partial }^2\phi }{\partial t^2}-\left(J_b+m\frac{E{I}_b}{G{A}_s}\right)\frac{{\partial }^4\phi }{\partial x^2\partial t^2}+\frac{m{J}_b}{G{A}_s}\frac{{\partial }^4\phi }{\partial t^4}+mc\frac{{\partial }^3\psi }{\partial x\partial t^2}=\frac{\partial q}{\partial x}.$$ E17

In a similar way, substituting Eqs. (8) and (9) into (14) yields

$$\frac{\partial {\psi }_s}{\partial x}=-\frac{E{I}_w}{G{I}_s}\frac{{\partial }^2\vartheta }{\partial x^2}+\frac{J_w}{G{I}_s}\frac{{\partial }^2\vartheta }{\partial t^2}.$$ E18

By inserting Eqs. (7), (9) and (10) into (15) one finds

$$E{I}_w\frac{{\partial }^4\vartheta }{\partial x^4}-G{I}_t\frac{{\partial }^2\vartheta }{\partial x^2}+J_t\frac{{\partial }^2\vartheta }{\partial t^2}-\left(J_w+J_t\frac{E{I}_w}{G{I}_s}\right)\frac{{\partial }^4\vartheta }{\partial x^2\partial t^2}+\frac{J_w}{G{I}_s}\frac{{\partial }^4\vartheta }{\partial t^4}+mc\frac{{\partial }^{3}w}{\partial x\partial t^2}=\frac{\partial \mu }{\partial x}.$$ E19

Furthermore, ψ in (17) can be split into $${\psi }_t+{\psi }_s$$ and the later term can be expressed with (18). Similar substitution can be done for $$w=w_b+w_s$$ in (19), where w _s is given with (16). Thus, taking into account that $$\phi =\partial w_b/\partial x$$ and $$\vartheta =\partial {\psi }_t/\partial x$$ , Eqs. (17) and (19) after integration per x read

$$\begin{array}{l}E{I}_b\frac{{\partial }^4{w}_b}{\partial x^4}+m\frac{{\partial }^2{w}_b}{\partial t^2}-\left(J_b+m\frac{E{I}_b}{G{A}_s}\right)\frac{{\partial }^4{w}_b}{\partial x^2\partial t^2}+\frac{m{J}_b}{G{A}_s}\frac{{\partial }^4{w}_b}{\partial t^4}\\ +mc\left(\frac{{\partial }^2{\psi }_t}{\partial t^2}-\frac{E{I}_w}{G{I}_s}\frac{{\partial }^4{\psi }_t}{\partial x^2\partial t^2}+\frac{J_w}{G{I}_s}\frac{{\partial }^4{\psi }_t}{\partial t^4}\right)=q\end{array}$$ E20

$$\begin{array}{l}E{I}_w\frac{{\partial }^4{\psi }_t}{\partial x^4}-G{I}_t\frac{{\partial }^2{\psi }_t}{\partial x^2}+J_t\frac{{\partial }^2{\psi }_t}{\partial t^2}-\left(J_w+J_t\frac{E{I}_w}{G{I}_s}\right)\frac{{\partial }^4{\psi }_t}{\partial x^2\partial t^2}\\ +\frac{J_w}{G{I}_s}\frac{{\partial }^4{\psi }_t}{\partial t^4}+mc\left(\frac{{\partial }^2{w}_b}{\partial t^2}-\frac{E{I}_b}{G{A}_s}\frac{{\partial }^4{w}_b}{\partial x^2\partial t^2}+\frac{J_b}{G{A}_s}\frac{{\partial }^4{w}_b}{\partial t^4}\right)=\mu .\end{array}$$ E21

After solving Eqs. (20) and (21) the total deflection and twist angle are obtained by employing (16) and (18), i.e.

$$w=w_b+w_s=w_b-\frac{E{I}_b}{G{A}_s}\frac{{\partial }^2{w}_b}{\partial x^2}+\frac{J_b}{G{A}_s}\frac{{\partial }^2{w}_b}{\partial t^2}+f\left(t\right)$$ E22

$$\psi ={\psi }_t+{\psi }_s={\psi }_t-\frac{E{I}_w}{G{I}_s}\frac{{\partial }^2{\psi }_t}{\partial x^2}+\frac{J_w}{G{I}_s}\frac{{\partial }^2{\psi }_t}{\partial t^2}+g\left(t\right)$$ E23

where f(t) and g(t) are integration functions, which depend on initial conditions.

The main purpose of developing differential equations of vibrations (20) and (21) is to get insight into their constitution, position and role of the stiffness and mass parameters, and coupling, which is realized through the inertia terms. If the pure torque T _t is excluded from the above theoretical consideration, it is obvious that the complete analogy between bending and torsion exists, (Pavazza, 1991).

Application of Eqs. (20) and (21) is limited to prismatic girders. For more complex problems, like ship hull, the finite element method is on disposal.

The shape functions of beam finite element for vibration analysis have to satisfy the following consistency relations for harmonic vibrations obtained from Eqs. (22) and(23), (Senjanović, 1990)

$$w=w_b+w_s=\left(1-{\omega }^2\frac{J_b}{G{A}_s}\right)w_b-\frac{E{I}_b}{G{A}_s}\frac{{\text{d}}^2{w}_b}{\text{d}{x}^2}$$ E24

$$\psi ={\psi }_t+{\psi }_s=\left(1-{\omega }^2\frac{J_w}{G{I}_s}\right){\psi }_t-\frac{E{I}_w}{G{I}_s}\frac{{\text{d}}^2{\psi }_t}{\text{d}{x}^2}.$$ E25

5. Beam finite element

The properties of a finite element for the coupled horizontal and torsional vibration analysis can be derived from the total element energy. It consists of the strain energy, the kinetic energy, the work of the distributed external lateral load, q, and the torque, μ, and the work of the boundary forces. Thus, according to (Senjanović, 1990, Senjanović & Grubišić, 1991),

$$\begin{array}{c}{E}_{tot}=\frac{1}{2}{\displaystyle \underset{0}{\overset{l}{\int }}\left[E{I}_b{\left(\frac{{\partial }^2{w}_b}{\partial x^2}\right)}^2+G{A}_s{\left(\frac{\partial w_s}{\partial x}\right)}^2+E{I}_w{\left(\frac{{\partial }^2{\psi }_t}{\partial x^2}\right)}^2+G{I}_s{\left(\frac{\partial {\psi }_s}{\partial x}\right)}^2+G{I}_t{\left(\frac{\partial {\psi }_t}{\partial x}\right)}^2\right]\text{d}x}\\ +\frac{1}{2}{\displaystyle \underset{0}{\overset{l}{\int }}\left[m{\left(\frac{\partial w}{\partial t}\right)}^2+J_b{\left(\frac{{\partial }^2{w}_b}{\partial x\partial t}\right)}^2+2mc\frac{\partial w}{\partial t}\frac{\partial \psi }{\partial t}+J_w{\left(\frac{{\partial }^2{\psi }_t}{\partial x\partial t}\right)}^2+J_t{\left(\frac{\partial \psi }{\partial t}\right)}^2\right]\text{d}x}\\ -{\displaystyle \underset{0}{\overset{l}{\int }}(qw+\mu \psi )}\text{d}x+{(}_Q^w,\end{array}$$ E26

where l is the element length. Since the beam has four displacements, $$w,\phi ,\psi ,\vartheta$$ , a two-node finite element has eight degrees of freedom, i.e. four nodal shear-bending and torsion-warping displacements respectively, Fig. 5,

$$\left\{U\right\}=\left\{\begin{array}{c}w(0)\\ \phi (0)\\ w(l)\\ \phi (l)\end{array}\right\},\left\{V\right\}=\left\{\begin{array}{c}\psi (0)\\ \vartheta (0)\\ \psi (l)\\ \vartheta (l)\end{array}\right\}$$ E27

Therefore, the basic beam displacements, w _b and ψ _t , can be presented as the third-order polynomials

$$\begin{array}{l}{w}_b=⟨a_k⟩\left\{{\xi }^k\right\},{\psi }_t=⟨d_k⟩\left\{{\xi }^k\right\},k=\mathrm{0,}\mathrm{1,}\mathrm{2,}3,\\ \xi =\frac{x}{l},⟨...⟩={\left\{...\right\}}^T.\end{array}$$ E28

Furthermore, satisfying alternately the unit value for one of the nodal displacement {U} and zero values for the remaining displacements, and doing the same for {V}, it follows that:

$$\begin{array}{l}{w}_b=⟨w_{bi}⟩\left\{U\right\},w_s=⟨w_{si}⟩\left\{U\right\},w=⟨w_i⟩\left\{U\right\},\\ {\psi }_t=⟨{\psi }_{ti}⟩\left\{V\right\},{\psi }_s=⟨{\psi }_{si}⟩\left\{V\right\},\psi =⟨{\psi }_i⟩\left\{V\right\},i=\mathrm{1,}\mathrm{2,}\mathrm{3,}\mathrm{4,}\end{array}$$ E29

where w _bi , w _si , w _i and ψ _ti , ψ _si , ψ _i are the shape functions specified below by employing relations (24) and (25)

$$\begin{array}{l}{w}_{bi}=⟨a_{ik}⟩\left\{{\xi }^k\right\},w_{si}=⟨b_{ik}⟩\left\{{\xi }^k\right\},w_i=⟨c_{ik}⟩\left\{{\xi }^k\right\}\\ {\psi }_{ti}=⟨d_{ik}⟩\left\{{\xi }^k\right\},{\psi }_{si}=⟨e_{ik}⟩\left\{{\xi }^k\right\},{\psi }_i=⟨f_{ik}⟩\left\{{\xi }^k\right\}\end{array}$$ E30

$$\left[a_{ik}\right]=\frac{1}{\alpha \left(\alpha +12\beta \right)}\left[\begin{array}{cccc}\alpha +6\beta & 0& -3\alpha & 2\alpha \\ -4\beta \left(\alpha +3\beta \right)l& \alpha \left(\alpha +12\beta \right)l& -2\alpha \left(\alpha +3\beta \right)l& {\alpha }^{2}l\\ 6\beta & 0& 3\alpha & -2\alpha \\ -2\beta \left(\alpha -6\beta \right)l& 0& -\alpha \left(\alpha -6\beta \right)l& {\alpha }^{2}l\end{array}\right]$$ E31

Figure 5.

Beam finite elemet

$$\begin{array}{l}{b}_{i0}=-\left(1-\alpha \right)a_{i0}-2\beta a_{i2}\\ b_{i1}=-\left(1-\alpha \right)a_{i1}-6\beta a_{i3}\\ b_{i2}=-\left(1-\alpha \right)a_{i2}\\ b_{i3}=-\left(1-\alpha \right)a_{i3},\end{array}$$ E32

$$\left[c_{ik}\right]=\left[a_{ik}\right]+\left[b_{ik}\right],i=\mathrm{1,2,3,4,}k=\mathrm{0,1,2,3}$$ E33

$$\alpha =1-{\omega }^2\frac{J_b}{G{A}_s},\beta =\frac{E{I}_b}{G{A}_s{l}^2}$$ E34

Constitution of torsional matrices $$\left[d_{ik}\right]$$ , $$\left[e_{ik}\right]$$ and $$\left[f_{ik}\right]$$ is the same as $$\left[a_{ik}\right]$$ , $$\left[b_{ik}\right]$$ and $$\left[c_{ik}\right]$$ , but parameters α and β have to be exchanged with

$$\eta =1-{\omega }^2\frac{J_w}{G{I}_s},\gamma =\frac{E{I}_w}{G{I}_s{l}^2}$$ E35

according to (25). By substituting Eqs. (29) (26) one obtains

$$E_{tot}=\frac{1}{2}{\left\{\begin{array}{c}U\\ V\end{array}\right\}}^T\left[\begin{array}{cc}{k}_{bs}& 0\\ 0& k_{ws}+k_t\end{array}\right]\left\{\begin{array}{c}U\\ V\end{array}\right\}+\frac{1}{2}{\left\{\begin{array}{c}\dot{U}\\ \dot{V}\end{array}\right\}}^T\left[\begin{array}{cc}{m}_{sb}& m_{st}\\ m_{ts}& m_{tw}\end{array}\right]\left\{\begin{array}{c}\dot{U}\\ \dot{V}\end{array}\right\}-{\left\{\begin{array}{c}q\\ \mu \end{array}\right\}}^T\left\{\begin{array}{c}U\\ V\end{array}\right\}-{\left\{\begin{array}{c}P\\ R\end{array}\right\}}^T\left\{\begin{array}{c}U\\ V\end{array}\right\},$$ E36

where, assuming constant values of the element properties, $${\left[k\right]}_{bs}=\left[E{I}_b{\displaystyle \underset{0}{\overset{l}{\int }}\frac{{\text{d}}^2{w}_{bi}}{\text{d}{x}^2}\frac{{\text{d}}^2{w}_{bj}}{\text{d}{x}^2}\text{d}x}+G{A}_s{\displaystyle \underset{0}{\overset{l}{\int }}\frac{\text{d}{w}_{si}}{\text{d}x}\frac{\text{d}{w}_{sj}}{\text{d}x}\text{d}x}\right]$$ – bending-shear stiffness matrix,

$${\left[k\right]}_{ws}=\left[E{I}_w{\displaystyle \underset{0}{\overset{l}{\int }}\frac{{\text{d}}^2{\psi }_{ti}}{\text{d}{x}^2}\frac{{\text{d}}^2{\psi }_{tj}}{\text{d}{x}^2}\text{d}x}+G{I}_s{\displaystyle \underset{0}{\overset{l}{\int }}\frac{\text{d}{\psi }_{si}}{\text{d}x}\frac{\text{d}{\psi }_{sj}}{\text{d}x}\text{d}x}\right]$$ – warping-shear stiffness matrix,

$${\left[k\right]}_t=\left[G{I}_t{\displaystyle \underset{0}{\overset{l}{\int }}\frac{\text{d}{\psi }_{ti}}{\text{d}x}\frac{\text{d}{\psi }_{tj}}{\text{d}x}\text{d}x}\right]$$ – torsion stiffness matrix,

$${\left[m\right]}_{sb}=\left[m{\displaystyle \underset{0}{\overset{l}{\int }}{w}_i{w}_j\text{d}x}+J_b{\displaystyle \underset{0}{\overset{l}{\int }}\frac{\text{d}{w}_{bi}}{\text{d}x}\frac{\text{d}{w}_{bj}}{\text{d}x}\text{d}x}\right]$$ – shear-bending mass matrix,

$${\left[m\right]}_{tw}=\left[J_t{\displaystyle \underset{0}{\overset{l}{\int }}{\psi }_i{\psi }_j\text{d}x}+J_w{\displaystyle \underset{0}{\overset{l}{\int }}\frac{\text{d}{\psi }_{ti}}{\text{d}x}\frac{\text{d}{\psi }_{tj}}{\text{d}x}\text{d}x}\right]$$ – torsion-warping mass matrix,

$${\left[m\right]}_{st}=\left[mc{\displaystyle \underset{0}{\overset{l}{\int }}{w}_i{\psi }_j\text{d}x}\right],{\left[m\right]}_{ts}={\left[m\right]}_{st}^T$$ – shear-torsion mass matrix,

$$\begin{array}{l}\left\{q\right\}=\left\{{\displaystyle \underset{0}{\overset{l}{\int }}q{w}_j\text{d}x}\right\}–\text{shearloadvector}\text{,}\\ \left\{\mu \right\}=\left\{{\displaystyle \underset{0}{\overset{l}{\int }}\mu {\psi }_j\text{d}x}\right\}–\text{torsionloadvector}\text{,}\\ i,j=\mathrm{1,}\mathrm{2,}\mathrm{3,}4.\end{array}$$ E37

The vectors {P} and {R} in Eq. (36) represent the shear-bending and torsion-warping nodal forces, respectively,

$$\left\{P\right\}=\left\{\begin{array}{c}-Q(0)\\ M(0)\\ Q(l)\\ -M(l)\end{array}\right\},\left\{R\right\}=\left\{\begin{array}{c}-T(0)\\ B_w(0)\\ T(l)\\ -B_w(l)\end{array}\right\}$$ E38

The above matrices are specified in Appendix A, as well as the load vectors for linearly distributed loads along the finite element, i.e.

$$q=q_0+q_1\xi ,\mu ={\mu }_0+{\mu }_1\xi .$$ E39

The total element energy has to be at its minimum. Satisfying the relevant conditions

$$\frac{\partial E_{tot}}{\partial \left\{U\right\}}=\left\{0\right\},\frac{\partial E_{tot}}{\partial \left\{V\right\}}=\left\{0\right\}$$ E40

and employing Lagrange equations of motion, the finite element equation yields

$$\left\{f\right\}=\left[k\right]\left\{\delta \right\}+\left[m\right]\left\{\ddot{\delta }\right\}-{\left\{f\right\}}_{q\mu },$$ E41

where

$$\begin{array}{l}\left\{f\right\}=\left\{\begin{array}{c}P\\ R\end{array}\right\},{\left\{f\right\}}_{q\mu }=\left\{\begin{array}{c}q\\ \mu \end{array}\right\},\left\{\delta \right\}=\left\{\begin{array}{c}U\\ V\end{array}\right\}\\ \left[k\right]=\left[\begin{array}{cc}{k}_{bs}& 0\\ 0& k_{ws}+k_t\end{array}\right],\left[m\right]=\left[\begin{array}{cc}{m}_{sb}& m_{st}\\ m_{ts}& m_{tw}\end{array}\right].\end{array}$$ E42

It is obvious that coupling between the bending and torsion occurs through the mass matrix only, i.e. by the coupling matrices [m] _st and [m] _ts .

6. Contribution of transverse bulkheads to the hull stiffness

This problem for container ships is extensively analyzed in (Senjanović et al., 2008b), where torsional modulus of ship cross-section is increased proportionally to the ratio of bulkhead strain energy and strain energy of corresponding hull portion. The bulkhead is considered as an orthotropic plate with very strong stool (Szilard, 2004). Bulkhead strain energy is determined for the given warping of cross-section as a boundary condition. The warping causes bulkhead screwing and bending. Here, only the review of the final results is presented. Bulkhead deflection (axial displacement) is given by the following formula, Fig. 6:

$$u\left(y,z\right)=-y\left\{\left(z-d\right)+\left[1-{\left(\frac{y}{b}\right)}^2\right]\frac{z^2}{H}\left(2-\frac{z}{H}\right)\right\}{\psi }^{\prime },$$ E43

where H is the ship height, b is one half of bulkhead breadth, d is the distance of warping centre from double bottom neutral line, y and z are transverse and vertical coordinates, respectively, and $${\psi }^{\prime }$$ is the variation of twist angle.

Figure 6.

Shape of bulkhead deformation

The bulkhead grillage strain energy includes vertical and horizontal bending with contraction and torsion (Senjanović et al., 2008b).

$$U_g=\frac{1}{1-{\nu }^2}[\frac{116H^3}{35b}{i}_y+\frac{32b^3}{105H}{i}_z+\frac{8Hb}{75}\nu \left(i_y+i_z\right)+\frac{143Hb}{75}\left(1-\nu \right)i_t]E{{\psi }^{\prime }}^2$$ E44

where i _y , i _z and i _t are the average moments of inertia of cross-section and torsional modulus per unit breadth, respectively. The stool strain energy is comprised of the bending, shear and torsional contributions

$$U_s=\left[\frac{12h^2{I}_{sb}}{b}+72\left(1+\nu \right)\frac{h^2}{b^3}\frac{I_{sb}^2}{A_s}+\frac{9b{I}_{st}}{10\left(1+\nu \right)}\right]E{{\psi }^{\prime }}^2$$ E45

where I _sb , A _s and I _st are the moment of inertia of cross-section, shear area and torsional modulus, respectively. Quantity h is the stool distance from the inner bottom, Fig. 7.

Figure 7.

Longitudinal section of container ship hold

The equivalent torsional modulus yields, Fig. 7

$$I_t^{\ast }=\left[1+\frac{a}{l_1}+\frac{4\left(1+\nu \right)C}{I_t{l}_0}\right]I_t,C=\frac{U_g+U_s}{E{{\psi }^{\prime }}^2},$$ E46

where a is the web height of bulkhead girders (frame spacing), l ₀ is the bulkhead spacing, $$l_1=l_0-a$$ is the net length, and C is the energy coefficient.

7. Contribution of engine room structure to the hull stiffness

Ultra Large Container Ships are characterized by relatively short engine room structure with length of about a half of ship breadth. Its complex deformation is illustrated in a case of a 7800 TEU container ship, Fig. 8. The deck shear deformation is predominant, while hold transverse bulkhead stool is exposed to bending. Due to shortness of the engine room, its transverse bulkheads are skewed but somewhat less pronounced than warping of the hold bulkheads. Warping of the transom is negligible, and that is an important fact when specifying boundary conditions in vibration analysis.

Figure 8.

Deformation of 7800 TEU container ship aft structure

7.1. Stiffness of engine room structure

A short engine room structure can be considered either as a closed segment with relevant stiffness or as an open segment with increased stiffness due to deck contribution (Pedersen, 1985). The latter simulation in fact gives results which agree better with 3D FEM results, than the former one (Pedersen, 1983). Deck contribution to hull stiffness can be determined by energy approach, as it is done in the case of transverse bulkheads (Senjanović et al., 2008b). Such a beam model is consistent at global level of energy balance, and that is sufficient for application in ship hydroelastic analysis, where proper natural frequencies and mode shapes of dry hull are required.

In the case of short engine room, torsion induces distortion of cross-section while hull bending is negligible. Solution of that complex problem is described here by employing the energy balance approach and concept of the effective stiffness due to reason of simplicity. A closed hull segment is considered as open one with deck influence. For that purpose let us determine deck strain energy. All quantities related to closed and open cross-section are designated by $${(.)}^{\ast }$$ and $${(.)}^{\circ }$$ , respectively

As it can be seen in Fig. 8, the upper deck is exposed to large deformation, while the double bottom in-plane deformation is quite small. The relative axial displacement of the internal upper deck boundaries, with respect to double bottom, is result of their warping

$$U=\left|U_D\right|+\left|U_B\right|=\left(\left|w_D\right|+\left|w_B\right|\right){{\psi }^{\prime }}_t$$ E47

It causes deck in-plane (membrane) deformation. The problem can be solved in an approximate analytical way by considering deck as a beam. Its horizontal anti-symmetric deflection consists of pure bending and shear contribution, Fig. 9. The former is assumed in the form

$$u_b=\frac{y}{2b}\left[3-{\left(\frac{y}{b}\right)}^2\right]U_b,$$ E48

which satisfies relevant boundary conditions: $$u_b\left(0\right)=0$$ and $${u^{″}}_b\left(0\right)=0$$ , where $$U_b$$ is the boundary bending deflection. Shear deflection depends on bending deflection

$$u_s=-\frac{EI}{GA}\frac{d^2{u}_b}{d{y}^2}=2\left(1+\nu \right){\left(\frac{a}{b}\right)}^2\frac{y}{b}{U}_b,$$ E49

where the internal deck cross-section area, $$A=2at$$ , its moment of inertia, $$I=\frac{2}{3}{a}^{3}t$$ , and the relation $$E=2\left(1+\nu \right)G$$ , are taken into account, Fig. 9. Total deflection is obtained by summing up its constitutive parts, Eqs. (48) and (49). Relation between total boundary deflection and the bending boundary deflection reads

$$U=\left[1+2\left(1+\nu \right){\left(\frac{a}{b}\right)}^2\right]U_b$$ E50

Figure 9.

Deck deformation and double bottom rotation, a)-bird view, b)-lateral view

The total internal deck strain energy consists of the bending and shear contributions

$$E_1=\frac{1}{2}EI{{\displaystyle \underset{-b}{\overset{b}{\int }}\left(\frac{d^2{u}_b}{d{y}^2}\right)}}^{2}dy+\frac{1}{2}GA{{\displaystyle \underset{-b}{\overset{b}{\int }}\left(\frac{d{u}_s}{dy}\right)}}^{2}dy$$ E51

By substituting Eqs. (48) and (49) into (51), one finds

$$E_1=4\left(1+\nu \right)Gt{\left(\frac{a}{b}\right)}^3\left[1+2\left(1+\nu \right){\left(\frac{a}{b}\right)}^2\right]U_b^2$$ E52

Finally, by taking into account Eqs. (47) and (50), yields

$$E_1=\frac{4\left(1+\nu \right)Gt{\left(\frac{a}{b}\right)}^3}{1+2\left(1+\nu \right){\left(\frac{a}{b}\right)}^2}{\left(\left|w_D\right|+\left|w_B\right|\right)}^2{{\psi }^{\prime }}_t^2$$ E53

On the other hand, total energy of the closed hull segment can be obtained by summing up energy of open segment and the deck strain energy, i.e.

$$E_{tot}^{\ast }=E_w^{\circ }+E_t^{\circ }+E_1^{}-E_{\mu }^{}$$ E54

where

$$E_w^{\circ }=\frac{1}{2}{\displaystyle \underset{-a}{\overset{a}{\int }}{B}_w^{\circ }{{\psi }^{″}}_t}dx,E_t^{\circ }=\frac{1}{2}{\displaystyle \underset{-a}{\overset{a}{\int }}{T}_t^{\circ }{{\psi }^{\prime }}_t}dx,E_{\mu }^{}={\displaystyle \underset{-a}{\overset{a}{\int }}{\mu }_x^{}\psi }dx.$$ E55

Within a short span 2a, constant value of $${{\psi }^{\prime }}_t$$ (as for deck) can be assumed, so that second term in Eq. (26) by inserting $$T_t^{\circ }$$ from Eqs. (7), leads to

$$E_t^{\circ }=G{I}_t^{\circ }a{{\psi }^{\prime }}_t^2.$$ E56

$$E_t^{\circ }$$ and $$E_1$$ in (54) can be unified into one term since both depend on

$${{\psi }^{\prime }}_t^2$$

$$E_t^{\circ }+E_1=Ga{\tilde{I}}_t{{\psi }^{\prime }}_t^2$$ E57

where

$${\tilde{I}}_t=\left(1+C\right)I_t^{\circ },C=\frac{E_1}{E_t^{\circ }}$$ E58

$${\tilde{I}}_t$$ is the effective torsional modulus which includes both open cross-section and deck effects.

Engine room structure is designed in such a way that the hold double skin continuity is ensured and necessary decks are inserted between the double skins. Strain energy is derived for the first (main) deck and for the others it can be assumed that their strain energy is proportional to the deck plating volume, V, and linearly increasing deformation with the deck distance from inner bottom, h, Fig. 9, since the double bottom is much stiffer than decks. In that way the coefficient C, Eq. (58b), by employing (53) and (56), reads

$$C=\frac{{\displaystyle \sum E_i}}{E_t}=\frac{4\left(1+\nu \right)t_1{\left(\frac{a}{b}\right)}^3{\left(\left|w_D\right|+\left|w_B\right|\right)}^{2}k}{\left[1+2\left(1+\nu \right){\left(\frac{a}{b}\right)}^2\right]I_t^{\circ }a}$$ E59

where

$$k={\displaystyle \sum \frac{V_i}{V_1}}{\left(\frac{h_i}{h_1}\right)}^2.$$ E60

In the above consideration distortion of cross-sections is not included and that is subject of further investigation.

7.2. Torsion of segmented girder

Let us consider a girder consisted of three segments, Fig. 10. The end segments are open and the middle one is closed, so that the girder is symmetric with respect to the z axis. Each segment is specified in its local coordinate system. The properties of the middle and end segments are designated by $${(.)}^{\ast }$$ and $${(.)}^{\circ }$$ , respectively. The relevant expressions for displacements and sectional forces are listed below (Senjanović et al., 2009a, 2010a):

$$\begin{array}{l}u=w\frac{\text{d}{\psi }_t}{\text{d}x}=w\left(\frac{A_1}{l}+A_2\alpha \text{sh}\alpha x+A_3\alpha \text{ch}\alpha x+{{\psi }^{\prime }}_p\right),\\ T_t=G{I}_t\left(\frac{A_1}{l}+A_2\alpha \text{sh}\alpha x+A_3\alpha \text{ch}\alpha x+{{\psi }^{\prime }}_p\right),\\ T_w=-G{I}_t\left(A_2\alpha \text{sh}\alpha x+A_3\alpha \text{ch}\alpha x\right)-E{I}_w{{\psi }^{‴}}_p,\\ T=G{I}_t\left(\frac{A_1}{l}+{{\psi }^{\prime }}_p\right)-E{I}_w{{\psi }^{‴}}_p,\end{array}$$ E61

$$B_w=G{I}_t\left(A_2\text{ch}\alpha x+A_3\text{sh}\alpha x\right)+E{I}_w{{\psi }^{″}}_p,$$ E62

where $${\psi }_p$$ represents particular solution of differential equation and coefficient $$\alpha$$ yields

$$\alpha =\sqrt{\frac{G{I}_t}{E{I}_w}}.$$ E63

The symbols $$A_i$$ and $$B_i$$ are used for the integration constants of the closed and open segments. The girder is loaded with torque M _t at the ends, while $${\mu }_x=0$$ . The ends are fixed against warping.

The boundary and compatibility conditions in the considered case, yield

$$\begin{array}{l}{\psi }^{\ast }\left(a\right)={\psi }^{\circ }\left(0\right),{\psi }_t^{\ast }{}^{\prime }\left(a\right)={\psi }_t^{\circ }{}^{\prime }\left(0\right),\\ T^{\ast }\left(a\right)=T^{\circ }\left(0\right),B_w^{\ast }\left(a\right)=B_w^{\circ }\left(0\right),\\ u^{\circ }\left(l^{\circ }\right)=\mathrm{0,}{T}^{\circ }\left(l^{\circ }\right)=M_t.\end{array}$$ E64

Figure 10.

Torsion of segmented girder

From the third and last conditions (63) one finds

$$A_1=\frac{M_{t}a}{G{\tilde{I}}_t^{}},B_1=\frac{M_t{l}^{\circ }}{G{I}_t^{\circ }}.$$ E65

The remaining four conditions (63) lead to the system of algebraic equations (Senjanović et al., 2010a) and its analytical solution reads:

$$A_3=\frac{D_{A3}}{D},B_2=\frac{D_{B2}}{D},B_3=\frac{D_{B3}}{D},$$ E66

where

$$\begin{array}{l}{D}_{A3}=\frac{M_t}{G}\left[\left(1-\frac{I_t^{\circ }}{{\tilde{I}}_t}\right)\text{ch}\beta l^{\circ }-1\right],\\ D_{B2}=-\frac{M_t}{G}\text{sh}\alpha a\left[\left(1-\frac{{\tilde{I}}_t}{I_t^{\circ }}\right)\text{ch}\beta l^{\circ }+\frac{{\tilde{I}}_t}{I_t^{\circ }}\right],\\ D_{B3}=\frac{M_t}{G}\left[\left(1-\frac{{\tilde{I}}_t}{I_t^{\circ }}\right)\text{sh}\alpha a\text{sh}\beta l^{\circ }-\frac{\alpha }{\beta }\text{ch}\alpha a\right],\\ D=I_t^{\circ }\alpha \text{ch}\alpha a\text{ch}\beta l^{\circ }+{\tilde{I}}_t\beta \text{sh}\alpha a\text{sh}\beta l^{\circ }.\end{array}$$ E67

8. Numerical procedure for vibration analysis

A thin-walled girder is modelled with a set of beam finite elements. Their assemblage in the global coordinate system, performed in the standard way, results in the matrix equation of motion, which may be extended by the damping forces

$$\left[K\right]\left\{\Delta \right\}+\left[C\right]\left\{\dot{\Delta }\right\}+\left[M\right]\left\{\ddot{\Delta }\right\}=\left\{F(t)\right\},$$ E68

where [K], [C] and [M] are the stiffness, damping and mass matrices, respectively; $$\left\{\Delta \right\},\left\{\dot{\Delta }\right\}\text{and}\left\{\ddot{\Delta }\right\}$$ are the displacement, velocity and acceleration vectors, respectively; and {F(t)} is the load vector.

In case of natural vibration {F(t)} = {0} and the influence of damping is rather low for the most of the structures, so that the damping forces may be ignored. Assuming

$$\left\{\Delta \right\}=\left\{\phi \right\}{\text{e}}^{i\omega t},$$ E69

where $$\left\{\phi \right\}$$ and ω are the mode vector and natural frequency respectively, Eq. (67) leads to the eigenvalue problem

$$\left(\left[K\right]-{\omega }^2\left[M\right]\right)\left\{\phi \right\}=\left\{0\right\}$$ E70

,

which may be solved by employing different numerical methods (Bathe, 1996) The basic one is the determinant search method in which ω is found from the condition

$$\left|\left[K\right]-{\omega }^2\left[M\right]\right|=0$$ E71

by an iteration procedure. Afterwards, $$\left\{\phi \right\}$$ follows from (69) assuming unit value for one element in $$\left\{\phi \right\}$$ .

The forced vibration analysis may be performed by direct integration of Eq. (67), as well as by the modal superposition method. In the latter case the displacement vector is presented in the form

$$\left\{\Delta \right\}=\left[\phi \right]\left\{X\right\}$$ E72

,

where $$\left[\phi \right]=\left[\left\{\phi \right\}\right]$$ is the undamped mode matrix and {X} is the generalised displacement vector. Substituting (71) into (67), the modal equation yields

$$\left[k\right]\left\{X\right\}+\left[c\right]\left\{\dot{X}\right\}+\left[m\right]\left\{\ddot{X}\right\}=\left\{f(t)\right\}$$ E73

,

where

$$\begin{array}{l}\left[k\right]={\left[\phi \right]}^T\left[K\right]\left[\phi \right]\text{–modalstiffnessmatrix}\\ \left[c\right]={\left[\phi \right]}^T\left[C\right]\left[\phi \right]\text{–modaldampingmatrix}\\ \left[m\right]={\left[\phi \right]}^T\left[M\right]\left[\phi \right]\text{–modalmassmatrix}\\ \left\{f(t)\right\}={\left[\phi \right]}^T\left\{F(t)\right\}\text{–modalloadvector}\text{.}\end{array}$$ E74

The matrices [k] and [m] are diagonal, while [c] becomes diagonal only in a special case, for instance if [C] = α ₀ [M] + β ₀ [K], where α ₀ and β ₀ are coefficients (Senjanović, 1990).

Solving (72) for undamped natural vibration, [k] = [ω ² m] is obtained, and by its backward substitution into (72) the final form of the modal equation yields

$$\left[{\omega }^2\right]\left\{X\right\}+2\left[\omega \right]\left[\zeta \right]\left\{\dot{X}\right\}+\left\{\ddot{X}\right\}=\left\{\phi (t)\right\}$$ E75

,

where

$$\begin{array}{l}\left[\omega \right]=\left[\sqrt{\frac{k_{ii}}{m_{ii}}}\right]\text{–naturalfrequencymatrix}\\ \left[\zeta \right]=\left[\frac{c_{ij}}{2\sqrt{(k_{ii}{m}_{ii})}}\right]\text{–relativedampingmatrix}\\ \left\{\phi (t)\right\}=\left\{\frac{f_i(t)}{m_{ii}}\right\}\text{–relativeloadvector}\text{.}\end{array}$$ E76

If [ζ] is diagonal, the matrix Eq. (74) is split into a set of uncoupled modal equations.

If vibration excitation is of periodical nature it can be split into harmonics, and the structure response for each of them is determined in the frequency domain. In a case of general or impulsive excitation the vibration problem has to be solved in the time domain.

Several numerical methods are available for this purpose, as for instance the Houbolt, the Newmark and the Wilson θ method (Bathe, 1996), as well as the harmonic acceleration method (Lozina, 1988, Senjanović, 1984).

It is important to point out that all stiffness and mass matrices of the beam finite element (and consequently those of the assembly) are frequency dependent quantities, due to coefficients α and η in the formulation of the shape functions, Eqs. (34) and (35). Therefore, for solving the eigenvalue problem (69) an iteration procedure has to be applied. As a result of frequency dependent matrices, the eigenvectors are not orthogonal. If they are used in the modal superposition method for determining forced response, full modal stiffness and mass matrices are generated. Since the inertia terms are much smaller than the deformation ones in Eqs. (24) and (25), the off-diagonal elements in modal stiffness and mass matrices are very small compared to the diagonal elements and can be neglected.

It is obvious that the usage of the physically consistent non-orthogonal natural modes in the modal superposition method is not practical, especially not in the case of time integration. Therefore, it is preferable to use mathematical orthogonal modes for that purpose. They are created by the static displacement relations yielding from Eqs. (24) and (25) with $$\omega =0$$ , that leads to $$\alpha =\eta =1$$ . In that case all finite element matrices, defined with Eqs. (37) and in Appendix A, can be transformed into explicit form, Appendix B.

9. Cross-section properties of thin-walled girder

Geometrical properties of a thin-walled girder include cross-section area A, moment of inertia of cross-section I _b , shear area A _s , torsional modulus I _t , warping modulus I _w and shear inertia modulus I _s . These parameters are determined analytically for a simple cross-section as pure geometrical properties (Haslum & Tonnessen, 1972, Pavazza, 1991, 2005, Vlasov, 1961).

However, determination of cross-section properties for an open multi-cell cross-section, as for instance in case of ship structures, is quite a difficult task. Therefore, the strip element method is applied for solving this statically indetermined problem (Cheung, 1976). That is well-known and widely used theory of thin-walled girders, which is only briefly described here. Firstly, axial node displacements are calculated due to bending caused by shear force, and due to torsion caused by variation of twist angle. Then, shear stress in bending τ _b , shear stress due to pure torsion τ _t , shear and normal stresses due to restrained warping τ _w and σ _w , respectively, are determined. Based on the equivalence of strain energies induced by sectional forces and calculated stresses, it is possible to specify cross-section properties in the same formulation as presented below. Furthermore, those formulae can be expressed by stress flows, i.e. stresses due to unit sectional forces (Senjanović & Fan, 1992, 1993).

Shear area:

$$A_s=\frac{Q^2}{{\displaystyle \underset{A}{\int }{\tau }_b^{2}dA}}=\frac{1}{{\displaystyle \underset{A}{\int }{g}_b^{2}dA}},g_b=\frac{{\tau }_b}{Q}$$ E77

.

Torsional modulus:

$$I_t=\frac{T_t^2}{{\displaystyle \underset{A}{\int }{\tau }_t^{2}dA}}=\frac{1}{{\displaystyle \underset{A}{\int }{g}_t^{2}dA}},g_t=\frac{{\tau }_t}{T_t}$$ E78

.

Shear inertia modulus:

$$I_s=\frac{T_w^2}{{\displaystyle \underset{A}{\int }{\tau }_w^{2}dA}}=\frac{1}{{\displaystyle \underset{A}{\int }{g}_w^{2}dA}},g_w=\frac{{\tau }_w}{T_w}$$ E79

.

Warping modulus:

$$I_w=\frac{B_w^2}{{\displaystyle \underset{A}{\int }{\sigma }_w^{2}dA}}=\frac{1}{{\displaystyle \underset{A}{\int }{f}_w^{2}dA}},f_w=\frac{{\sigma }_w}{B_w};I_w={\displaystyle \underset{A}{\int }{w}^{2}dA}$$ E80

.

The above quantities are not pure geometrical cross-section properties any more, since they also depend on Poisson's ratio as a physical parameter.

The mass parameters can be expressed with the given mass distribution per unit length, m, and calculated cross-section parameters, i.e.

$$J_b=\frac{m}{A}{I}_b,J_t^0=\frac{m}{A}{I}_p,J_w=\frac{m}{A}{I}_w$$ E81

.

where $$I_p=I_{by}+I_{bz}$$ is the polar moment of inertia of cross-section.

10. Illustrative numerical examples

For the illustration of the procedure related to engine room effective stiffness determination, 3D FEM analysis of ship-like pontoon has been undertaken. The 3D FEM model is constituted according to 7800 TEU container ship with main dimensions $$L_{pp}xBxH=319x42.8x24.6$$ m, and detailed desciption given in (Tomašević, 2007). The complete hydroelastic analysis of the same ship has been performed.

Stiffness properties of ship hull are calculated by program STIFF, based on the theory of thin-walled girders (STIFF, 1990), Fig. 11.

Figure 11.

Program STIFF – warping of ship cross-section

Influence of the transverse bulkheads is taken into account by using the equivalent torsional modulus for the open cross-sections instead of the actual values, i.e. $$I_t^{*}=2.4I_t$$ . This value is applied for all ship-cross sections as the first approximation.

10.1. Analysis of ship-like segmented pontoon

Torsion of the segmented pontoon of the length L = 300 m, with effective parameters is considered. Torsional moment M _t = 40570 kNm is imposed at the pontoon ends. The pontoon is considered free in the space and the problem is solved analytically according to the formulae given in Section 4. The following values of the basic parameters are used: $$a=10.1$$ m, $$b=19.17$$ m, $$t_1=0.01645$$ m, $$w_D=-221$$ m², $$w_B=267$$ m², $$I_t^{\circ }=14.45$$ m⁴, $$k=1.894$$ . As a result $$C=22.42$$ , Eq. (59), and accordingly $${\tilde{I}}_t^{}=338.4$$ m⁴, Eq. (58a), are obtained. Since $${\tilde{I}}_t^{}=0.36I_t^{\ast }$$ , effect of the short engine room structure on its torsional stiffness is obvious.

Figure 12.

Deformation of segmented pontoon, lateral and bird view

Figure 13.

Lateral, axial, bird and fish views on deformed engine room superelement

Figure 14.

Twist angles of segmented pontoon

The 3D FEM model of segmented pontoon is made by commercial software package SESAM and consists of 20 open and 1 closed (engine room) superelement. The pontoon ends are closed with transverse bulkheads. The shell finite elements are used. The pontoons are loaded at their ends with the vertical distributed forces in the opposite directions, generating total torque M _t = 40570 kNm. The midship section is fixed against transverse and vertical displacements, and the pontoon ends are constrained against axial displacements (warping). Lateral and bird view on the deformed segmented pontoon is shown in Fig. 12, where the influence of more rigid engine room structure is evident. Detailed view on this pontoon portion is presented in Fig. 13. It is apparent that segment of very stiff double bottom and sides rotate as a “rigid body”, while decks and transverse bulkheads are exposed to shear deformation. This deformation causes the distortion of the cross-section, Fig. 13.

Twist angles of the analytical beam solution and that of 3D FEM analysis for the pontoon bottom are compared in Fig. 14. As it can be noticed, there are some small discrepancies between $${\psi }_{\left(1+2\right)D}$$ and $${\psi }_{3D,bottom}$$ , which are reduced to a negligible value at the pontoon ends

Fig. 14 also shows twist angle of side structure and the difference $$\delta ={\psi }_{3\text{D}\text{,bottom}}-{\psi }_{3\text{D}\text{,side}}$$ represents distortion angle of cross-section which is highly pronounced. As it is mentioned before, the problem will be further investigated.

10.2. Validation of 1D FEM model

The reliability of 1D FEM analysis is verified by 3D FEM analysis of the considered ship. For this purpose, the light weight loading condition of dry ship with displacement Δ=33692 t is taken into account. The equivalent torsional stiffness of the engine room structure, as well as equivalent stiffness of fore and aft peaks is not taken into account in this example for the time being. However, it will be done in the next step of investigation. The lateral and bird view of the first dominantly torsional and second dominantly horizontal mode of the wetted surface, determined by 1D model, is shown in Fig. 15.

.

Figure 15.

The first and second mode, lateral and bird view, light weight, 1D modelTwist angles of segmented pontoon

The first and second 3D dry coupled natural modes of the complete ship structure are shown in Fig. 16. They are similar to that of 1D analysis for the wetted surface. Warping of the transverse bulkheads, which increases the hull torsional stiffness, is evident.

The first four corresponding natural frequencies obtained by 1D and 3D analyses are compared in Table 1.

Mode no.	Vert.		Horiz. + tors.		Mode no.
1D	3D	1D	3D
1	7.35	7.33	4.17	4.15	1(H0 + T1)
2	15.00	14.95	7.34	7.40	2(H1 + T2)
3	24.04	22.99	12.22	12.09	3(H2 + T3)
4	35.08	34.21	15.02	16.22	4(H3 + T4)

Table 1.

Dry natural frequencies, light weight, ω _i [rad/s]

Figure 16.

The first and second mode, lateral and bird view, light weight, 3D modelThe first and second mode, lateral and bird view, light weight, 1D modelTwist angles of segmented pontoon

Quite good agreement is achieved. Values of natural frequencies for higher modes are more difficult to correlate, since strong coupling between global hull modes and local substructure modes of 3D analysis occurs.

10.3. Hydroelastic response of large container ship

Transfer functions of torsional moment and horizontal bending moment at the midship section, obtained using 1D structural model, are shown in Figs. 17 and 18, respectively. The angle of 180° is related to head sea. They are compared to the rigid body ones determined by program HYDROSTAR. Very good agreement is obtained in the lower frequency domain, where the ship behaves as a rigid body, while large discrepancies occur at the resonances of the elastic modes, as expected.

Figure 17.

Transfer function of torsional moment, χ=120°, U=25 kn,

Figure 18.

Transfer function of horizontal bending moment, χ=120°, U=25 kn, x=155.75 m from APTransfer function of torsional moment, χ=120°, U=25 kn,

11. Conclusion

Ultra large container ships are quite elastic and especially sensitive to torsion due to large deck openings. The wave induced response of such ships should be determined by using mathematical hydroelastic models which are consisted of structural, hydrostatic and hydrodynamic parts.

In this chapter the methodology of ship hydroelastic analysis is briefly described, and the role of structural model is discussed. After that, full detail description of the sophisticated beam structural model, which takes shear influence on torsion, as well as contribution of transverse bulkheads and engine room structure to the hull stiffness, is given. Numerical procedure for vibration analysis is also described and determination of ship cross-section properties is explained. The developed theories are illustrated through the numerical examples which include analysis of torsional response of a ship-like segmented pontoon, free vibration analysis of a large container ship and comparison with the results obtained using 3D FEM model, and complete global hydroelastic analysis of a container ship.

It is shown that the used sophisticated beam model of ship hull, based on the advanced thin-walled girder theory with included shear influence on torsion and a proper contribution of transverse bulkheads and engine room structure to its stiffness, is a reasonable choice for determining wave load effects. However, based on the experience, stress concentration in hatch corners calculated directly by the beam model is underestimated. This problem can be overcome by applying substructure approach, i.e. 3D FEM model of substructure with imposed boundary conditions from beam response. In any case, 3D FEM model of complete ship is preferable from the viewpoint of determining stress concentration. Concerning further improvements of the beam model, the distortion induced by torsion is of interest.

The illustrative numerical example of the 7800 TEU container ship shows that the developed hydroelasticity theory, utilizing sophisticated 1D FEM structural model and 3D hydrodynamic model, is an efficient tool for application in ship hydroelastic analyses. The obtained results point out that the transfer functions of hull sectional forces in case of resonant vibration (springing) are much higher than in resonant ship motion.

13. Appendix A – consistent finite element properties (frequency dependent formulation)

The stiffness and mass matrices, Eqs. (37), are expressed with one or two integrals, which can be classified in three different types. For general notation of shape functions

$$g_i=⟨g_{ik}⟩\left\{{\xi }^k\right\},i=\mathrm{1,2,3,4};k=\mathrm{0,1,2,3}$$ E82

where $$q_{ik}$$ are coefficients and $$\xi =x/l$$ , one finds the solutions of integrals in the following form:

$$\begin{array}{l}{I}_{ij}^{\left(0\right)}\left(g_{ik},g_{jk}\right)={\displaystyle \underset{0}{\overset{l}{\int }}{g}_i{g}_j\text{d}x}=⟨g_{ik}⟩{\displaystyle \underset{0}{\overset{l}{\int }}\left\{{\xi }^k\right\}⟨{\xi }^k⟩\text{d}x}\left\{g_{jk}\right\}\\ =l[g_{i0}{g}_{j0}+\frac{1}{2}\left(g_{i1}{g}_{j0}+g_{i0}{g}_{j1}\right)+\frac{1}{3}\left(g_{i2}{g}_{j0}+g_{i1}{g}_{j1}+g_{i0}{g}_{j2}\right)\\ \text{+}\frac{\text{1}}{\text{4}}\left(g_{i0}{g}_{j3}+g_{i1}{g}_{j2}+g_{i2}{g}_{j1}+g_{i3}{g}_{j0}\right)+\frac{1}{5}\left(g_{i1}{g}_{j3}+g_{i2}{g}_{j2}+g_{i3}{g}_{j1}\right)\\ \text{+}\frac{\text{1}}{\text{6}}\left(g_{i2}{g}_{j3}+g_{i3}{g}_{j2}\right)+\frac{1}{7}{g}_{i3}{g}_{j3}]\end{array}$$ E83

$$\begin{array}{l}{I}_{ij}^{\left(1\right)}\left(g_{ik},g_{jk}\right)={\displaystyle \underset{0}{\overset{l}{\int }}\frac{\text{d}{g}_i}{\text{d}x}\frac{\text{d}{g}_j}{\text{d}x}\text{d}x}=⟨g_{ik}⟩{\displaystyle \underset{0}{\overset{l}{\int }}\left\{k{\xi }^{k-1}\right\}⟨k{\xi }^{\text{k-1}}⟩\text{d}x}\left\{g_{jk}\right\}\\ =\frac{1}{l}[g_{i1}{g}_{j1}+g_{i1}{g}_{j2}+g_{i2}{g}_{j1}+g_{i1}{g}_{j3}+g_{i3}{g}_{j1}\\ \text{+}\frac{\text{4}}{\text{3}}{g}_{i2}{g}_{j2}+\frac{3}{2}\left(g_{i2}{g}_{j3}+g_{i3}{g}_{j2}\right)+\frac{9}{5}{g}_{i3}{g}_{j3}]\end{array}$$ E84

$$\begin{array}{l}{I}_{ij}^{\left(2\right)}\left(g_{ik},g_{jk}\right)={\displaystyle \underset{0}{\overset{l}{\int }}\frac{{\text{d}}^2{g}_i}{\text{d}{x}^2}\frac{{\text{d}}^{\text{2}}{g}_j}{\text{d}{x}^2}\text{d}x}=⟨g_{ik}⟩{\displaystyle \underset{0}{\overset{l}{\int }}\left\{k\left(k-1\right){\xi }^{k-2}\right\}⟨k\left(k-1\right){\xi }^{k-2}⟩\text{d}x}\left\{g_{jk}\right\}\\ =\frac{4}{l^3}\left[g_{i2}{g}_{j2}+\frac{3}{2}\left(g_{i2}{g}_{j3}+g_{i3}{g}_{j2}\right)+3g_{i3}{g}_{j3}\right].\end{array}$$ E85

Thus, the finite element properties can be written in the following systematic way suitable for coding.Stiffness matrices

$$\begin{array}{l}{\left[k\right]}_{bs}=E{I}_b\left[I_{ij}^{\left(2\right)}\left(a_{ik},a_{jk}\right)\right]+G{A}_s\left[I_{ij}^{\left(1\right)}\left(b_{ik},b_{jk}\right)\right]\\ {\left[k\right]}_{ws}=E{I}_w\left[I_{ij}^{\left(2\right)}\left(d_{ik},d_{jk}\right)\right]+G{I}_s\left[I_{ij}^{\left(1\right)}\left(e_{ik},e_{jk}\right)\right]\\ {\left[k\right]}_t=G{I}_t\left[I_{ij}^{\left(1\right)}\left(d_{ik},d_{jk}\right)\right]\end{array}$$ E86

Mass matrices

$$\begin{array}{l}{\left[m\right]}_{sb}=m\left[I_{ij}^{\left(0\right)}\left(c_{ik},c_{jk}\right)\right]+J_b\left[I_{ij}^{\left(1\right)}\left(a_{ik},a_{jk}\right)\right]\\ {\left[m\right]}_{tw}=J_t\left[I_{ij}^{\left(0\right)}\left(f_{ik},f_{jk}\right)\right]+J_w\left[I_{ij}^{\left(1\right)}\left(d_{ik},d_{jk}\right)\right]\\ {\left[m\right]}_{st}=m_c\left[I_{ij}^{\left(0\right)}\left(c_{ik},f_{jk}\right)\right],{\left[m\right]}_{ts}={\left[m\right]}_{st}^T\end{array}$$ E87

Load vectors

$$\begin{array}{l}\left\{q\right\}=l\left\{q_0\left(c_{i0}+\frac{1}{2}{c}_{i1}+\frac{1}{3}{c}_{i2}+\frac{1}{4}{c}_{i3}\right)+q_1\left(\frac{1}{2}{c}_{i0}+\frac{1}{3}{c}_{i1}+\frac{1}{4}{c}_{i2}+\frac{1}{5}{c}_{i3}\right)\right\}\\ \left\{\mu \right\}=l\left\{{\mu }_0\left(f_{i0}+\frac{1}{2}{f}_{i1}+\frac{1}{3}{f}_{i2}+\frac{1}{4}{f}_{i3}\right)+{\mu }_1\left(\frac{1}{2}{f}_{i0}+\frac{1}{3}{f}_{i1}+\frac{1}{4}{f}_{i2}+\frac{1}{5}{f}_{i3}\right)\right\}\end{array}$$ E88

14. Appendix B – simplified finite element properties, from appendix A (frequency independent formulation)

Stiffness matrices:

$${\left[k\right]}_{bs}=\frac{2E{I}_b}{\left(1+12\beta \right)l^3}\left[\begin{array}{cccc}6& 3l& -6& 3l\\ & 2\left(1+3\beta \right)l^2& -3l& \left(1-6\beta \right)l^2\\ & & 6& -3l\\ Sym.& & & 2\left(1+3\beta \right)l^2\end{array}\right]$$ E89

$${\left[k\right]}_{ws}=\frac{2E{I}_w}{\left(1+12\gamma \right)l^3}\left[\begin{array}{cccc}6& 3l& -6& 3l\\ & 2\left(1+3\gamma \right)l^2& -3l& \left(1-6\gamma \right)l^2\\ & & 6& -3l\\ Sym.& & & 2\left(1+3\gamma \right)l^2\end{array}\right]$$ E90

$${\left[k\right]}_t=\frac{G{I}_t}{30{\left(1+12\gamma \right)}^{2}l}\left[\begin{array}{cccc}36& 3\left(1-60\gamma \right)l& -36& 3\left(1-60\gamma \right)l\\ & 4\left(1+15\gamma +360{\gamma }^2\right)l^2& -3\left(1-60\gamma \right)l& -\left(1+60\gamma -720{\gamma }^2\right)l^2\\ & & 36& -3\left(1-60\gamma \right)l\\ Sym.& {}^{}& & 4\left(1+15\gamma +360{\gamma }^2\right)l^2\end{array}\right]$$ E91

Mass matrices:

$${\left[m\right]}_{sb}={\left[m\right]}_s+{\left[m\right]}_b$$ E92

$${\left[m\right]}_s=\frac{ml}{420{\left(1+12\beta \right)}^2}\left[\begin{array}{cccc}156+3528\beta +20160{\beta }^2& \left(22+462\beta +2520{\beta }^2\right)l& 54+1512\beta +10080{\beta }^2& -\left(13+378\beta +2520{\beta }^2\right)l\\ & \left(4+84\beta +504{\beta }^2\right)l^2& \left(13+378\beta +2520{\beta }^2\right)l& -\left(3+84\beta +504{\beta }^2\right)l^2\\ & & 156+3528\beta +20160{\beta }^2& -\left(22+462\beta +2520{\beta }^2\right)l\\ Sym.& & & \left(4+84\beta +504{\beta }^2\right)l^2\end{array}\right]$$ E93

$${\left[m\right]}_b=\frac{J_b}{30{\left(1+12\beta \right)}^{2}l}\left[\begin{array}{cccc}36& \left(3-180\beta \right)l& -36& \left(3-180\beta \right)l\\ & \left(4+60\beta +1440{\beta }^2\right)l^2& \left(-3+180\beta \right)l& -\left(1+60\beta -720{\beta }^2\right)l^2\\ & & 36& \left(-3+180\beta \right)l\\ Sym.& & & \left(4+60\beta +1440{\beta }^2\right)l^2\end{array}\right]$$ E94

$${\left[m\right]}_{tw}={\left[m\right]}_t+{\left[m\right]}_w$$ E95

$${\left[m\right]}_t=\frac{J_{t}l}{420{\left(1+12\gamma \right)}^2}\left[\begin{array}{cccc}156+3528\gamma +20160{\gamma }^2& \left(22+462\gamma +2520{\gamma }^2\right)l& 54+1512\gamma +10080{\gamma }^2& -\left(13+378\gamma +2520{\gamma }^2\right)l\\ & \left(4+84\gamma +504{\gamma }^2\right)l^2& \left(13+378\gamma +2520{\gamma }^2\right)l& -\left(3+84\gamma +504{\gamma }^2\right)l^2\\ & & 156+3528\gamma +20160{\gamma }^2& -\left(22+462\gamma +2520{\gamma }^2\right)l\\ Sym.& & & \left(4+84\gamma +504{\gamma }^2\right)l^2\end{array}\right]$$ E96

$${\left[m\right]}_w=\frac{J_w}{30{\left(1+12\gamma \right)}^{2}l}\left[\begin{array}{cccc}36& \left(3-180\gamma \right)l& -36& \left(3-180\gamma \right)l\\ & \left(4+60\gamma +1440{\gamma }^2\right)l^2& \left(-3+180\gamma \right)l& -\left(1+60\gamma -720{\gamma }^2\right)l^2\\ & & 36& \left(-3+180\gamma \right)l\\ Sym.& & & \left(4+60\gamma +1440{\gamma }^2\right)l^2\end{array}\right]$$ E97

$$\begin{array}{l}{\left[m\right]}_{st}=\frac{mlc}{420\left(1+12\beta \right)\left(1+12\gamma \right)}\times \\ \left[\begin{array}{cccc}156+1764\beta +1764\gamma +20160\beta \gamma & \left(22+252\beta +210\gamma +2520\beta \gamma \right)l& 54+756\beta +756\gamma +10080\beta \gamma & -\left(13+168\beta +210\gamma +2520\beta \gamma \right)l\\ \left(22+210\beta +252\gamma +2520\beta \gamma \right)l& \left(4+42\beta +42\gamma +504\beta \gamma \right)l^2& \left(13+210\beta +168\gamma +2520\beta \gamma \right)l& -\left(3+42\beta +42\gamma +504\beta \gamma \right)l^2\\ 54+756\beta +756\gamma +10080\beta \gamma & \left(13+168\beta +210\gamma +2520\beta \gamma \right)l& 156+1764\beta +1764\gamma +20160\beta \gamma & -\left(22+252\beta +210\gamma +2520\beta \gamma \right)l\\ -\left(13+210\beta +168\gamma +2520\beta \gamma \right)l& -\left(3+42\beta +42\gamma +504\beta \gamma \right)l^2& -\left(22+210\beta +252\gamma +2520\beta \gamma \right)l& \left(4+42\beta +42\gamma +504\beta \gamma \right)l^2\end{array}\right]\end{array}$$ E98

$${\left[m\right]}_{ts}={\left[m\right]}_{st}^T$$ E99

Load vectors:

$$\left\{q\right\}=\frac{q_{0}l}{12}\left\{\begin{array}{c}6\\ l\\ 6\\ -l\end{array}\right\}+\frac{q_{1}l}{60\left(1+12\beta \right)}\left\{\begin{array}{c}9+120\beta \\ \left(2+30\beta \right)l\\ 21+240\beta \\ -\left(3+30\beta \right)l\end{array}\right\}$$ E100

$$\left\{\mu \right\}=\frac{{\mu }_{0}l}{12}\left\{\begin{array}{c}6\\ l\\ 6\\ -l\end{array}\right\}+\frac{{\mu }_{1}l}{60\left(1+12\gamma \right)}\left\{\begin{array}{c}9+120\gamma \\ \left(2+30\gamma \right)l\\ 21+240\gamma \\ -\left(3+30\gamma \right)l\end{array}\right\}$$ E101

Stiffness ratios:

$$\beta =\frac{E{I}_b}{G{A}_s{l}^2},\gamma =\frac{E{I}_w}{G{I}_s{l}^2}.$$ E102

Acknowledgments

This investigation is carried out within the EU FP7 Project TULCS (Tools for Ultra Large Container Ships) and the project of Croatian Ministry of Science, Education and Sports Load and Response of Ship Structures.