Static model of a snowboard undergoing a carved turn: validation by full-scale test

The purpose of this study was to define a method for the validation of a numerical model representing a snowboard structure undergoing the conditions of a carved turn. A static load bench was developed to expose a snowboard to in-situ conditions. The deformed shape of the structure was measured with the use of retro-reflective markers, whose positions in space were tracked by six cameras and determined by triangulation. The experimental set-up was idealized in a finite element model, representing the composite structure and the loading environment. The model was validated by comparing the measured and computed displacement fields. The congruence between the two deformed surfaces was expressed by statistical means and constitutes the target function for optimization frameworks. Additionally, the contact pressure at the ground interface was experimentally assessed with the use of pressure measurement tape and compared with the numerical predictions.


Introduction
Structural analysis applied to thin-walled structures has been widely used in the past, especially in the aerospace industry. Nonlinear investigation methods applied to composite materials have been widely documented [1], and automation techniques have made them affordable yet applicable to other industrial sectors [2]. Regardless of the application, the reliability of the numerical representations must be assessed through a stringent validation procedure, to control the quality of the engineering solutions. The development of sports equipment, and more particularly composite snowboard structures, is taken here as an illustration of how the manufacturing environment can be numerically represented by creating a 'digital twin' of the product.
Previous research on the behavior of snowboard structures have proposed models to represent the stiffness at rest [3,4]. In addition, analytical simulations of the dynamics of a carved turn have been developed [5] and extensive knowledge is available to represent the on-snow conditions [6]. In these aforementioned studies, the validation of the numerical models consisted of comparing scalar values of the maximum displacements under bending and torsional loading independently. However, this approach is insufficient since it does not guarantee the validity of the entire deflected shape. Furthermore, it does not account for combined problems including contact formulations, as required in the present investigation. The conditions of a carved turn have been previously investigated to qualify the interface loads [7] and the behavior of a snowboard set on a carving angle [8]. Based on this work, the present paper proposes a methodology to validate a mathematical model by comparing the global displacement fields.

Experimental set-up
To simulate the loading conditions occurring on a snowboard during a carved turn, a static load bench was developed and measurements of deformations were performed on a simplified snowboard structure. The experimental set-up is described within this section.

Local coordinate system and conventions
The local coordinate system relative to the snowboard ( O xyz ) is defined as Cartesian (Fig. 1), such that the origin is positioned half-way between the centers of the binding locations ( B R and B F ), and coincident with the horizontal plane (H), averaged from the four contact point locations ( P i with i = 1, … , 4 ) on the bottom surface side of the snowboard. The lines P 1 P 2 and P 3 P 4 are referred to as 'contact limits'. The z-axis is normal to the plane (H), with the positive direction pointing towards the top surface of the snowboard. The global x-axis runs along the length of the snowboard, defined as the projection of B R B F onto (H). The y-axis runs across the width according to the right-hand rule. As a convention, the front tip of the snowboard (referred to as 'nose') points towards the positive x, whereas the rear tip (referred to as 'tail') points towards the negative x. The present conventions and reference systems are used throughout this paper unless mentioned otherwise.

Specimen investigated
The prototype consisted of a simplified snowboard structure manufactured according to the current industrial standards [9]. The structure featured a 5.55 mm constant thickness wood core made out of ash ( 60 mm-wide strip along the length of the snowboard, centered in the width) and poplar (remaining area), sandwiched between two composite skins made out of epoxy reinforced fiberglass (each 0.665 mm thick). The stacking sequence and corresponding fiber orientations are shown in Sect. 3.2. The snowboard geometry was defined as symmetric twin-tip (two planes of symmetry: xz-plane and yz-plane) and the sidecut had a constant radius of 7.05 m . The prototype's geometrical parameters are reported in Table 1 and illustrated in Fig. 1.

Static load bench
The test bench developed in-house is represented in Fig. 2. It consisted of two independent structural profiles (d) connected to the snowboard binding inserts via six screws each. Each profile was equipped with a support beam (i) that could be loaded with weight rings (j). The position of the weight rings was adjusted along the support beams to bring the system to the desired tilt angle, ranging from 40 • to 60 • . The angular position was measured via the reference profiles (e). The contact took place between the edge of the snowboard and a rigid, flat contact plate (g) made out of 1.5 mm thick aluminum.
Once the weights were in position, the angular position of the system was set with a threaded rod (c) and adjustment screws (a) until static equilibrium was achieved. Thus, the threaded rod was no longer transmitting neither axial nor bending load.
A schematic representation of the forces and moments applied to the snowboard is shown in Fig. 3. The application point of the input force P is represented at the center of the  It must be noted that the represented static equilibrium slightly differs from the static equilibrium of a snowboard undergoing an in-situ carved turn. According to [6], the input force induced by the skier or snowboarder is normal to the ski/snowboard surface (the force components in the xy-plane being small enough to be neglected). In the present study an artificial lateral force is generated at the binding interface, equal to the y-component of the input force P and proportional to the sine of the tilt angle : This loading condition including the extra force component is correspondingly regarded for the numerical representation of the experiment and throughout this paper. The impacts on the validation procedure are discussed in Sect. 6.

Measurement strategy
A motion tracking system (Vicon [10] and Vero v2.2 cameras) was used to capture the deformed shape of the structure. The measurement process consisted of positioning spherical retro-reflective markers on the top surface of the snowboard (such that they did not contribute to a stiffening of the structure) and measuring their position in space by triangulation, with the use of six cameras (see Fig. 4). The absolute positioning error of the measurement system is reported to be 0.15 mm for static measurements [11].
Due to the relatively small curvature of the deformed structure, the direct measurements of the marker positions were assumed to represent the deformation of the midsurface of the laminate ( max ≈ 5 × 10 −5 mm −1 leading to ∼ 0.06 % positioning error). The angular positions of the front and rear bindings with respect to the horizontal plane were measured during the experiment via the reference profiles (see Fig. 2), on which extremities two Vicon markers were positioned.

Evaluation scheme
The reference stance (distance between both bindings) was set to s = 600 mm for the present measurements. Symmetrical loading conditions were defined: the total loading weight was equally split on the front and rear bindings and The reaction force to the ground R is shown with its two components R y and R z Fig. 4 Front view of the load bench and the tested snowboard prototype equipped with 29 Vicon markers. The distance between both bindings denoted s (stance) is represented positioned such that the tilt angles on both bindings were similar. Twelve measurement sets were defined with different combinations of the total loading weight P (20, 40 and 60 kg ) and tilt angles (ranging from 39 • to 59 • ). The numbering of the measurement sets was established according to the performance order of the measurements (see Table 2).

Contact pressure
The contact pressure at the interface between the snowboard and the contact plate was experimentally determined for the measurement set exhibiting the highest predicted contact pressure (set A). A Prescale pressure measurement film [12] was inserted during the experiment between the edge of the snowboard and the contact plate (2-10 MPa measurement range). When pressure is applied, the measurement film turns red, the intensity of the color varying with the applied pressure. The measurement error of such films is reported to be approximately 10-15% for low contact pressure gradients [13]. It must be noted that the actual contact sizes are often overestimated by Prescale films [14].
The resulting colored film was scanned and the color densities converted into pressure units (MPa). A smoothing spline was fitted through the scattered measurements and the pressure profile along the longitudinal direction x could be assessed.

Finite element modeling
A numerical simulation of the experimental set-up was developed, representing the snowboard structure and the loading conditions occurring on the static load bench. The mathematical model and the verification procedure are described within this section.

Numerical representation
The snowboard prototype geometry and structure were idealized in a Finite Element Model (FEM) representing the experimental conditions (see Sect. 2.3). General-purpose shell elements (3 and 4 nodes) accounting for finite membrane strain were used to represent the snowboard. To describe the composite characteristics of the specimen, the laminated structure was represented by several layers consisting of an orthotropic material and accounting for each lamina's representative thickness. Hence, the stress and strain distributions within each layer could be accessed.
The tilt angles and the vertical loading P (see Fig. 3) were introduced via two reference nodes R R and R F (see Fig. 5), positioned relative to the snowboard at the center of mass of the system which includes the static bench and the weight rings. These reference nodes were coupled with the binding interface nodes (six nodes corresponding to the insert positions for each binding) via distributed coupling constraints, where the displacements and rotational degrees of freedom of the reference nodes were coupled to the average displacements of the binding interface nodes. In this way, the tilt angle and the input weight could be defined independently for the front and rear bindings, according to the experimental set-up (see Table 2).
To simulate the interaction between the edge of the snowboard and the contact plate, a rigid body plane was introduced with a linear pressure overclosure contact formulation. A friction coefficient of 0.4 was assumed between the composite snowboard and the rigid aluminum plate. Boundary conditions constraining the three translations were considered to eliminate all rigid body degrees of freedom from the model.

Material representation
In the snowboard prototype investigated, the thickness was constant over the entire structure (Sect. 2.2) and the stacking was defined as per Table 3. The material 0 • orientation corresponds to the longitudinal direction and a positive material orientation corresponds to a positive (counter-clockwise) rotation around the element normal (see Fig. 5).
The materials used were represented with linear elastic properties, evaluated in-house by laboratory testing and completed by additional literature references as follows: The modulus of elasticity of ash and poplar woods as well as E-fiberglass/epoxy samples was statistically determined by three-point bending tests. The values of Poisson's ratio Table 2 Evaluation scheme. Twelve samples were defined with different combinations of the input weight (P) and tilt angles, reported for the rear binding ( R ), front binding ( F ) and the average over both values ( ) and shear modulus were extrapolated from [15] for the wood samples, and from [16] for E-fiberglass/epoxy samples. The material properties used for the present investigation are reported in Table 4.

Analysis details
A static stress analysis was performed considering geometrical nonlinearities. The loading was applied incrementally within each analysis step, where the nonlinear equilibrium equations were solved using Newton-Raphson iterations. Model verification checks were performed according to [17]: a mesh convergence study was conducted, and advanced characteristics such as contact formulation and point load introduction were verified by representative benchmarks. The output results monitored were the nodal displacements and the contact pressure at the ground interface.
With an element size of approximately 10 × 10 mm , the FEM contained 4500 elements and over 30,000 degrees of freedom. The analysis required approximately 5 minutes utilizing the finite element software Abaqus® Academic Research release 6.13-2, on a 4-cores processor Intel® i7-6700 @3.40 GHz platform and a maximum RAM memory of 16 GB.

Validation strategy
To qualify and validate the FEM predictions, the numerical displacement fields were observed and compared with the experiment. First, the experimental outliers were detected and excluded from the measurements (Sect. 4.1). Then, the experimental displacement fields were superimposed with the numerical predictions (Sect. 4.2) and their equivalence    The polynomial coefficients of S were defined as the leastsquare fit of all markers positions, such that the sum of the squares of the vertical deviations was minimized:

Experimental outliers
For each marker position m i , the local normal direction to the interpolated surface S could be expressed as: The absolute distance along n i was computed for each marker and reported as w i , such that the distribution w = {w i ∈ ℝ | i = 1, … , k} represents the quality of the polynomial fitting. The root mean square error (RMSE) was computed as a scalar indicator to assess the fitting quality of the measurements: To identify possible outliers from the experiment (due to the partial obstruction of some markers by the surrounding structure), a recursive check was implemented where each marker was successively excluded from the measurement set, and the polynomial interpolation S performed anew. The changes in the residual RMSE were observed, and every marker showing a decrease in the RMSE higher than five percent were considered as outlier and excluded from the measurement set: (2) S(x, y) = A + Bx + Cy + Dxy + Ex 2 + Fy 2 + Gxy 2 + Hx 2 y + Ix 3 where RMSE (j) is the response value when excluding the marker j, and RMSE (all) the initial residual error considering all markers.

Superimposition of the displacement fields
As there exists no common reference coordinate system between the measurements and the numerical representation, both displacement fields were superimposed to enable their comparison. The method defining the optimal rigid body transformation to achieve the least-square superimposition is described within this section. The point cloud of the deformed FE nodes was regarded as the reference one. An isometric transformation was applied to bring it in the local snowboard coordinate system (see Sect. 2.1). In this way, the x-and y-coordinates described a unique node position. Only the portion of the surface included between the contact limits ( P 1 P 2 and P 3 P 4 ) was regarded as relevant for the following steps: Let A = {a i ∈ ℝ 3 | i = 1, … , n} represent the 3D point set of the FE results, where a i = a ix a iy a iz T represents the deformed nodal coordinates and n the total number of nodes. The polynomial surface S representing the experimental measurements was expressed as a point cloud resents the surface position at the FE nodal coordinates, such that b ix = a ix , b iy = a iy and b iz = S(a ix , a iy ).
To ensure correspondence between the numerical deformed results A ∈ ℝ 3×n and the interpolated experimental positions B ∈ ℝ 3×n , the elements of these two matrices were ranked in corresponding order.
The least-square rigid body transformation between the two point clouds could be expressed as follows [18]. With h = 1 1 ⋯ 1 T ∈ ℝ n , the corresponding matrices A and B are related with a translation t ∈ ℝ 3 and a 3D rotation R ∈ ℝ 3×3 by: If we compute the centroids (a x , a y , a z ) and (b x , b y , b z ) of the two sets A and B , respectively, and translate them in space so that the centroids coincide with the origin of the global coordinate system, A ′ and B ′ represent the translated point clouds and are defined as: where K = I − 1 n hh T is the n × n normalization matrix and I ∈ ℝ n×n the identity matrix.
The rigid body transformation can, therefore, be rewritten so that it depends only on R: The estimation of R results in the following minimization problem: where || ⋅ || Fro denotes the Frobenius norm of a matrix, defined as the square root of the sum of the absolute squares of its elements. One solution of this problem uses singular value decomposition [19], so that: in which U ∈ ℝ 3×3 is the matrix of left-singular vectors (corresponding to a rotation or reflection), ∈ ℝ 3×3 is a diagonal matrix containing the singular values (corresponding to a scaling) and V ∈ ℝ 3×3 is the matrix of right-singular vectors (corresponding to another rotation or reflection). The best-fit rotation matrix is then obtained from Eq. 12: The corresponding translation becomes: The homogeneous transformation is then applied to the point cloud of the measured marker positions, such that: , k} t h e t r a n sformed point cloud of the measured marker positions, h ′ = 1 1 ⋯ 1 T ∈ ℝ k and k the total number of markers.

Congruence between the experimental and numerical representations
Once the superimposition of the experimental and numerical displacement fields has been performed, their congruence can be assessed as follows: The absolute distance between a marker position and the deformed FE surface can be approximated by the distance along the normal direction of the closest finite element, respectively, for each marker of the measurement sets. Let n (1) i T represent the spatial coordinates of the three closest FE nodes to the marker m trans i . The two vectors v (12) i and v (23) i represent the local tangent plane to the deformed FE surface: The vector normal to the deformed FE surface is defined as v N i ∈ ℝ 3 , such that: The absolute positioning difference q i between the marker position m trans i and the numerical FE surface is calculated along v N i , as: Let q = {q i ∈ ℝ | i = 1, ⋯ , k} be the dataset of the absolute positioning differences q i , k being the total number of measured markers. The accuracy of the FEM predictions of displacements was evaluated through the statistical distribution of q : the mean absolute positioning difference q i , the corresponding 90%-quantile ( Q 90 ), the maximum value of the absolute positioning difference max(q i ) and the RMSE were computed.
To quantify the accuracy of the FE deformations predictions as a scalar relative positioning error, the mean absolute positioning difference q i was expressed for each measurement set as a percentage of the maximum longitudinal bending deflection u max , such that:

Results
This section presents the comparison between the experimental results and the predictions, to assess the validity of the numerical model. The quality of the polynomial fitting of the measurements is first assessed (Sect. 5.1), and then the displacement fields are compared (Sect. 5.2). The interface contact pressure is finally evaluated in a qualitative comparison (Sect. 5.3).

Interpolation of experimental measurements
For each measurement set, the polynomial surface S interpolated from the measured marker positions was represented (Fig. 6). Less than 15% of the markers were excluded from each measurement set by applying the outlier criterion defined in Eq. 6.
The interpolated polynomial surface S was used to superimpose the experimental measurements with the FE deformed node clouds. The direct comparison of the z-coordinate of the FE node cloud with S, the interpolated experimental surface ( z DIFF = z FEM − z S ), reveals the quality of the (18) fitting with regard to the FEM deformations (Fig. 7). The largest discrepancies were identified around the load application areas for all measurement sets ( B R and B F ), where local deflections under the bindings could not be captured in the experiment. The fitting showed a better accuracy along the contact edge (up to 0.3 mm absolute difference), whereas larger discrepancies were observed on the opposite edge (up to 1.4 mm between the binding locations).

Assessment of deformation predictions
The validity of the numerical predictions was assessed by direct computation of the normal distances between the experimental marker positions and the deformed FEM. The statistical distribution of q , the dataset of the absolute positioning differences and pos , the final relative positioning error are given in Table 5 for each measurement set.
The absolute positioning differences are presented by means of a box plot of all the measurement sets (Fig. 8), where the medians of the distributions (red strips) are graphically represented, together with the first and third quartiles (Q1 and Q3 corresponding to the box boundaries) and the sample maxima/minima (whiskers) excluding outliers (crosses). The outlier boundaries were set to (Q1 − 1.5 IQR) and (Q3 + 1.5 IQR) , with IQR denoting the interquartile range.
The mean absolute positioning difference ranged from 0.398 mm (measurement set A) to 0.756 mm (measurement set D), with an average over all measurement sets of

Assessment of the contact pressure
The measurement of the contact pressure along the edge line of the snowboard was performed for the measurement set A. The results exhibited four distinct pressure locations, denoted as CPR (contact point rear area), BR (rear binding area), CPF (contact point front area), and BF (front binding area). The raw measurements are shown on Fig. 9 together with the smoothing splines.
Over the rear half of the contact line, the mean value of the contact pressure was 2.31 MPa, about 10% higher than over the front half (2.09 MPa). The pressure distribution in the BF area showed a small discontinuity towards the front side of the area. The maximum pressure measured on the CPR area was 4.43 MPa, and the maximum value in the BR area was 2.56 MPa (peak ratio of 1.7 for the rear location). The maximum pressure measured on the CPF area was 3.42 MPa, and the maximum value in the BF area was 2.16 MPa (peak ratio of 1.6 for the front location).
The averaged width of the contact pressure line was approximated by direct measuring on the marked Prescale film. The line was 1.1 mm wide at the contact points, and 0.7 mm wide at the binding locations. The linear pressure output from the FEM was converted into a surface pressure accordingly, and the results are shown in Fig. 10. The peak ratios observed numerically were 1.9 for the rear location and 2.0 for the front location.

Discussion
The experimental measurements of the displacement fields were in good agreement with the numerical predictions (Sect. 5.2). Over the 12 experimental cases, the measured markers were on average within 0.559 mm to the FEM, with

Fig. 10
Numerical distribution of the contact pressure between the snowboard and the contact plate, for the measurement set A. The binding locations ( B R , B F ) and the contact points locations ( P 1 , P 3 ) are given below the chart as reference 90% of the measurements below 1.200 mm ( Table 5). The mean relative positioning error of the FEM predictions with respect to the maximal longitudinal bending deflections was established at pos = 2.57%. The RMSE of the final positioning error q i was shown to be a good scalar indicator to qualify the accuracy of the model. This indicator is suggested to be used as a target function to be minimized for the optimization of the model parameters (such as material and contact properties), for example in a sensitivity study framework.
The distribution of the experimental pressure along the edge of the snowboard showed a good qualitative correlation with the numerical predictions (Sect. 5.3): the location and spread of the pressure peaks were similar (Figs. 9, 10). The high-pressure peaks situated at the contact points were confirmed by the experiment and were 1.6 to 2.0 times higher than in the binding areas. The irregular distribution of the contact pressure suggests that there exist opportunities for further optimization of snowboard designs. For future works, more accurate experimental data on the interface contact pressure could be determined by the use of an edge load static bench, such as described in [7].
It was mentioned that the static equilibrium studied in this paper differs from the real conditions of a carved turn: an additional lateral force was induced at the binding locations (Fig. 3). It can be stated that since the model is validated including an extra force component, it may be applied to investigate similar load cases featuring simplified external load combinations, including an input force normal to the snowboard only. The validation can be generally extended to similar load cases as long as the assumptions underlying the shell element formulation remain valid (small transverse shear deformation, negligible through-thickness stresses).

Conclusion and future work
A method for the static validation of a numerical model representing a snowboard undergoing a carved turn was implemented, comparing the predicted global displacements with a full-scale laboratory test. Both deformed shapes were superimposed and the accuracy of the model could be quantified with a single scalar indicator, to be minimized for the calibration of the model parameters. Additionally, the contact pressure along the edge of the snowboard was qualitatively correlated with the experiment and showed a strong discontinuous distribution.
Various extensions to this work could be pursued. A dynamic approach to the modeling of a carved turn could be regarded, and the validation procedure extended to the assessment of modal parameters. Thus, the influence of geometrical and structural parameters on the dynamic performance characteristics of the snowboard such as chatter could be further investigated. The idealization of mechanical snow properties could also be implemented including an in-situ validation procedure, to simulate the real on-snow behavior. Additionally, the present numerical model could be used for inverse optimization frameworks, where the design could be steered by the desired structural response under operation, rather than the undeformed properties.