A nonlinear joint model for large-amplitude vibrations of initially curved panels: reduced-order modelling and experimental validation

This study conducts an extensive theoretical-experimental investigation into the nonlinear dynamical response of a base-excited initially curved panel clamped at two ends via bolted joints. A new distributed nonlinear joint stiﬀness model is proposed capable of exhibiting interconnected eﬀects at various states of the contact interface. More speciﬁcally, the proposed model allows the control of the nonlinear softening behaviour of the panel through controlling the displacement threshold for micro-slip and new stick states, as well as the rate at which the micro-slip region develops. Due to the inclined angle of the clamping supports, the panel is slightly curved in its fastened arrangement. Hence, the panel is modelled as a shallow shell using Donnell’s nonlinear shallow shell theory, taking into account von Kármán strain nonlinearities, while retaining all in-plane and out-of-plane displacements. Structural damping is considered via use of the Kelvin-Voigt model. Taking into account the work of the nonlinear joint stiﬀness model, the partial diﬀerential equations governing the in-plane and out-of-plane motions of the panel are derived using the generalised Hamilton’s principle, which are then discretised into a reduced-order model using a two-dimensional Galerkin modal decomposition approach. For the experimental part, two nominally identical panels are considered and several tests are conducted through base excitation in the primary resonance region and the backbone curves are obtained directly via the Phase Lock Loop method. Extensive comparisons are conducted between experimental results and theoretical predictions for both backbone curves and time histories of the panel midpoint velocity and displacement. It is shown that the theoretical predictions are in very good agreement with the experimental results, with the proposed joint stiﬀness model being capable of predicting the signiﬁcant variability in the experimental results


Introduction
Having an accurate mathematical model is an essential need in structural dynamics in order to predict the dynamical characteristics of mechanical structures, as the industrial trend has been to shorten the design process and reduce the test and prototyping costs.However, developing accurate models has proven to be challenging in structural dynamics [1].The modelling process of an assembly is even more challenging when it takes place in the realm of nonlinear dynamics, where the traditional linear models and methods fail.There are various types of nonlinearities that could arise in mechanical assemblies such as geometric, inertial, material, and contact nonlinearities [2,3].Only geometric and contact nonlinearities are relevant to this study, and hence the focus of the introduction will be on these types of nonlinearities.
There are also different approaches to modelling nonlinear systems.The first approach is to use nonparametric methods [4,5,6] in which a model is simultaneously estimated and characterised during the modelling process.In contrast, parametric methods [7,8,9] utilise pre-selected or pre-defined models whose unknown parameters are estimated during the process.Each of these approaches has its advantages and disadvantages [10].
Thin-walled structures [11,12,13] are a group of mechanical systems that are very likely to show nonlinear behaviour.These structures can easily undergo large deformations which triggers geometric nonlinearities.Additionally, contact nonlinearities are also important if an assembly of thin-walled structures is being examined.Hence, due to the presence of complex nonlinearities, analysing the nonlinear large-amplitude vibrations of thin-walled structures requires accurate mathematical models and advanced numerical methods, and can be computationally expensive.However, investigating the dynamics of such structures is crucial for the design and optimisation of many engineering applications, such as aerospace structures, pressure vessels, and shells used in chemical processing [14,15,16].
There are different approaches to modelling thin-walled structures including geometric nonlinearities.
One approach is to use beam models to study such structures [17].The most common method for modelling geometric nonlinearities in beams is through the use of nonlinear Euler-Bernoulli beam theory which accounts for the large deflection of the beam structure.In this theory, the assumptions of small deformation and strain are relaxed, and the beam's curvature is included in the strain-displacement relationship [18].Another approach is to consider the structure as a shallow shell and utilise shell theories for modelling and analysing such structures [19,20,21,22,23].The curvature of the structure could be intentionally made during the manufacturing process, or may be caused due to geometric imperfections or by an initial pre-load applied to the shell resulting in a pre-stressed structure.The most common example of the geometric nonlinearity of a curved shallow-shell is the snap-through instability [24,25,26], where a shallow shell suddenly snaps from one stable configuration to another due to a small increase in the excitation amplitude.Wrinkling [27], bifurcation [28], and limit-point instability [29] are other examples of geometric nonlinearity in shallow shells.
There are different models and theories used to describe the behaviour of shallow shells, each with its own strengths and limitations [30,31,32,33].Classical shell theory assumes that the shell is a thin-walled structure with negligible transverse shear deformation [34].It is based on the Kirchhoff-Love hypothesis and provides a simplified approach to analyse the static and dynamic behaviour of shells.There are also firstorder and higher-order shear deformation theories which include the effects of transverse shear deformation and allows for the analysis of moderately thick shells [35,36].However, for thin shell structures, the added benefits of using such theories are negligible, especially considering the added computational cost.Instead, Donnell's non-linear shallow shell theory [37] is perfectly suitable for analysing thin shell structures of shallow curvature which is the case in this study.This approach retains the von Kármán nonlinearities associated with transverse vibrations and accounts for the curvature of the shell profile in the strains [30].
Alijani et al. [38] investigated the nonlinear forced vibrations of FGM doubly curved shallow shells with a rectangular base using Donnell's nonlinear shallow-shell theory, assuming the shell is simply supported.
Yazdi [39] used the same theory to estimate the nonlinear vibrations of simply supported doubly curved cross-ply shells; the geometric nonlinearities were modelled utilising the von-Kármán theory.Amabili [40] conducted a detailed comparison of different shell theories for analysing large amplitude vibrations of circular cylindrical shells.
Structure-related nonlinearities aside, non-ideal boundary conditions play a significant role in the dynamics of thin-walled assemblies.To simplify the analysis, much of the literature assumes ideal boundary conditions.Some studies, on the other hand, have taken a more realistic approach and investigated the effect of non-ideal boundary conditions on the dynamics of thin plates and shallow shells.Sobota and Seffen [41] employed a Föppl-von Kármán analytical model to identify stable shapes of an isotropic axisymmetrical shallow shells of arbitrary polynomial shapes; additionally, they explored how the boundary conditions affect the critical geometry for achieving bistable inversion.Wang et al. [42] carried out an analytical and experimental investigation on the dynamics of thin-walled panels with non-ideal boundary conditions.They described the deformation of the panels based on Kirchhoff-Love's assumptions, considering the positioning variation and clamping force in the assembly process.The dynamics of a laminated plate with viscoelastic supports was investigated by Yang and Yang [43].They considered the geometrical nonlinearity of the plate and used a Kelvin-Voigt model to simulate the viscoelastic behaviour of the flexible boundary conditions.
Their results showed it is essential to take the effects of non-ideal boundary conditions into account particularly in the nonlinear transient analysis of a laminated plate.Despite all these efforts focusing on the effect of non-ideal boundary conditions, there still exists a gap in addressing nonlinear behaviour imposed by bolted joints in the boundaries of shallow shells, using a distributed joint model.
Joints are an important part of modern mechanical assemblies, enabling engineers to build complex structures by assembling simpler parts together through jointed interfaces.Bolted joints provide a great flexibility, particularly in maintenance of such assembled structures.However, they add more complexity to the modelling and analysis of mechanical structures due to the contact nonlinearity in the joint interfaces [44].
Although there are numerous methods that are widely used to measure accurate asperities of surfaces using high-tech equipment [45], the dynamics of jointed interfaces is yet to be fully understood.This is in part due to the fact that not all the joint interfaces have similar asperity distributions, not even structures that look identical at macro scale.On the other hand, joint interfaces are often inaccessible [46], and therefore, what happens in large contact interfaces cannot be measured directly.There have been attempts to capture some useful measurements internal to an interface such as measuring contact pressure distribution using pressure films [46,47], or using strain gauges [48] and smart piezoelectric bolts [49] for the direct measurement of the real-time bolt tension.These measurements can give significant insight about the physics of bolted joints.However, such measurements themselves may bring additional uncertainties to the actual behaviour of joints [46].Hence, modelling the jointed interfaces in assembled structures is based upon the measured response of the structure at the macro scale, rather than the micro scale dynamics of the joint interface.
Considering other factors such as variations in bolt pre-loads [50], contact pressure distribution [51,52], varying amplitude-dependent contact area [53], bolts loosening during dynamical working conditions [54], and environmental conditions such as humidity and temperature [55], one can see that there is a large uncertainty in predicting the dynamics of jointed structures [56,57].
In this study, the nonlinear dynamics of a thin panel which is doubly clamped via bolted joints is examined extensively.For such a configuration, there are two major sources of nonlinearities: first the geometric nonlinearities arising from large-amplitude vibrations of the panel resulting in mid-plane stretching, and second the nonlinearities associated with the joint behaviour.Hence, the panel is modelled via Donnell's nonlinear shallow shell theory taking into the geometric nonlinearities due to the combined effect of initial curvature and mid-plane stretching.Furthermore, a new joint stiffness model is proposed to account for stiffness variation of the boundary conditions during large-amplitude oscillations, as well as the precompression that can be applied to the panel during assembly.Several experiments are conducted on two nominally identical panels resulting in time histories of the panel's midpoint velocity and displacements, from which the backbone curves are constructed.Extensive comparisons are conducted between experimental and theoretical results to assess the accuracy and robustness of the proposed model.The paper is organised through the following sections: the experimental set-up and procedure is detailed in Section 2, followed by a detailed model development in Section 3. The proposed joint model is examined in more detail in Section 4, discussing the effect of different joint model parameters and the effect of pre-compression.A thorough comparison between experimental and theoretical results is conducted in Section 5 followed by a concluding summary of the key findings in Section 7.

Experimental setup and procedure
The benchmark structure examined in this study is a doubly clamped initially curved thin plate shown in Fig. 1 which was originally designed for a Tribomechadynamics Research Challenge (TRChallenge) proposed in 2021 by the International Committee on Jointed Structures1 .The purpose of this challenge was to evaluate the capability of progressive modelling tools to predict the dynamic behaviour of mechanical structures with high levels of nonlinearity in the absence of prior knowledge of experimental readings from the structure.The benchmark structure was designed in partnership with industrial researchers to ensure that the design considerations mirrored industrial best practices and that the obstacles encountered in modelling the benchmark system represented the challenges encountered by industrial researchers.
The structure is composed of a thin plate of 1.5 mm thickness referred to as the panel in this paper, with each end sandwiched between a pillar and a blade at an inclined angle.The distance between the two pillars is 300 mm and the width of the panel is 100 mm.A set of six M6 bolts and washers are used to clamp the panel ends between the blades and pillars.The support is clamped to the ground via a set of M10 bolts.The panel is arched in its fastened arrangement due to the flat contact surfaces with the pillars being intentionally mismatched by an angle of around 2 degrees about the y-axis on each side.This is to ensure that all panels have very similar curved profiles when assembled, irrespective of their initial natural curvature.The inclination of the panel can be seen in Fig. 1(b) where a horizontal yellow centre line is added to highlight the panel curvature.The support and blades are composed of hardened steel, while the panel is made of stainless steel.All the results output of the TRChallenge and the different prediction approaches used by the participants can be found in [58].The benchmark structure shown in Fig. 2(a) was manufactured and tested at the University of Stuttgart.
A unique assembly procedure was followed to tighten the M6 bolts at a torque level of 17 Nm in order to minimise the variability in the joint system.More specifically, the bolts were tightened in a specific order as shown in Fig. 3 for every assembly-disassembly procedure.Further details of the assembly procedure are described in [58].The support structure is clamped to a slip table of a large shaker to apply a base excitation to the system in the z direction, i.e. the transverse direction of the panel.is measured using a Multi-Point Vibrometer (MPV).Finally, a set of PCB-accelerometers are used to track the acceleration of support, blade and shaker table to ensure a good transmissibility of the energy through the system components.For data acquisition and control, dSPACE ControlDesk and a MicoLabBox is used as a rapid control prototyping system.The details of the experimental setup are shown in Fig. 4(a).
Phase-resonant testing, also known as a Phase Lock Loop (PLL), is used to identify the lowest frequency nonlinear mode.This approach is based on the Single Nonlinear Mode Theory, that provides the theoretical foundation for this method, and is a useful and reliable technique for the characterisation of nonlinear mechanical systems [59,60].This controller is a force control strategy whereby the distinct harmonics of To ensure repeatability of the measurements, a number of direct repeats in addition to repetitions that included disassembly and reassembly in between, were carried out.Furthermore, two nominally identical panels as shown in Fig. 2(b) are used in the testing, and due to the panel symmetry, it was possible to test the two sides of each panel.Due to the initial natural curvature of the panels, a convex and concave side are obtained for each panel.In this paper, the convex side will be referred to as side "a" and the concave side will be referred to as side "b".This results in four different configurations of the benchmark system summarised in the Table 1.Several preliminary tests were conducted to characterise the base excitation, the uniformity of the bolted joint contact interfaces using pressure films, and the accuracy of the output signals of different sensing systems used in the testing protocol [58].
In this study, only the results obtained for configurations 1 and 3 (based on Table 1) are analysed and compared to the theoretical predictions as explained in more detail later in Section 5.

Model development
As shown in Fig. 1(b), the panel in the test set-up has an initial arc-shape curvature due to the angle of the clamping pillars.Hence, the curved panel is modelled as a shallow shell, using Donnell's nonlinear doubly curved shallow shell theory.It should be noted that the panel has a shallow curvature with a radius of curvature of more than 8 m and a thickness of 1.5 mm, hence satisfying Donnell's shallow shell theory  requirement of h/R << 1.A general doubly curved shallow shell configuration is shown in Fig. 5(a) in which the shell is described in an orthogonal curvilinear coordinate system (x; y; z) with the origin O at one edge.The curvilinear coordinates are defined as x = R x θ x and y = R y θ y , with R x and R y denoting the principal radii of curvature and θ x and θ y representing the angular coordinates.The length of the shell edges are shown by a and b and the thickness is denoted by h.u(x, y, t), v(x, y, t), and w(x, y, t) show the shell's middle surface displacements in the x, y, and z directions.According to Donnell's nonlinear shallow shell theory, the nonlinear strain-displacement relationships can be formulated as In what follows, the nonlinear partial differential equations governing the motions of the panel are derived utilising the generalised Hamilton's principle, while considering a Kelvin-Voigt energy dissipation mechanism and a nonlinear joint model for the clamped ends.The stress-strain equations, taking into account the Kelvin-Voigt energy dissipation mechanism, can be written as yy + σ (vs) yy , in which E, η, and ν denote the Young's modulus, material damping coefficient, and Poisson's ratio, respectively.Superscripts (e) and (vs) are used to distinguish between the elastic and viscous components of the stresses, respectively.
The variation of the strain energy of the panel, based on Donnell's shallow shell theory, can be formulated as in which δ denotes the variational operator.
The kinetic energy of panel, taking into account the motion of the base (i.e.z 0 sin(Ωt)), can be written as where ρ denotes the mass density, and ∂ t denotes ∂/∂t.
The virtual dissipative work of the viscous components of the stresses can be expressed as σ (vs) xx δε xx + σ (vs) yy δε yy + τ (vs) xy δγ xy dz dy dx. ( The boundary conditions for the two clamped ends (i.e.x = 0 and a) are modelled such that they allow slipping in the x direction.More specifically, all rotational and translational degrees of freedom are fixed, with the exception of the displacement in the x direction.Furthermore, it is assumed that the structure experiences an axial compressive load, P x , at the two boundaries once the bolts are tightened.In fact, for such a structure, it is highly likely that tightening the bolts induces a pre-compressive or pre-tensile force.
Here, this residual load is modelled as a pre-compression.
The virtual work of the axial compressive load P x in the two movable clamped ends can be expressed as in which δ D represents the Dirac Delta function.
In order to model micro-slips in the boundaries, a distributed nonlinear joint model, J (u| x=0,a ), is proposed in this study.This model focuses only on the stiffness of the joints at the two boundaries via use of distributed nonlinear springs.The joint model can be expressed as in which k 1 u| x=0,a represents the linear stiffness of the joint (with k 1 being the linear stiffness coefficient and u| x=0,a being the in-plane displacement in the x direction at x = 0 and a) which is multiplied by a nonlinear function g given by This piecewise nonlinear function, along with its parameters (i.e. s, α, and c), is explained and examined in detail in the next section (i.e.Section 4) to maintain the flow of the model development in this section.
Here, we assume that in the assembled configuration, the structure is experiencing a constant precompression.Concurrently, the axial stiffness of the clamped boundary is characterised by a nonlinear distributed spring which is a function of the axial displacement, and is at its neutral state in the assembled configuration.Hence, the axial compressive load is applied prior to attaching the distributed nonlinear springs to the boundaries (i.e.assuming k 1 = 0) as shown in Fig. 6.To do this, a static analysis is required to be conducted to determine the new static equilibrium configuration of the structure as a result of the application of the compressive load which displaces the boundaries slightly in the x direction with the magnitudes u s (x = 0, y) and u s (x = a, y) in the two clamping ends.The formulation for the neutral configuration for the proposed joint function can then be expressed as J u(x = 0, y, t) − u s (x = 0, y) and J u(x = a, y, t) − u s (x = a, y) for the two boundaries.With this explanation, the virtual work of the proposed joint model at the two clamped ends can be formulated as Substituting Eqs. ( 3)-( 6) and Eq. ( 9) into the generalised Hamilton's principle and simplifying the instances of σ (e) xx + σ (vs) xx , σ (e) yy + σ (vs) yy , and τ (e) xy + τ (vs) xy to σ xx , σ yy , and τ xy , respectively, the equations of motion for u, v, and w can be derived as Expanding the stress terms and conducting the integration over the thickness yields the following nonlinearly coupled partial differential equations governing the in-plane and out-of-plane motions of the panel To be able to analyse the response of the system numerically, Eqs. ( 13)-( 15) need to be discretised into a set of nonlinear ordinary differential equations (ODEs) via use of the Galerkin method.More specifically, a modal decomposition method is utilised here where the mid-plane displacements are approximated through finite series expansions consisting of suitable spatial trial functions multiplied by unknown time-dependent generalised coordinates.
The types of boundary conditions are clamped at two edges (i.e.x = 0 and a) and free at the other two edges (i.e.y = 0 and b), noting that the clamped edges allow for motion in the x direction.Hence, the mid-plane displacements are defined via the following series expansions in which u m,n , ûm,n , v m,n , vm,n , w m,n , and ŵm,n are unknown generalised coordinates which are functions of time.The spatial trial functions Γ and Φ are the eigenfunctions for a linear free-free and clamped-clamped beam, respectively.For an arbitrary positive integer k, Γ k (y/b) and Φ k (x/a) can be expressed as where γ k = (cosh λ k − cos λ k )/(sinh λ k − sin λ k ) and λ k is kth root of the frequency equation cos λ cosh λ − 1 = 0.The selected trial functions satisfy the geometric boundary conditions of the panel.The terms ûm,n , vm,n , and ŵm,n are added to improve convergence and allow the solver to better satisfy the dynamic boundary conditions.Furthermore, given that the panel is free at two edges, i.e. y = 0 and b, there are rigid-body modes involved in the series expansions, namely the first term in the u (i.e. for n = 0) and v displacements, and the first two terms in the w displacement.It should be noted that w m,−1 is associated with the symmetric rigid-body transverse mode in the y direction while w m,0 is associated with the nonsymmetric one.Furthermore, Eq. ( 16) provides the general series expansions for the displacements for the boundary conditions of axially movable clamped in the x direction and free in the y direction.However, only the symmetric modes are considered in the discretisation process and later in the numerical simulations.
To ensure fully converged results, a total of 41 symmetric modes are retained in the series expansion, including 12 out-of-plane modes: w 1,−1 , w 1,1 , w 3,−1 , w 3,1 , w 5,−1 , w 5,1 , w 7,−1 , w 7,1 , ŵ1,1 , ŵ3,1 , ŵ5,1 , and ŵ7,1 , 13 in-plane modes in the x direction: u 1,0 , u 1,2 , u 3,0 , u 3,2 , u 5,0 , u 5,2 , u 7,0 , u 7,2 , u 9,0 , û1,2 , û3,2 , û5,2 , and û7,2 , and 16 in-plane modes in the y direction: v5,1 , v7,1 , v9,1 , and v11,1 .This results in a 41-degree-of-freedom discretised model of the panel and ensures converged results.Substituting the assumed displacements into the equations of motion, i.e.Eqs. ( 13)-( 15), and applying a two-dimensional Galerkin technique results in a set of 41 second-order nonlinearly coupled ODEs, which is first recast into a set of double-dimensional first-order ODEs and then solved through use of a pseudo-arclength continuation technique [61].It should be noted that during the discretisation procedure, the joint function cannot be integrated in a closed form, due to presence of generalised coordinates in the cos function argument and the absolute function as well as the piece-wise nature of the function.Instead, the joint function is kept as an unknown function and integrated numerically while retaining sufficient number of terms; the function is then calculated during the numerical simulation depending on the value of the in-plane It should be highlighted here once again that u s (x, y) in Eq. ( 13) is the static in-plane displacement in the x direction as a result of the application of the compressive load P x in the boundaries.Hence, in order to solve the equations of motion, i.e.Eqs. ( 13)-( 15), when both P x and k 1 have non-zero values, a static simulation should first be conducted for the desired value of P x while setting k 1 = 0.This allows obtaining the new equilibrium configuration of the structure through calculating the static displacements in the x, y, and z directions, i.e. u s (x, y), v s (x, y), and w s (x, y).Having obtained u s (x, y), the equations of motion are then solved for the desired values of P x and k 1 to obtain the backbone curves and time histories.
A note should be made here about the computations of the backbone curves.The proposed model is a 41-degree-of-freedom model with 12 out-of-plane (i.e.transverse) modes; however, the experimental results are reported for the transverse velocity of the panel midpoint, rather than each mode separately.Hence, in this study, in order to obtain the backbone curves of the panel's midpoint displacement and velocity, the corresponding frequency responses are initially obtained at various base acceleration levels.Then for each frequency response, the excitation frequency corresponding to a 90-degree phase difference between the force and displacement signals is extracted and the displacement and velocity amplitudes at that frequency are obtained, allowing for constructing the displacement and velocity backbone curves.

Joint modelling for large contact interfaces
In this section, we describe the stiffness function utilised to represent the frictionally clamped ends in the panel.It is assumed that the clamped ends are fixed in the out-of-plane z direction and the in-plane y direction while experiencing micro-slips in the axial (in-plane x) direction.This assumption can be justified as follows.In the z direction, the stiffness of the blocks used to clamp the panel is significantly higher than the lateral stiffness of the panel itself.Thus, the boundary condition can be considered rigid in the z  direction.In the x direction, the stretching resulting from the clamped ends generates a restoring force that induces local micro-slips.The stretching force component in the y direction is negligible compared to that in the x direction.Consequently, the boundary condition in the y direction is also considered fixed.
The typical graphs of the functions g(u) and J (u) are illustrated in Fig. 7 showing different regions in the contact interface.As seen, the boundary condition exhibits a dominantly linear behaviour for small displacements in the axial direction.This is shown as initial stick state in Fig. 7(b).However, as the displacement amplitude increases, small slip regions form within the contact interface.This leads In this paper, we propose a novel stiffness function designed for large contact interfaces, aiming to accurately capture the described behaviour.As depicted in Fig. 6(b), the boundary condition in the axial direction is modelled using a nonlinear spring.We consider the restoring force of these springs as the product of a linear stiffness k 1 u and a nonlinear piece-wise displacement function g(u) as in Eq. ( 8).The proposed joint stiffness function J (u) has the capability to capture both linear and nonlinear behaviours of a contact interface by incorporating a suitable g(u) function.The nonlinear g(u) function introduces adaptability to the proposed stiffness function, allowing for the generation of diverse behaviours within the contact interface.In this study, it is assumed that bolt shank pinning does not occur due to the fact that a clearance was designed between the panel and the bolt shanks in the benchmark tests to avoid pinning; however, the model can be modified via replacing the part for |u| > s with a suitable function to account for  pinning if the axial displacement at which pinning occurs is known.The g(u) function presented in Eq. ( 8) incorporates three control parameters: s, α, and c.These parameters exhibit interconnected effects at various states of the contact interface.s and α serve as the control parameters for determining the displacement threshold for micro-slip and new stick states, respectively as seen in Figs. 8 and 9, respectively.Furthermore, through changing both s and α simultaneously, the initial stick and the micro-slip states could be altered without affecting the final stick state (i.e. the latter hardening region in the backbone curve) as shown in Fig. 10.Parameter c controls the stiffness of the new stick state as shown in Fig. 11.In all these figures, P x is set to zero to highlight the effect of each joint function parameter.
Apart from the g(u) function parameters, the axial pre-compression caused by the bolted joints, i.e.P x , plays an important role in the nonlinear resonant behaviour of the panel as well.The effect of P x on the backbone curves of the panel is shown in Fig. 12.As seen, the softening behaviour of the panel becomes stronger as P x is increased.It should be noted increasing P x slightly increases ω 1,1 as well.In this study, the nonlinear fundamental natural frequency ω nl on the horizontal axis of backbone curves is normalised with respect to the linear natural frequency ω 1,1 for each case.Additionally, when P x is non-zero, the transverse displacement can be written as w(x, y, t) = w s (x, y) + w d (x, y, t) where w s denotes the static transverse displacement as a result of the application of P x and w d represents the transverse oscillation around the new deflected configuration.In all the backbone curves reported in this study, w d denotes panel's midpoint transverse displacement while ẇ represents panel's midpoint transverse velocity.
Finally, parameters P x and c can be used to control the transverse displacement asymmetry in the backbone curves.More specifically, as already shown, increasing the value of c reduces the initial softening behaviour while increasing P x strengthens the initial softening behaviour.Two cases are examined in Fig. 13, one with c = 1.5 and P x = 0 while the other with c = 4.5 and P x = 870 N/m.As seen, the transverse velocity backbone curve is almost identical for both cases but for the transverse displacement, the latter case shows a stronger asymmetry; this behaviour is shown in the time history plots for the velocity and displacement as well.Hence, P x and c allow for tuning the displacement backbone curves without impacting the velocity backbone response, which will be used in the next section to calibrate the joint function parameters as well as P x for different benchmark test cases.

Theoretical-experimental results and discussions
This section provides a detailed comparison between the theoretical and experimental results.A total of six benchmark test cases are selected for further analysis in this section, namely three benchmark tests for configuration 1 and three tests for configuration 3.For brevity, these test cases are referred to as A1-A3 and C1-C3, based on the naming in Table 2.Each experimental test is conducted at 46 different levels of base acceleration/frequency for both upsweep and downsweep cases.To construct the backbone curves for each case, first the transverse velocity time history is extracted for the last 0.25 s at each base acceleration/frequency level to ensure steady state results; the peak amplitude of the time history is then used to obtain one point on the backbone curve.This process is repeated for all base acceleration/frequency levels to obtain the backbone curve for one benchmark case, and then repeated for all six test cases.The time history for the transverse displacement is obtained via integrating the velocity time signal assuming that the displacement is zero when the velocity is maximum.More specifically, for each excitation frequency point on the experimental backbone curve, the last 0.25 seconds of the velocity signal is integrated to ensure a steady-state response and the zeroing out of the displacement at maximum velocity is done only once at the first peak of the velocity signal.The backbone curves for the transverse displacement are then constructed by extracting the maximum and minimum displacements amplitudes in the steady state time response.It is worth noting that for initially curved structures excited in the nonlinear regime, the steady state transverse displacement is not symmetric around the equilibrium configuration with the amplitude in the direction of the initial curvature always being smaller than that in the opposite direction of the curvature.
The values listed in Table 3 are used for numerical calculations.It should be noted that although the distance between the two clamping pillars is 300 mm, the length of the panel is considered to be 303 mm due to the fillet radius of the pillars at each clamping end which slightly increase the effective length of the panel.
The test results for all six cases are illustrated in Fig. 14, showing the backbone curves for the panel midpoint transverse velocity.It should be noted that the horizontal axis shows the normalised frequency ratio, with ω nl representing the nonlinear fundamental natural frequency and ω 1,1 denoting the linear counterpart.
Figure 14(a) shows the upsweep backbone curves for all cases while sub-figures (b)-(g) show both upsweep and downsweep results for each case separately.As seen, there is a significant variability between the results, with the results for cases A2, A3, and C1 being quite close and forming a band.Case A1 shows the weakest nonlinear softening behaviour among all cases examined here while case C3 shows the strongest softening behaviour.Here we have assumed that the parameters listed in Table 3 are the same for all cases, given that panels 1 and 2 are nominally identical, and that the examined configurations 1 and 3 are both for the same side of the panels (see Table 1).Hence, the variability in the results is considered to be due to the variability in the stiffness of the joints at the two clamping ends for each test as well as the pre-compression exerted by the joints on the panel in the assembly process.This variability is modelled via the proposed joint stiffness model J (u) = k 1 ug(u) explained in detail in Section 4 together with a pre-compression P x .It should be noted that the value of k 1 is assumed to be constant for all cases, and only the parameters of the nonlinear function g(u) and the pre-compression P x are varied to predict the experimentally observed behaviour.
Before discussing the comparison between experimental and theoretical results, it is worth explaining in more detail how P x and the parameters of the function g(u) are selected for each case.There are four parameters that need to be calibrated for each case, three g(u) parameters, i.e. s, α, and c, and the pre-compression P x .The calibration process is done in two stages; in stage 1, P x is set to zero and g(u) parameters are adjusted to obtain the best match between experimental and theoretical results, considering only the velocity backbone curve.Then, the values of c and P x are adjusted (as already explained in Fig. 13) to match the experimental and theoretical displacement backbone curves.
The initial calibration stage of g(u) parameters (assuming P x = 0) for test case C1 is shown in Fig. 15.Having examined the effect of different g(u) parameters in detail in Section 4, the values are adjusted here to obtain the best match between theoretical and experimental results.First, the value of k 1 is set to 3.0 × 10 9 N/m 2 , resulting in a linear natural frequency of 104.8 Hz which falls within the reported range of 103 Hz ±1% [58].Determining the value of k 1 is itself an iterative process which will be detailed after going through the joint function parameters calibration process.The parameters of function g(u) are calibrated for this case within four steps, with each row of sub-figures in Fig. 15 corresponding to one step.In step 1, i.e.Figs.15(a)-(c), g(u) is set to 1 to compare the numerical and experimental results.As seen, the experimentally observed softening behaviour is much stronger than the one predicted by the model assuming The way this is interpreted is that the stiffness of the joints in the clamping ends is reduced at large-amplitude oscillations due to micro slips, as already discussed in Section 4. Hence, in step 2, i.e. the   the displacement backbone curves will be considered and the values of P x and c are adjusted to obtain a match between those as well.For instance, it will be shown later that for case C1, the best match is obtained by changing c and P x to 2.4 and 1000 N/m.
The important point about this "manual" calibration process is that, through a deeper understanding of the physical effect of each g(u) parameter, the initial calibration of the parameters could be done within only a few steps.Then, the values of P x and c can be calibrated in the second stage through again only a few steps knowing how they impact the displacement backbone curves.It is worth highlighting this could be done by optimisation, but the physical insight would be lost and the number of steps might increase significantly.Now the process for obtaining the parameter k 1 is explained in more details.k 1 indicates the stiffness of the joint at low-amplitude (i.e.linear) oscillations.As such, its magnitude can be obtained through a linear analysis and comparison to experimental linear natural frequencies; however, we allowed for ±2% deviation in the natural frequency for calibration of k 1 due to different measurement uncertainties (such as uncertainties in the initial curvature profile and amplitude) and then utilised the nonlinear displacement backbone curve, and the magnitude of asymmetry between positive and negative branches, to obtain a final value for k 1 , i.e. 3.0 × 10 9 N/m 2 .This value is then kept fixed for all cases examined in this section.
The values of the parameters reported in Table 4 are near-optimal values that are obtained through manual tweaking of the parameters in several steps for each case, based on the physical meaning of each parameter as already explained above.If the parameters are determined via an extensive optimisation, the values might be slightly different, but we do not expect the match between the results to change much.The end goal here is not to ensure fully optimised joint parameters for a perfect match with experimental results, but rather to show that via tuning the physically meaningful parameters, the model could work well or very well in capturing the experimental observations.Hence, the focus has been on preserving the physical interpretation of the parameters to allow for an informed tweaking of these parameters, which could not only save time in calibrating the model, but also allows for improvement of the model.Furthermore, the main advantage of such a model is in its capability to predict the uncertain response of a jointed structure before having test results, and for that, the physical meaning of the parameters and their expected range become very important.Regarding an acceptable range for the parameters listed in Table 4, the key point that should be considered is that the parameters act in pairs, meaning that s and α could be changed in a way to affect only certain regions of the backbone curve, and same goes for c and P x , as already explained in detail in Section 4. Hence, if an optimisation procedure is to be used, the following ranges for the parameters could  (b)       mentioned before, the calibrated control parameters for each case are listed in Table 4. Again, a very good agreement between the experimental and theoretical results is observed showcasing the flexibility of the proposed joint function.The comparisons between the experimentally measured and theoretically predicted time histories for cases A3 and C1 are shown in Figs.20 and 21, respectively.As seen, there is an excellent match between the results for both cases.
The last two cases, i.e.C2 and C3, are examined in Fig. 22 for the backbone comparison and in Figs.23 and 24 for the time history comparison.As seen in the left column of Fig. 22, corresponding to case C2, there is again a very good agreement between experimental and theoretical velocity and displacement backbone curves.For case C3, i.e. the right column of Fig. 22, the model slightly overestimates the maximum amplitude of the midpoint transverse velocity and displacement at relatively large vibration amplitudes in the latter hardening region.Looking at the time histories, the match is quite good for case C2 (Fig. 23), while slight differences appear for both velocity and displacement signals for case C3 (Fig. 24).
It is worth highlighting that case C3 corresponds to the strongest softening behaviour in all the six cases examined in this study.Such strong softening behaviour, which is much stronger than the softening behaviour caused by the geometrically nonlinear behaviour of the initially curved panel, could indicate larger slipping in the joints compared to previous cases.Hence, it is possible that for this case, due to higher levels of slipping in the boundaries, pinning might be occurring at large oscillation amplitudes, where the discrepancy between theory and experiment is observed.However, as already mentioned in Section 4, pinning is not considered in the proposed joint model due to the fact that a clearance was designed between the panel and the bolt shanks in the benchmark tests to avoid pinning.Furthermore, adding a pinning capability to the joint model, although very much possible, will add more unknown parameters and makes the model as well as the calibration process more complex.The proposed model is considered a first step towards development of a new comprehensive joint model, and as such it has some limitations.However, as shown in this section, the model does an overall excellent job in capturing the experimental results with high accuracy, and allows for modelling the significant variability in the joint behaviour through physically meaningful parameters.
Furthermore, as seen in Table 4, the variability in parameters s and α is quite small, and mostly c and P x are varied to obtain a good match with experimental results, which shows the robustness of the proposed model.

Further discussion on joint interface parameters
This section further examines the effect of different joint interface parameters on the backbone curves, focusing mostly on c and P x parameters and their variability among different cases summarised in Table 4.
Reviewing the experimental results in Fig. 14 , it is seen that the velocity backbone curves for cases A2, A3, and C1 are relatively close to each other, yet the c and P x parameter values listed in Table 4 are quite different.
As already highlighted in Section 4, c and P x can be varied such that the velocity backbone curves remain mostly unchanged and only displacement backbone curves are varied, which is the main reason here for the variability of these two parameters for the mentioned cases.In other words, although the velocity backbone curves look similar, the change in the parameters is required to achieve the best match for the displacement backbone results as well.Considering that, one could expect a relatively good match if the parameters for case A2 are used for case C1 as well; but the match with the experimental results could be further improved by tweaking c and P x values.To further clarify this point, Fig. 25  case C1, but not so well for A3.However, even for case C1, the negative displacement backbone curve, i.e.
Fig. 25(c), shows visible differences between theory and experiment.Hence, it can be concluded that while parameters of case A2 can be used for case C1, the match between the results can be further improved by adjusting c and P x values.
To showcase the effect of parameter P x in the proposed joint model, the results for case C2, i.e. the left column in Fig. 22, are reconstructed here with and without P x parameter.More specifically, the calibration process is repeated with P x set to 0, and then the results of the model without P x are compared with those of the model with P x and the experimental results, as shown in Fig. 26.It is visible that the model which accounts for the pre-compression in the joint interface better captures the complex behaviour of the system and that the predictions of the model without P x are not as good as those of the model with P x , highlighting the significance of this parameter.
To further examine the performance of the proposed joint model, two additional cases for configuration 1 were examined (A4 and A5), but the results were not included in Section 5 for brevity.The experimental backbone curves for all cases are depicted in Fig. 27, showing that the velocity backbone curves of A4 and A5 are relatively close to those of A1 and A3, respectively.Furthermore, the calibrated parameters for A4 and A5 are listed in Table 5 whose values are in line with other configuration 1 cases.Comparing the calibrated parameters for A1 and A4 as well as those for A3 and A5 shows that only c and P x values are adjusted to obtain a good match with experimental results; more specifically, c is reduced to account for stronger softening behaviour and P x is reduced to account for weaker asymmetry in the displacement backbone curves for both A4 and A5 cases.10.0 0.9 0.9 0 A5

Concluding remarks
In this study the nonlinear large-amplitude vibration of an initially curved panel was examined theoretically and experimentally.For the theoretical part, a new joint stiffness model was proposed to account for the effect of micro-slips on the nonlinear resonant response.The joint model is capable of capturing the effects of transition from an initial stick state to micro-slip, and finally to a new stick state in the joint contact interface.In addition, the effect of axial pre-compression that can be caused by the bolted joints during assembly was accounted for as well.This model was used in the context of Donnell's nonlinear shallow shell theory to develop a nonlinear panel model capable of capturing geometric nonlinearities due to mid-plane stretching and initial curvature as well as joint stiffness nonlinearities.A 41-dof reduced-order model was constructed through use of the Galerkin modal decomposition method, which was then solved utilising a continuation method.For the experimental part, two nominally identical panels were base excited in the primary resonance region and the backbone curves for various cases were obtained through the phase lock loop method.Extensive comparisons were conducted between theoretical and experimental results, showcasing the capabilities of the proposed model.
Examining the effect of different joint parameters showed that the proposed model is capable of significantly altering the nonlinear behaviour in the primary resonance region without affecting the linear natural frequency.Furthermore, examining the effect of pre-compression showed that it is an effective parameter for controlling the asymmetry in the displacement backbone curves.It was shown that once the effect of the pre-compression and joint function parameters is well understood, a two-stage process can be followed to calibrate their values for each experimental test case.In the first stage, P x is set to zero and the values of g(u) parameters are calibrated within a few steps to reproduce the experimental velocity backbone curve.In the second stage, the values of pre-compression P x and joint parameter c are adjusted to match the experimental displacement backbone curves.Each control parameter has a physical meaning which affects the nonlinear resonance response in a predictable way, hence simplifying the calibration process.
Examining the experimental results showed that there is a significant variability between the responses of different cases, even for tests on the same panel.Given that the initial curvature magnitude was considered to be the same for all cases examined in this study, the source of the variability in the results was assumed to be the joint behaviour, including the joint stiffness and the axial pre-compression.In fact, uncertainty in the joint contact interfaces could play a significant role in the variability observed in the dynamic response.This uncertainty stems from various factors, including the deformation of asperities and the assembly process.The interaction between the contacting surfaces in a joint introduces variations in the distribution and magnitude of contact forces, leading to deviations in the overall behaviour of the system.It was shown that for most test cases, i.e.A1-A3 as well as C1 and C2, the predictions of proposed model were in excellent agreement with the experimental results.For the test case C3, i.e. the case with strongest softening behaviour, the model predictions were still in good agreement with the experiments with slight discrepancies between the results near peak amplitude oscillations.It was shown that the pre-compression parameter, i.e.P x , is required in order to achieve a good match with experimental results, especially for cases with strong softening behaviour such as C2 and C3.
Overall, it is shown that the proposed distributed joint stiffness model works very well in reproducing such scattered set of experimental results.Hence, this is considered a solid initial step towards development of a new comprehensive model for capturing uncertain behaviour of distributed joints in thin-walled structures.

Figure 1 :
Figure 1: (a) Design of the benchmark structure.(b) Top view.

Figure 3 :
Figure 3: Specific order of tightening the bolts.Four sensing systems are used in the experimental setup and described in the sketch in Fig. 4(b).A single Laser Doppler Vibrometer (LDV) pointed to the right-side blade is used to obtain the absolute response velocity of the support.A Differential LDV with two sensing heads (with one pointed to the middle of the panel and the second pointed to the left-side blade) is used to measure the relative response velocity of the panel to the support.The velocity response from 15 points distributed along the panel as shown in Fig. 2(a)
Fig. 4(b)) are brought into phase resonance using a Proportional Integral Derivative (PID) controller where the excitation frequency is adjusted to keep the phase lag equal to 90 degrees.A stepped voltage input is applied to the shaker to vary the base excitation amplitude.

Figure 5 :
Figure 5: (a) Schematic of a general doubly curved shallow shell in an orthogonal curvilinear coordinate system.(b) Schematic of the arch panel used for experiments, represented as a shallow shell.
in which ∂ x and ∂ y denote ∂/∂ x and ∂/∂ y, respectively.It should be noted that γ xy represents the engineering shear strain.The strain-displacement relationships in Eq. (1) are for a doubly curved shell.However, the panel examined in this study is curved only in the x direction as shown in the schematic in Fig.5(b).This means that R y = ∞ which results in w/R y = 0.In what follows, the equations of motion of the panel are derived assuming curvature only in the x direction and hence neglecting the term w/R y in ε yy .Furthermore, it should be noted that the panel is under base excitation of the form z 0 sin(Ωt) in the z direction, with z 0 and Ω being the amplitude and frequency of the base excitation, respectively.

Figure 6 :
Figure 6: The procedure for modelling the assembled panel, accounting for the pre-compression applied by the bolted joints and the joint stiffness: (a) applying axial compressive load in the boundaries in the absence of any constraints in the x direction, moving the structure to a new equilibrium configuration; (b) representing the joint stiffness via distributed nonlinear springs, ensuring the springs are in their neutral configuration.

Figure 12 :
Figure 12:The effect of axial pre-compression P x on the (a) transverse velocity and (b, c) maximum and minimum transverse displacement backbone curves.g function parameters are set to s = 14 µm, α = 0.9, and c = 1.5.

Figure 13 :
Figure 13:The combined effect of parameters c and P x on the (a) transverse velocity and (b, c) maximum and minimum transverse displacement backbone curves.(d,e) time histories corresponding to the peak amplitude on the velocity backbone curve.s = 14 µm and α = 0.9 for both cases.

Figure 17 :Figure 18 :
Figure 17: Time histories of the panel midpoint transverse (a) velocity and (b) displacement for case A1 at peak recorded oscillations (i.e.peak amplitude of the left column backbone curves in Fig. 16).[ ]: experimental results; [ ]: proposed model's predictions.(a)

Figure 16
Figure 16 shows the comparison between the experimental and theoretical backbone curves, with the left column showing results for case A1 and the right column showing those for case A2. [ ] and [ ] symbols illustrate the upsweep and downsweep experimental results, respectively, while the solid line [ ] shows the result obtained via the proposed model.Furthermore, the backbone curves are shown for the panel midpoint

Figure 20 :Figure 21 :
Figure 20: Time of the panel midpoint transverse (a) velocity and (b) displacement for case A3 at peak recorded oscillations (i.e.peak amplitude of the left column backbone curves in Fig. 19).[ ]: experimental results; [ ]: proposed model's predictions.(a)

Figure 23 :Figure 24 :
Figure 23: Time histories of the panel midpoint transverse (a) velocity and (b) displacement for case C2 at peak recorded oscillations (i.e.peak amplitude of the left column backbone curves in Fig. 22).[ ]: experimental results; [ ]: proposed model's predictions.(a)

Figure 25 :
Figure 25: Experimental backbone curves for cases A2, A3, and C1, and theoretical backbone curves based on the model calibrated for case A2; backbone curves are shown for the panel's midpoint (a) transverse velocity, (b) maximum transverse displacement, and (c) minimum transverse displacement.Symbols show upsweep experimental results, and solid line shows the theoretical predictions.

Table 1 :
Description of test configurations.

Table 2 :
Different benchmark tests examined in Section 5 and their short names.

Table 3 :
Geometric dimensions and material characteristics of the stainless steel panel.

Table 4 :
Calibrated joint function parameters and axial pre-compression for each benchmark test case.
). [ ]: experimental results; [ ]: proposed model's predictions.maximumtransversevelocityaswellas maximum and minimum transverse displacements.It is seen that the agreements between the experimental results and the theoretical predictions are very good for velocity and displacements.The experimental results for case A2 shows a stronger softening behaviour compared to case A1; this change in the softening behaviour can be robustly modelled via the proposed joint stiffness function and by adjusting P x and g(u) parameters.The comparisons between the time histories at peak recorded oscillations for A1 and A2 test cases are shown in Figs.17 and 18, respectively.It should be noted that the time histories are plotted for two periods of stead-state oscillations with t n denoting the normalised time with respect to the period of oscillation.As seen, the agreement between the experimental and theoretical results is very good for both velocity and displacement signals.Furthermore, it is seen that transverse displacement signal is nonsymmetric with respect to the equilibrium configuration which is a typical characteristics of initially curved structures.The theoretical and experimental backbone curves for cases A3 and C1 are shown in Fig.19.As is constructed showing the experimental backbone results for cases A2, A3, and C1 and the theoretical predictions based the model calibrated for case A2.As seen, the parameters for case A2 work generally well in capturing the experimental results for Experimental backbone curves (upsweep) for the panel midpoint velocity for all test cases.
0 Figure 26: Experimental and theoretical backbone curves for the case C2; backbone curves are shown for the panel's midpoint (a) transverse velocity, (b) maximum transverse displacement, and (c) minimum transverse displacement.[]/[]: upsweep/downsweep experimental results; solid line: theoretical predictions using the model with P x , calibrated for case C2 (as listed in Table4); dashed line: theoretical predictions using a model without P x , i.e.P x = 0, and c = 0.8 (with s and α remaining unchanged).

Table 5 :
Calibrated joint function parameters and axial pre-compression for cases A4 and A5.