1 Introduction

Since their first use, composites have achieved properties that single materials cannot obtain. Combining two or more materials enhances its performance, leading to stiffness and strength increase and weight reduction. This potential improvement makes composites an aerospace-suitable material, as weight-saving and superior strength are crucial factors for an aircraft’s structure. However, these materials do not present the same behaviour once an impact occurs, although they generally do not show visible delamination or other surface damages. This event enormously varies their damage tolerance and resistance, not being noticeable. Among other reasons, stiffeners usage aims to increase composite shell's impact resistance and improve its dynamic behaviour. Even though several pieces of research have already been made to analytically model the low-velocity impact behaviour of composite plates, i.e., scenarios with impact energy ranging up to 30 J, there are no significant previous studies proposing models for the behaviour of stringer-reinforced composite panels under those conditions.

Therefore, this study seeks to develop an analytical model that can predict the dynamic behaviour of the super stringers assembly. The most representative studies concerning composite plate impact modelling are presented below.

Song et al. [1] studied failure and dynamic behaviour of stringer-stiffened curved composite panels under hail ice impact at two distinct locations within the panel. The damage phenomenon was investigated through experimentally verified finite element methods, and the dynamic response unveiled a stress wave propagation. Numerical simulations brought precise predictions of failure mode and its propagation, mainly delamination, the total delamination area against impact velocity, and the debonding between skin and stringer flange.

Rout and Karmakar [2] also developed finite element models to study the low-velocity impact response of laminated stiffened cylindrical shells. The formulation was based on the first-order shear deformation theory. The modified Hertzian contact law and Newmark’s time integration algorithm were used to compute the contact force and solve the time-dependent equations. Within this study, some conclusions emerged, such as the dependency of contact force magnitude to the stiffener and impact velocity and the effect of delamination and its location in displacement histories.

Sun et al. [3] undertook experimental studies of the effect of stiffener damage caused by low-velocity impact on compressive buckling and failure modes of T-stiffened composite panels. Tests showed that the stiffener damage and stiffener-panel debonding got severe as impact energy rose, whereas panel damage was minor than Barely Visible Impact Damage (BVID). Experimental results from the compression after impact (CAI) tests revealed that due to damage propagation after impact, compressive stiffness exhibited a decreasing tendency with the increase of compressive load, with a maximum of 44% reduction compared to the original panel.

Sun et al. [4] performed experiments and finite element modelling of low-velocity impacts against metallic and composite foam-cored sandwich panels. Different combinations of face sheets materials and core density distribution showed that the contact force history's dominant plate is the front skin. If used, stepwise graded core densities affect the contact force history and energy absorption capacity, being this effect found to diminish as the impact energy increases. The increase in the front core layer's density increased the maximum contact force while reducing the indentation magnitude.

Pernas-Sánchez et al. [5] researched ice impact modelling on steel and composite panels. Numerical methodology and experimental validation were first developed to understand soft body impactors' behaviour, focusing on the projectile model capability to transfer the dynamic load's effects to rigid steel panels. The validated ice material model was then used in conjunction with a simplified composite model to reproduce intralaminar and interlaminar failure of laminate plates, which proved viable to predict delaminated areas.

Gong and Lam [6] developed an analytical method to predict low-velocity impact behaviour on stringer stiffened composite plates. The contact force was predicted using the linearised Hertzian contact law, while the dynamic behaviour was modelled through the spring-mass approach on a comparable property flat panel. Also, FE models were developed using LS-DYNA, in which the plate and the stringer were explicitly modelled, while the interaction between the striker and the stiffened plate was obtained though Hertzian contact law. Both approaches were reasonably accurate with respect to the available literature.

Khalili et al. [7] studied the dynamic response of a thin smart curved composite panel subjected to a low-velocity transverse impact embedded with shape-memory wires. Their work was based on the linear Hertzian contact model, which is linearised for the impact analysis of the curved composite panel. The curved panel's governing equations are provided by the first-order shear deformation theory and solved by the Fourier series. Sing and Mahajan [8] developed an improved analytical model to predict the response of composite laminate under low velocity large mass impact is proposed. In their work, a spring-mass system was used to represent the contact, bending, shear and membrane stiffnesses of the laminate structure. The stiffness of the springs was evaluated by dividing the laminate into a damaged central region with degraded stiffness surrounded by an undamaged area.

Schiffer et al. [9] investigated analytical models to predict the transient elastic response of fully clamped circular composite plates subject to high-velocity impact by a rigid spherical projectile. The models are based on first-order shear deformation plate theory and account for the effects of large deformations as well as propagation and reflection of flexural waves. Analytical predictions of plate deflection history and peak strain in the plates were found in good agreement with those obtained from detailed explicit FE simulations.

Gregori et al. [10] studied analytical and numerical models to simulate the perforation of ceramic-composite targets by small-calibre projectiles. The modified Bernoulli equation has been implemented in a new analytical model to simulate the interaction between projectile and ceramic tile. An energy formulation based on the wave propagation theory has been adopted to describe the composite backing's energy absorption.

Wagner et al. [11] developed a simplified empirical approach to describe a composite bolted joint's overall behaviour for explicit crash and impact simulations. Based on analytical equations, classical laminate theory and empirical factors, a trilinear force–displacement curve was proposed for bolted joints under in-plane loads. The empirical factors were derived from single-lap shear experiments of carbon fibre prepreg specimens (Hexcel M18-1/G939) with different fasteners.

Sharyat and Roshanfar [12] developed an analytical solution based on a new idea of superposition of two kinematic descriptions is presented for dynamic response analysis of a multi-layer/sandwich composite plate with point supports subjected to an eccentric low-velocity impact. Direct and virtual-work-based novel energy formulations were proposed for the problem that considers the indentation region's potential energy. The non-linear governing equations of motions are found based on minimisation of the total potential energy of the whole mass-plate system, including the inertia forces, employing the Ritz technique and transformation of the time-dependent non-linear system of governing equations to a non-linear algebraic one through a novel concept.

Arachchige and Ghasemnejad [13] developed an analytical approach to study the dynamic behaviour of curved variable-thickness composite plates under low-velocity impact conditions. In this research, first-order shear deformation theory and linearised Hertzian contact law were used, whereas the dynamic behaviour was obtained as a Spring-Mass model. The model was first validated with the available literature, leading to accurate results, and then used for analysing different laminate and thickness layout combinations.

Saleh and Soutis [14] reviewed the knowledge gap in the literature of 3D woven composites, suggesting opportunities for future research in this field and room for improvement in utilising Non-Destructive Techniques (NDT), such as Digital Image Correlation (DIC), Acoustic Emission (AE) and X-ray Computed Tomography (CT), for observing damage initiation and evolution in 3D woven composites that could be used to calibrate and evaluate analytical and numerical models.

However, none of this research has included details of structural components such as a change of thickness or stiffness in the analytical models. In this regard, there is still a gap in developing an advanced analytical model to represent a real composite structure, including details of Omega stringers and variable stiffness against impact loading. The proposed method in this paper uses an analytical force function based on material properties, the mass of the shell/striker and the impact velocity to predict the impact force response. This force function will be used to depict contact-force histories for different cases of impactor masses and velocity. This force function is used to represent contact-force histories for different cases of impactor masses and velocity. The models will be based on the Hertzian contact model, which is linearised for the composite panel's impact analysis with simply supported, clamped boundary conditions and a combination of them. This function can be converted to force–displacement to measure the amount of absorbed energy.

2 Analytical Impact Modelling

As a result of composites anisotropy and large Young's modulus to shear modulus ratio, the shear effects throughout the thickness have to be considered. Therefore, plate behaviour will be modelled according to the first-order shear deformation theory (FSDT) in this study. Then, the governing displacement field can be expressed as in Eqs. (1−3).

$$\mathrm{u}\left(\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{t}\right)={\mathrm{u}}^{0}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)+\mathrm{z}{\uppsi }_{\mathrm{x}}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)$$
(1)
$$\mathrm{v}\left(\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{t}\right)={\mathrm{v}}^{0}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)+\mathrm{z}{\uppsi }_{\mathrm{y}}(\mathrm{x},\mathrm{y},\mathrm{t})$$
(2)
$$\mathrm{w}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)={\mathrm{w}}^{0}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)$$
(3)

For a simply supported rectangular plate of dimensions \(\mathrm{a}\times \mathrm{b}\), a model was developed a system of equations that defines [12, 13], based on the FSDT, the motion field \(\left(\mathrm{u},\mathrm{v},\mathrm{w},{\uppsi }_{\mathrm{x}},{\uppsi }_{\mathrm{y}}\right)\) of the plate. The assumption is based on a perfect bonding between stringer and skin, and there is no debonding at interfaces. Besides, for specially orthotropic laminates (i.e., \({\mathrm{A}}_{16}={\mathrm{A}}_{26}={\mathrm{D}}_{16}={\mathrm{D}}_{26}={\mathrm{B}}_{\mathrm{ij}}=0\)), and by neglecting in-plane accelerations and rotary inertias, these plate equations are reduced to Eqs. (4, 5).

$${\mathrm{D}}_{11}\frac{{\partial }^{2}{\uppsi }_{\mathrm{x}}}{\partial {\mathrm{x}}^{2}}+{\mathrm{D}}_{66}\frac{{\partial }^{2}{\uppsi }_{\mathrm{x}}}{\partial {\mathrm{y}}^{2}}+\left({\mathrm{D}}_{12}+{\mathrm{D}}_{66}\right)\frac{{\partial }^{2}{\uppsi }_{\mathrm{y}}}{\partial \mathrm{x}\partial \mathrm{y }}-{\mathrm{k}}_{\mathrm{sh}}{\mathrm{A}}_{55}\left({\uppsi }_{\mathrm{x}}+\frac{\partial \mathrm{w}}{\partial \mathrm{x}}\right)=0$$
(4)
$$\left({\mathrm{D}}_{12}+{\mathrm{D}}_{66}\right)\frac{{\partial }^{2}{\uppsi }_{\mathrm{x}}}{\partial \mathrm{x}\partial \mathrm{y}}+{\mathrm{D}}_{66}\frac{{\partial }^{2}{\uppsi }_{\mathrm{y}}}{\partial {\mathrm{x}}^{2}}+{\mathrm{D}}_{22}\frac{{\partial }^{2}{\uppsi }_{\mathrm{y}}}{\partial {\mathrm{y}}^{2}}-{\mathrm{k}}_{\mathrm{sh}}{\mathrm{A}}_{44}\left({\uppsi }_{\mathrm{y}}+\frac{\partial \mathrm{w}}{\partial \mathrm{y }}\right)=0$$
(5)
$${\mathrm{k}}_{\mathrm{sh}}{\mathrm{A}}_{55}\left(\frac{\partial {\psi }_{\mathrm{x}}}{\partial \mathrm{x}}+\frac{{\partial }^{2}\mathrm{w}}{\partial {\mathrm{x}}^{2}}\right)+{\mathrm{k}}_{\mathrm{sh}}{\mathrm{A}}_{44}\left(\frac{\partial {\psi }_{\mathrm{y}}}{\partial \mathrm{y}}+\frac{{\partial }^{2}\mathrm{w}}{\partial {\mathrm{y}}^{2}}\right)+\mathrm{q}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)=\rho\mathrm{h}\frac{{\partial }^{2}\mathrm{w}}{\partial {\mathrm{t}}^{2}}$$
(6)

where \(\mathrm{h}\) is the plate thickness and \({\mathrm{k}}_{\mathrm{sh}}\) is the shear correction factor commonly assumed to be \({\uppi }^{2}/12\) [12]. Within the framework of this study, some non-specially orthotropic laminates are analysed through the use of Eqs. (4, 5), which may be done provided that the components \({\mathrm{D}}_{16}\) and \({\mathrm{D}}_{26}\) of the laminate, \(\mathrm{D}\) matrix is at least an order of magnitude smaller than the other \({\mathrm{D}}_{\mathrm{ij}}\) coefficients. The system of differential Eqs. (4, 5), along with the boundary conditions in Eqs. (7, 8), define the solution of the vertical displacement \(\mathrm{w}\) and shear rotations, \({\uppsi }_{\mathrm{x}}\) and \({\uppsi }_{\mathrm{y}}\), for a simply supported rectangular plate of dimensions \(\mathrm{a x b}\) under the load case \(\mathrm{q}\).

$$\mathrm{w}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)={\mathrm{M}}_{\mathrm{x}}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)=0\mathrm{ for x}=0 \ \&\mathrm{ x}=\mathrm{a}$$
(7)
$$\mathrm{w}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)={\mathrm{M}}_{\mathrm{x}}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)=0\mathrm{ for x}=0 \ \&\mathrm{ x}=\mathrm{b}$$
(8)

Given Eqs. (4, 5, 6, 7), the motion of the plate when a known load \(\mathrm{q}\) is applied can be obtained using a double expansion in Fourier series of \(\mathrm{w}\), \({\uppsi }_{\mathrm{x}}\) and \({\uppsi }_{\mathrm{y}}\), and the load \(\mathrm{q}\) as well as per Eqs. (9, 10).

$${\psi }_{\mathrm{x}}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)={\sum }_{\mathrm{m}}{\sum }_{\mathrm{n}}{\mathrm{X}}_{\mathrm{mn}}\left(\mathrm{t}\right)\mathrm{cos}\left(\frac{\mathrm{m}\pi \mathrm{x}}{\mathrm{a}}\right)\mathrm{sin}\left(\frac{\mathrm{n}\pi \mathrm{y}}{\mathrm{b}}\right)$$
(9)
$${\psi }_{\mathrm{y}}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)={\sum }_{\mathrm{m}}{\sum }_{\mathrm{n}}{\mathrm{Y}}_{\mathrm{mn}}\left(\mathrm{t}\right)\mathrm{sin}\left(\frac{\mathrm{m}\pi \mathrm{x}}{\mathrm{a}}\right)\mathrm{cos}\left(\frac{\mathrm{n}\pi \mathrm{y}}{\mathrm{b}}\right)$$
(10)
$$\mathrm{w}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)={\sum }_{\mathrm{m}}{\sum }_{\mathrm{n}}{\mathrm{W}}_{\mathrm{mn}}\left(\mathrm{t}\right)\mathrm{sin}\left(\frac{\mathrm{m}\pi \mathrm{x}}{\mathrm{a}}\right)\mathrm{sin}\left(\frac{\mathrm{n}\pi \mathrm{y}}{\mathrm{b}}\right)$$
(11)
$$\mathrm{q}\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)={\sum }_{\mathrm{m}}{\sum }_{\mathrm{n}}{\mathrm{Q}}_{\mathrm{mn}}\left(\mathrm{t}\right)\mathrm{sin}\left(\frac{\mathrm{m}\pi \mathrm{x}}{\mathrm{a}}\right)\mathrm{sin}\left(\frac{\mathrm{n}\pi \mathrm{y}}{\mathrm{b}}\right)$$
(12)

For the case of a concentrated load at a point defined at (\({\mathrm{x}}_{\mathrm{c}},{\mathrm{y}}_{\mathrm{c}}\)), \({\mathrm{Q}}_{\mathrm{mn}}(\mathrm{t})\) can be expressed as:

$${\mathrm{Q}}_{\mathrm{mn}}\left(\mathrm{t}\right)=\frac{4\mathrm{F}(\mathrm{t})}{\mathrm{ab}}\mathrm{sin}\left(\frac{\mathrm{m}\pi }{\mathrm{a}}{\mathrm{x}}_{\mathrm{c}}\right)\mathrm{sin}\left(\frac{\mathrm{n}\pi }{\mathrm{b}}{\mathrm{y}}_{\mathrm{c}}\right)$$
(13)

Then, the dependency from time and position of each variable has been decoupled, allowing to express of the equation system in Eqs. (4, 5) as a linear differential system of \({\mathrm{X}}_{\mathrm{mn}}\), \({\mathrm{Y}}_{\mathrm{mn}}\) and \({\mathrm{W}}_{\mathrm{mn}}\), dependent only on the load applied \({\mathrm{Q}}_{\mathrm{mn}}\), and on the laminate properties.

$$\begin{aligned}&\left(\begin{array}{ccc}{\mathrm{D}}_{11}{\left(\frac{\mathrm{m}\pi }{\mathrm{a}}\right)}^{2}+{\mathrm{D}}_{66}{\left(\frac{\mathrm{n}\pi }{\mathrm{b}}\right)}^{2}+{\mathrm{k}}_{\mathrm{sh}}{\mathrm{A}}_{55}& \left({\mathrm{D}}_{12}+{\mathrm{D}}_{66}\right)\frac{\mathrm{m}\pi }{\mathrm{a}}\frac{\mathrm{n}\pi }{\mathrm{b}} & {\mathrm{k}}_{\mathrm{sh}}{\mathrm{A}}_{55}\frac{\mathrm{m}\pi }{\mathrm{a}}\\ & {\mathrm{D}}_{66}{\left(\frac{\mathrm{m}\pi }{\mathrm{a}}\right)}^{2}+{\mathrm{D}}_{22}{\left(\frac{\mathrm{n}\pi }{\mathrm{b}}\right)}^{2}+{\mathrm{k}}_{\mathrm{sh}}{\mathrm{A}}_{44} & {\mathrm{k}}_{\mathrm{sh}}{\mathrm{A}}_{44}\frac{\mathrm{n}\pi }{\mathrm{b}}\\ \mathrm{SIM}& & \left({\mathrm{k}}_{\mathrm{sh}}{\mathrm{A}}_{55}\right){\left(\frac{\mathrm{m}\pi }{\mathrm{a}}\right)}^{2}+\left({\mathrm{k}}_{\mathrm{sh}}{\mathrm{A}}_{44}\right){\left(\frac{\mathrm{n}\pi }{\mathrm{b}}\right)}^{2}\end{array}\right) \\&\left(\begin{array}{c}{\mathrm{X}}_{\mathrm{mn}}\\ {\mathrm{Y}}_{\mathrm{mn}}\\ {\mathrm{W}}_{\mathrm{mn}}\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ \mathrm{Q}-\rho \mathrm{h}\ddot{\mathrm{W}}\end{array}\right)\end{aligned}$$
(14)

One-to-one time-dependent relations of both \({\mathrm{X}}_{\mathrm{mn}}\) and \({\mathrm{Y}}_{\mathrm{mn}}\) with the vertical displacement field \({\mathrm{W}}_{\mathrm{mn}}\) can be obtained by linearly combining the first two equations of Eq. (14). Those relations allow the definition of the latter as the result of the Single Degree of Freedom (SDOF) differential equation [See Eq. (15)].

$${\ddot{\mathrm{W}}}_{\mathrm{mn}}\left(\mathrm{t}\right)+\frac{{\mathrm{L}}_{13}{\mathrm{K}}_{\mathrm{X}}+{\mathrm{L}}_{23}{\mathrm{K}}_{\mathrm{Y}}+{\mathrm{L}}_{33}}{\rho \mathrm{h}}{\mathrm{W}}_{\mathrm{mn}}\left(\mathrm{t}\right)={\mathrm{Q}}_{\mathrm{mn}}\left(\mathrm{t}\right)$$
(15)

The Eq. (15) is then analogue to the respective motion equation for an SDOF spring-mass system. Thus, it is possible to define the equivalent out-of-plane stiffness of the plate (i.e., the relation between the force applied in the centre and the deflection in that point), as in Eq. (16), where the definition of the equivalent Rayleigh–Ritz mass of the plate, \({\mathrm{M}}^{*}\) is based on the deformed shape function of the plate. Given that this research studies low-velocity impacts, the plate deformation will be dominated by the first vibration mode, and so, higher vibration modes may be neglected. Therefore, the equivalent Rayleigh–Ritz mass for a simply supported mass subjected to transverse is defined in Eq. (17). In our model, two webs and head have been projected to the skin to create a panel with variable stiffnesses. Then the panel is divided into sub-sections, and natural frequency has been calculated for each section. The summation of all these natural frequencies is equivalent to the natural frequency of the panel.

$$\mathrm{K}={\mathrm{M}}^{*}\cdot {\omega }_{\mathrm{mn}}^{2}={\mathrm{M}}^{*}\cdot \frac{{\mathrm{L}}_{13}{\mathrm{K}}_{\mathrm{X}}+{\mathrm{L}}_{23}{\mathrm{K}}_{\mathrm{Y}}+{\mathrm{L}}_{33}}{\rho\mathrm{h}}$$
(16)
$${\mathrm{M}}^{*}=\rho \cdot \mathrm{h}{\int }_{\mathrm{A}}{\left(\mathrm{sin}\left(\frac{\pi \mathrm{x}}{\mathrm{a}}\right)\mathrm{sin}\left(\frac{\pi \mathrm{y}}{\mathrm{b}}\right)\right)}^{2}\mathrm{dA}=\frac{\mathrm{M}}{4}$$
(17)

2.1 Stringer Stiffened Model

The geometry of super-stringer panels consists of a flat panel and an omega stringer, bonded together as illustrated in Fig. 1a. This work aims to develop a methodology to find the global stiffness of this structure using the aforementioned plate theory formulation, and so, be able to solve this impact problem using a double spring-mass system. For this purpose, and considering perfect bonding between plate and stringer, the assembly has been divided into five different parts: the single plate region without stringer (A), the stringer toes (B), the plate region in between the stringer toes (C), the stringer webs (D) and its wing (E). All of them extend along the panel's full length, so the dependency with the longitudinal coordinate y will be considered in the simple plate analysis. The damping effect is negligible during impact events, and it has a minor impact on penetration and acceleration histories. This is the main reason why we can simplify the mass-spring model and ignore the damping effect.

Fig. 1
figure 1

a Structure partitions and idealisation, b Stringer-stiffened panel equivalent spring model

Full size image

The proposed approach regards the different sections as simply supported plates whose supports are not fixed along but allow a vertical movement \(({\overline{\mathrm{w}}}_{\mathrm{B}}\) and \({\overline{\mathrm{w}}}_{\mathrm{C}}\) in Fig. 1a determined by the sections placed outwards of them. As a result, the analysis of a super-stringer panel is divided into five analysis for simply supported regular panels, which are to be called effective plates. Then, the stiffness figure for each effective plate (\({\mathrm{K}}_{\mathrm{i}}\)) are obtained as per the Eq. (16) and then combined to get the overall stiffness (\({\mathrm{K}}_{1}\)), considering this as the ratio between the applied force at the centre of region C and its displacement \({\mathrm{W}}_{\mathrm{mn}}\). The thickness and effective breadth used for each section are in Table 1, where the latter magnitude is defined as the nominal breadth times the respective \({\uplambda }_{\mathrm{i}}\) parameter. This parameter accounts for the rotation restriction at the interface regions, which causes the actual section stiffness to be higher than if simply supported conditions were actually happening. The determination of these effective length parameters is shown in Annexe A. Then, it is this effective section breadth that is used for the respective section’s natural frequency calculations, whereas, for mass calculations, the full breadth of each section as in Table 1 is used, as the total mass of the panel is to be accounted for in the problem dynamics.

Table 1 Effective plates thickness and width
Full size table

In this formulation, it is assumed that every region respective to the stringer (B, D and E) move equally for a given cross-section perpendicular to the stringer axis. So, they are considered as associated in parallel. Therefore, this group of springs is serially associated with those from sections A and C, given that they are assumed to be subjected to the same force. A schematic representation of the equivalent spring model is shown in Fig. 1b. However, as the stiffness obtained by Eq. (16) refers to the displacement at the centre of the plate, an equivalent section stiffness \({\mathrm{K}}_{{\mathrm{i}}_{\mathrm{eq}}}\) is defined so that the relation between the vertical displacement of a point different than the centre and the applied force is known [See Fig. 1b]. Then, the equivalent stiffness of the stringer-stiffened panel, \({\mathrm{K}}_{1}\) is obtained from the A to E’s section stiffness as in Eq. (18). The model is developed based on the 'Thin-walled' laminate assumption and can consider the summation of skin and stringer thicknesses. In this work, the thickness of the stringer and skin is assumed to be the same. The laminate consists of 10 plies with a thickness range of 0.3 to 0.6 mm, giving a plate thickness range between 3.0 to 6.0 mm.

$${\mathrm{K}}_{1}={\left[\left(\frac{\mathrm{sin}\left(\frac{\uppi {\mathrm{a}}_{\mathrm{A}}}{\mathrm{a}}\right)}{{\mathrm{K}}_{\mathrm{A}}}\right)+\left(\frac{\mathrm{sin}\left(\frac{\uppi {\mathrm{a}}_{\mathrm{B}}}{2{\mathrm{a}}_{\mathrm{B}}+{\mathrm{a}}_{\mathrm{C}}}\right) }{{\mathrm{K}}_{\mathrm{B}}+{\mathrm{K}}_{\mathrm{D}}+{\mathrm{K}}_{\mathrm{E}}}\right)+\left(\frac{1}{{\mathrm{K}}_{\mathrm{C}}}\right)\right]}^{-1}$$
(18)

To fully characterise the dynamic behaviour of the plate through a spring-mass model, the equivalent mass of the plate is determined through the Rayleigh–Ritz method. Equation (19) shows that the half-span normalised deflection shape is the deflection symmetric for the super-stringer midplane (\(\mathrm{x}=\mathrm{a}/2\)).

$$\frac{\mathrm{W}\left(\mathrm{x},\mathrm{y}\right)}{{\mathrm{W}}_{\mathrm{mn}}}=\left\{\begin{array}{cc}\frac{{\mathrm{K}}_{1}}{{\mathrm{K}}_{{\mathrm{A}}_{\mathrm{eq}}}}\mathrm{sin}\left(\frac{\pi \mathrm{x}}{\mathrm{a}}\right)\mathrm{sin}\left(\frac{\pi \mathrm{y}}{\mathrm{b}}\right)& 0<\mathrm{x}<{\mathrm{a}}_{\mathrm{A}}\\ \left(\stackrel{-}{\mathrm{w}}\left({\mathrm{a}}_{\mathrm{A}},\frac{\mathrm{b}}{2} \right)+\frac{{\mathrm{K}}_{1}\mathrm{sin}\left(\frac{\pi \left(\mathrm{x}-{\mathrm{a}}_{\mathrm{A}}\right)}{2{\mathrm{a}}_{\mathrm{B}}+{\mathrm{a}}_{\mathrm{C}}}\right)}{{\mathrm{K}}_{{\mathrm{B}}_{\mathrm{eq}}}+{\mathrm{K}}_{{\mathrm{D}}_{\mathrm{eq}}}+{\mathrm{K}}_{{\mathrm{E}}_{\mathrm{eq}}}}\right)\mathrm{sin}\left(\frac{\pi \mathrm{y}}{\mathrm{b}}\right)& {\mathrm{a}}_{\mathrm{A}}<\mathrm{x}<{\mathrm{a}}_{\mathrm{A}}+{\mathrm{a}}_{\mathrm{B}}\\ \left(\stackrel{-}{\mathrm{w}}\left({\mathrm{a}}_{\mathrm{A}}+{\mathrm{a}}_{\mathrm{B}},\frac{\mathrm{b}}{2}\right)+\frac{{\mathrm{K}}_{1}\mathrm{sin}\left(\frac{\pi \left(\mathrm{x}-{\mathrm{a}}_{\mathrm{A}}-{\mathrm{a}}_{\mathrm{B}}\right)}{{\mathrm{a}}_{\mathrm{C}}}\right)}{{\mathrm{K}}_{{\mathrm{C}}_{\mathrm{eq}}}}\right)\mathrm{sin}\left(\frac{\pi \mathrm{y}}{\mathrm{b}}\right)& {\mathrm{a}}_{\mathrm{A}}+{\mathrm{a}}_{\mathrm{B}}<\mathrm{x}<\frac{\mathrm{a}}{2}\end{array}\right.$$
(19)

Therefore, given that the thickness only depends on the coordinate x, as it is constant within sections but different in each one, the equivalent super stringer mass, \({\mathrm{M}}_{1}\) is defined as in Eq. (20). The definition of the laminate thickness across the plate is given in Eq. (21), wherein the interval respective to the stringer toe (\({\mathrm{a}}_{\mathrm{A}}<\mathrm{x}<{\mathrm{a}}_{\mathrm{A}}+{\mathrm{a}}_{\mathrm{B}}\)), an equivalent constant laminate thickness is defined so that, for a length equal to \({\mathrm{a}}_{\mathrm{B}}\), it weighs the same as the regions B, D and E. Once the equivalent effective stiffness and mass of the super stringer are known, the problem is solved in the same way as in the constant thickness laminate previously analysed.

$${\mathrm{M}}_{1}=2\uprho {\int }_{0}^{\mathrm{b}}{\int }_{0}^{\frac{\mathrm{a}}{2}}\mathrm{t}\left(\mathrm{x}\right){\left(\frac{\mathrm{W}\left(\mathrm{x},\mathrm{y}\right)}{{\mathrm{W}}_{\mathrm{mn}}}\right)}^{2}\mathrm{ dxdy}$$
(20)
$$\mathrm{t}(\mathrm{x})=\left\{\begin{array}{cc}{\mathrm{t}}_{\mathrm{pl}}& 0<\mathrm{x}<{\mathrm{a}}_{\mathrm{A}}\\ {\mathrm{t}}_{\mathrm{pl}}+\frac{\left({\mathrm{a}}_{\mathrm{B}}+\frac{{\mathrm{a}}_{\mathrm{D}}}{\mathrm{cos}\left(\uptheta \right)}+\frac{{\mathrm{a}}_{\mathrm{E}}}{2}\right)}{{\mathrm{a}}_{\mathrm{B}}}{\mathrm{t}}_{\mathrm{st}}& {\mathrm{a}}_{\mathrm{A}}<\mathrm{x}<{\mathrm{a}}_{\mathrm{A}}+{\mathrm{a}}_{\mathrm{B}}\\ {\mathrm{t}}_{\mathrm{pl}}& {\mathrm{a}}_{\mathrm{A}}+{\mathrm{a}}_{\mathrm{B}}<\mathrm{x}<\frac{\mathrm{a}}{2}\end{array}\right.$$
(21)

2.2 Impact Modelling

According to the spring-mass model, the interaction between the plate and the striker is governed by the effective contact stiffness, \({\mathrm{K}}_{2}^{*}\). However, the actual contact force is determined by the Hertzian contact law in Eq. (22).

$$\text{F}_{\text{c}}\left( \text{t}\cdot\right)= \text{K}_{2}{\cdot\left(\alpha \left(\mathrm{t}\right)\right)}^{\mathrm{p}}$$
(22)

In this study, the striker considered is spherical and made out of isotropic material. Thereby, the definition of the contact stiffness for the Hertzian improved law is in Eq. (23). Where \({\upnu }_{23}\) stands for the composite Poisson's ratio in the 23 planes, \({\mathrm{E}}_{2}\) for its transversal Young's Modulus and \(\mathrm{R}\) for the spherical impactor ratio [12, 13].

$${\mathrm{K}}_{2}=\frac{2}{3}\frac{{\mathrm{E}}_{2}}{\left(1-{\upnu }_{23}^{2}\right)}\sqrt{{\mathrm{R}}_{\mathrm{striker}}}$$
(23)

Nevertheless, to use the spring-mass model, the effective contact stiffness \({\mathrm{K}}_{2}^{*}\) must be obtained [Eq. (24)]. The effective contact stiffness \({\mathrm{K}}_{2}^{*}\) is defined to simultaneously minimise, for the actual and linearised problem, the difference between peak forces and between impulses given a contact duration \(\mathrm{T}\). The expression for \({\mathrm{K}}_{2}^{*}\) proposed is shown in Eq. (25), where p is the Hertzian contact exponent, found out to be 1.5 for composite targets.

$$\mathrm{F}\left(\mathrm{t}\right)={\mathrm{K}}_{2}^{*}\alpha \left(\mathrm{t}\right)={\mathrm{K}}_{2}^{*}\left({\mathrm{z}}_{2}\left(\mathrm{t}\right)-{\mathrm{z}}_{1}\left(\mathrm{t}\right)\right)$$
(24)
$${\mathrm{K}}_{2}^{*}=\sqrt{\pi }\Gamma \left(\frac{\mathrm{p}+1}{2}\right)\frac{2\Gamma \left(\frac{\mathrm{p}}{2}+1\right)+\sqrt{\pi }\Gamma \left(\frac{\mathrm{p}+1}{2}\right)}{4{\Gamma }^{2}\left(\frac{\mathrm{p}}{2}+1\right)+\pi {\Gamma }^{2}\left(\frac{\mathrm{p}+1}{2}\right)} {({\alpha }_{\mathrm{m}}^{\mathrm{BC}})}^{\mathrm{p}-1}{\mathrm{K}}_{2}$$
(25)

\(\Gamma \left(\mathrm{x}\right)\) stands for Euler's Gamma function and \({\alpha }_{\mathrm{m}}^{\mathrm{BC}}\) for the maximum contact deformation, which depends on the boundary conditions and follows Eqs. (26). The super indices f, ss and c stands for free, simply supported, and clamped edges, respectively.

$$\begin{array}{c}{\alpha }_{\mathrm{m}}^{\mathrm{f}}<{\alpha }_{\mathrm{m}}^{\mathrm{ss}}<{\alpha }_{\mathrm{m}}^{\mathrm{c}}\\ {\alpha }_{\mathrm{m}}^{\mathrm{f}}={\left(\frac{{\mathrm{M}}_{1}{\mathrm{M}}_{2}}{{\mathrm{M}}_{1}+{\mathrm{M}}_{2}}\right)}^{0.4}{\left(\frac{\left(\mathrm{p}+1\right){\mathrm{V}}^{2}}{2{\mathrm{K}}_{2}}\right)}^{\frac{1}{\mathrm{p}+1}}\\ {\alpha }_{\mathrm{m}}^{\mathrm{c}}={\left({\mathrm{M}}_{2}\right)}^{0.4}{\left(\frac{\left(\mathrm{p}+1\right){\mathrm{V}}^{2}}{2{\mathrm{K}}_{2}}\right)}^{\frac{1}{\mathrm{p}+1}}\end{array}$$
(26)

For this study, \({\alpha }_{\mathrm{m}}^{\mathrm{ss}}\) is defined as the mean value of the maximum contact deformation for clamped and free edges, respectively. Initial conditions of both striker and plate are given in Eq. (27), where V is the striker's speed.

$$\begin{array}{c}{\mathrm{z}}_{1}\left(0\right)={\mathrm{z}}_{2}\left(0\right)={\dot{\mathrm{z}}}_{1}\left(0\right)=0\\ {\dot{\mathrm{z}}}_{2}\left(0\right)=\mathrm{V}\end{array}$$
(27)

The resultant motion of both the striker and the plate is now obtained by solving the equivalent two-degree-of-freedom spring-mass model.

$$\begin{array}{c}{\mathrm{z}}_{1}\left(\mathrm{t}\right)={\mathrm{A}}_{1}\mathrm{sin}\left({\upomega }_{1}\mathrm{t}\right)+{\mathrm{A}}_{2}\mathrm{sin}\left({\upomega }_{2}\mathrm{t}\right)\\ {\mathrm{z}}_{2}\left(\mathrm{t}\right)={\mathrm{C}}_{1}{\mathrm{A}}_{1}\mathrm{sin}\left({\upomega }_{1}\mathrm{t}\right)+{\mathrm{C}}_{2}{\mathrm{A}}_{2}\mathrm{sin}({\mathrm{w}}_{2}\mathrm{t})\end{array}$$
(28)

where:

$$\begin{array}{c}{\upomega }_{\mathrm{1,2}}^{2}=\frac{1}{2}\left(\frac{{\mathrm{K}}_{1}+{\mathrm{K}}_{2}^{*}}{{\mathrm{M}}_{1}}+\frac{{\mathrm{K}}_{2}^{*}}{{\mathrm{M}}_{2}}\right)\mp \sqrt{\frac{1}{4}{\left(\frac{{\mathrm{K}}_{1}+{\mathrm{K}}_{2}^{*}}{{\mathrm{M}}_{1}}-\frac{{\mathrm{K}}_{2}}{{\mathrm{M}}_{2}}\right)}^{2}+\frac{{{\mathrm{K}}_{2}^{*}}^{2}}{{\mathrm{M}}_{1}{\mathrm{M}}_{2}}} \\ \begin{array}{cc}\begin{array}{cc}{C}_{1}=\frac{{\mathrm{K}}_{2}^{*}}{{\mathrm{K}}_{2}^{*}-{\upomega }_{1}^{2}{\mathrm{M}}_{2}},& {C}_{2}=\frac{{\mathrm{K}}_{2}^{*}}{{\mathrm{K}}_{2}^{*}-{\upomega }_{2}^{2}{\mathrm{M}}_{2}},\end{array}& \begin{array}{cc}{A}_{1}=\frac{\mathrm{V}}{{\upomega }_{1}({\mathrm{C}}_{1}-{\mathrm{C}}_{2})},& {A}_{2}=\end{array}\end{array}\frac{\mathrm{V}}{{\upomega }_{2}\left({\mathrm{C}}_{2}-{\mathrm{C}}_{1}\right)}\end{array}$$
(29)

Therefore, it is now possible to calculate the contact force history. For this research, the effective contact force \(\mathrm{F}\left(\mathrm{t}\right)\) is to be considered the actual contact force, as they are, by definition, similar, and the definition for the second one is much simpler. Then, the history of the contact force is given by Eq. (30).

$$\mathrm{F}\left(\mathrm{t}\right)={\mathrm{K}}_{2}^{*}\left({\mathrm{A}}_{1}\left({\mathrm{C}}_{1}-1\right)\mathrm{sin}\left({\upomega }_{1}\mathrm{t}\right)+{\mathrm{A}}_{2}\left({\mathrm{C}}_{2}-1\right)\mathrm{sin}\left({\upomega }_{2}\mathrm{t}\right)\right)$$
(30)

3 Validation of Theoretical Model

Before applying this model to a stringer-stiffened panel, it has been considered necessary to verify the procedure's accuracy in Sect. 3 compared to the literature results for a flat plate. In that regard, the research carried by Vaziri et al. [15] has been used for model validation where two different impactors initial velocity have been used. Composite material and impactor properties are shown in Table 2.

Table 2 Carbon Fibre Reinforced Polymer (CFRP) lamina properties and dimensions of the composite plate and the impactor
Full size table
Fig. 2
figure 2

Contact force history for different impact velocities: a 1.76 m/s, b 2.68 m/s

Full size image

As shown in Fig. 2a, b, good agreement is observed between the present study and the existent literature for the considered low energy impacts (\(\sim 9.5\mathrm{J}\) and \(\sim 22\mathrm{J}\)). It is important to highlight that the original author has filtered the results shown by Vaziri et al. [15] to reduce the high-frequency vibration. This is the reason for the difference in the behaviour within that spectrum. Thereby, the results obtained from the herein developed analysis for low-velocity impacts on composite plates represent reality. However, given that no failure mechanism has been considered, in case this latter event happens, the obtained solution would pose an upper bound for the actual contact force.

4 Finite Element Modelling

In order to prove the accuracy of this low fidelity method, the obtained results are compared with those from a finite element based dynamic and explicit model. As no experimental results are available for a low-velocity impact on a super stringer panel, the results from Vaziri et al. [14] are used to perform a mesh sensitivity analysis on a simply supported plate. The impact characteristics, plate and impactor dimensions and material properties are shown in Table 2. The striker is discretised with 8-noded brick elements (type: C3D8), whereas 8-noded continuum shell elements (type: SC8R) are used for the super stringer. To account for the interaction between the impactor and the plate, hard contact frictionless formulation is used. No viscous regularisation scheme has been used to avoid interference on the final results. Given that the kinetic energy is dominant in this analysis over the internal energy, no mass-scaling is used. Subsequently, the mesh sensitivity analysis outcome is shown in Fig. 3a, where it is proved that the model setup leads to accurate results to those available in the literature.

Fig. 3
figure 3

a Details of Finite Element Model (FEM) of super-stringer with an optimum element size of 2.5 mm, b Finite element mesh validation (V = 1.76 m/s), and c) Analytical super stringer method validation against finite element results (V = 3 m/s)

Full size image

Therefore, for a super stringer as described in Fig. 4, and using the model setup from the simply supported plate validation in Fig. 3a, the contact force history from the finite element model is shown in Fig. 3b. After the parameter \(\uplambda\) calibration, the present study's force history results show a good correlation with those for the finite element model, thus proving its accuracy. Parametric studies on the impact velocity, impactor mass, and laminate characteristics are then to be performed.

Fig. 4
figure 4

The geometry of the super-stringer panel. Dimensions in millimetres

Full size image

5 Results and Discussions

The resolution of this problem has been carried out in three sections: firstly, a study of the effect of the plate layup configuration, secondly, the effect of ply thickness and thirdly, the effect of the impactor mass and velocity. For the first analysis, six different symmetric layups shown in Table 3 have been applied to the panel. For this study, the stringer layup was set to that of laminate 4, as it is the most suitable for a stringer among the considered, given that it presents the maximum number of fibres along with it. Laminate 4 has more numbers of 90o orientation (fibre along the stringer), showing more resistance against buckling (in-plane load). The different configurations were then compared using the response in Eq. (31). The plate and impactor material properties are in Table 2, whereas the super-stringer dimensions are in Fig. 4. A spherical impactor of radius 0.01 m and weighing 6.14 kg is used, and the impact velocity is set to 3 m/s.

Table 3 Stacking sequences and layups
Full size table
$$\upxi =\frac{1}{1-{\left(\frac{\upomega }{{\upomega }_{\mathrm{mn}}}\right)}^{2}}$$
(31)

Therefore, for every possible plate layup, the gain response studied through a ply-thickness range from 0.3 to 0.6 mm, being their gain distribution shown in Fig. 5a. These results show two different optimum plate layups depending on the ply thickness considered, being the laminate 1 slightly better than laminates 5 and 6 for thin plies, whereas the latter showcases a superior performance for thicker laminates. As the improvement for thin plies provided by laminate one is minor, and the laminate 6 is better over a greater range, this one is considered the optimal plate laminate for subsequent analysis.

Fig. 5
figure 5

a Gain for several plate laminates and thicknesses. Stringer layup set to Laminate 4, b Super stringer contact force history for several ply thicknesses combination

Full size image

The contact force history has been studied for a super stringer with the laminate 6 in the plate.

A ply thickness range from 0.3 to 0.6 mm was applied to the stringer and the panel independently. Force–time curves for each combination have been plotted together in Fig. 5b, and gain results have been represented for each stringer and plate ply thickness in 3D graphs in Fig. 6a. Force history curves in Fig. 5b showcase a greater dependence on the plate ply thickness than that of the stringer. This results from the effect that each plate's stiffness has on the overall super stringer stiffness. As the stringer is severely stiffer in the out-of-plane response than the plate, the overall super-stringer stiffness is dominated by that of the plate [see Eq. 19]. Therefore, an increase in the stringer stiffness will be less effective than an increase in the plate stiffness. The stiffer the plate is, the lower the contact period is and the higher the peak force, given that the impulse remains.

Fig. 6
figure 6

Gain variation with a both stringer and plate ply thickness, b both impactor mass and initial velocity

Full size image

Gain curves showcase greater slopes for slenderer thicknesses in Fig. 6a, where a differentiated dependency with either the plate or the stringer ply thickness is appreciated. For a given stringer ply thickness, the gain shows an inverse garland-like distribution; for a fixed plate ply thickness, the distribution is much closer to a linear one, where the gain increases with the stringer ply thickness. Moreover, for thick plate thicknesses, the gain response presents a minimum value. For the configuration considered, this minimum value is obtained for plate and stringer ply thicknesses of 0.6 mm and 0.385 mm, respectively. Regarding the plate thickness variation, the slope of the curves increases as the plies become thicker; that is, for plate thicknesses from the first half of its range (from 0.3 mm to 0.45 mm, approximately), the gain has a significant decrease with thickness increase, compared to the second half.

Regarding these results, it is concluded that the thicker a laminate is, the more impact resistance it possesses. Nevertheless, given that for many applications, such as aerospace, weight is a crucial factor, further considerations have to be taken into account when deciding the optimal ply thickness. Therefore, the design parameters that quantify the quality of a configuration have to be defined by the designer according to design objectives and requirements.

Moreover, force–time curves variation with the mass and the impactor body's initial velocity are presented in Fig. 7a, b, respectively. It can be appreciated that increasing the striker mass also increases the peak force and contact period, given the larger impulse to be absorbed and the larger impactor inertia. On the other hand, the variation of the initial velocity of the striker showcases that the peak force grows with the velocity, whereas the contact time is slightly reduced.

Fig. 7
figure 7

a Super stringer contact force history for several impactor mass, b Super stringer contact force history for several impactors’ initial velocity

Full size image

To understand the effect these variables have on other problem parameters, the impactor mass and velocity have been plotted against the gain. Figure 6b gives an overview of this dependency, mainly dominated by the impactor velocity. Both the gain surface slope and value increased as the impact velocity is reduced. Moreover, the gain response is considerably less sensitive to the impactor mass. Considering these results, it can be concluded that the super stringer response shows a better performance for smaller and faster impactors. It is essential to notice that this analysis on the gain does not concern the material failure nor the magnitude of the force applied, but the relation between the applied force and the dynamic response of the plate.

6 Conclusion

In the present study, a low-velocity impact analysis of stringer-stiffened composite panels has been performed. An analytical approach based on FSDT for plate behaviour and spring-mass model for contact interaction is used. Contact force history and vibration gain response are calculated using a MATLAB® code, which is the design parameter. To prove the accuracy and validity of the developed method, a comparison with results from both the available in the literature and a self-developed finite element model is performed, showing a good correlation with the developed analytical tool, which is then used to analyse the impact performance of simply supported super-stringers. Firstly, the effect of the stringer stacking sequence for each plate laminate is investigated, showing the optimal combination in every case. A parametric study in which the influence of different ply thickness of both plate and stringer and plate stacking sequences (stringer layup is set to the optimal) is performed. For the best plate and stringer layup combination, another two parametric studies are conducted. At first, the gain response is represented against plate and stringer ply thickness, which vary separately. This study aims to identify the optimum ply thickness combination. Secondly, a parametric study is performed for the mass and the initial velocity of the impactor, using both the force history curves and the gain variation representation to understand their effect on the system behaviour. Results showcase that the layup, among the studied ones, which leads to a better performance regarding gain response is [0/0/90/90/90]s for the stringer setup and [45/-45/90/90/90]s for the plate. For these layups, the contact force history study exhibits a greater dependence on the plate ply thickness than that of the stringer. It is shown that the thicker the plate is, the lower the contact period and the higher the peak force are. Considering the gain response, its variation has a diverse tendency for the plate and the stringer. Curves decrease as the plate ply thickness increases, whereas they increase for slenderer stringer ply thicknesses. The minimum gain value corresponds to a plate and stringer ply thicknesses of 0.6 mm and 0.385 mm. Finally, from the mass and initial impactor velocity study, it is understood that both parameters make the peak force higher as they increase.

In contrast, they exhibit a differentiated behaviour for the contact period. While an increase in the impactor mass enlarges the contact duration, a larger impact velocity leads to a slight reduction. Concerning the gain response dependency on the impactor mass and velocity, it is proven that this magnitude is predominantly sensitive to the latter parameter. Finally, the developed models are capable to predict the low-velocity impact response of laminated composite panels with variable stiffnesses, layups and simply supported boundary conditions. However, the current model does not provide a complete representation of the complexity of real structures and is limited to an idealised geometry. There is still a gap in this field to introduce a robust mathematical model to represent a real composite structure, including details of sub-components and complex boundary conditions against impact loading.

7 Data Availability Statements

The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.