A fast water-slide layout approach using smooth surfaces, Saint Venant’s stationary water streamline, and passenger sliding interaction – theory, implementation, and validation

Water-slide manufacturers are constantly searching for innovative designs that maximize the rider’s fun while complying with safety requirements. However, the current state-of-the-art does not offer any simulation tool capable of aiding the design phase before manufacture in an affordable computational time. In this paper, we present a novel approach using smooth surface models together with Saint Venant’s filament theory to predict the water channel and, by this, a realistic estimation of passenger/water-slide interaction. The paper covers the theory and the experimental validation for the case of circular slide cross-sections. It is shown that the predicted water channel follows very well the measurement even without the influence of Saint Venant’s water-height terms and of the elliptical correction of the water cross-section for slanted water streamlines with respect to the centerline of the circular tube. Thus, for this kind of applications, a simple Bernoulli streamline model with proper water-slide friction parameters is shown to be sufficient. On the other hand, parameter identification of fluid-surface friction shows that three terms must be considered: a static one, a curvature-dependent correction, and a third correction factor for changes from low to large curvature (splashing effect). The results are also compared with fluid simulations using OpenFOAM. It is shown that the results of the simple model are not less accurate than the fluid-dynamics model, while needing only 1/100.000 of the computational time. The paper shows also the passenger/slide interaction by a multibody model including contact and passenger/fluid interactions. It is shown that the interaction between the passenger and water is significant for correct reproduction of real motion, showing that the water channel is relevant for computer-based water-slide layouts. As the water-channel parameters can be sufficiently modeled by Bernoulli’s streamline theory, the approach is independent of the slide cross-section and thus extendable to general cases. Also, as the computations run at 7/10 of real time, they can be used for online optimizations and passenger-specific actions on site.


Introduction
Water-slide manufacturers are constantly searching for innovative designs that maximize the rider's fun while complying with safety requirements. However, the current state-of-the-art does not offer any computational framework capable of aiding this search in an affordable computational time. There is a need for fast and reliable simulation tools allowing for an efficient computational approach to this fun/safety tradeoff, using, e.g., numerical optimization routines. To this end, three challenges need to be tackled: (1) A description of the geometry of the water-slide surface with an intuitive set of parameters of manageable size, (2) an efficient modeling of the dynamics of a sliding object, and (3) an efficient method for the approximate computation of the stationary water streamline on the slide and the interaction between water streamline and sliding objects.
In the currently available literature, the geometry of the water slide is either built as a concatenation of defined water-slide elements (e.g., straight or circular segments) with circular cross-section [1], or they are modeled using surface fitting with bivariate B-splines [2]. Whereas the first approach limits the variety of possible designs, the second approach leads to a large set of surface-related parameters, which create heavy computational load [3]. The dynamics of a sliding object is modeled either as a point-mass motion on the surface of the water slide using minimal coordinates [2,3] or as rigid-body motion under the influence of regularized contact forces in Cartesian coordinates [1]. The water streamline and its interaction with the sliding person is either not considered [2,3] or computed with a finite-volume method [4], which leads to very large computational times.
In this paper, we introduce a novel computational framework, which uses curve fitting with B-splines to model the differential geometry of the water-slide surface with circular cross-section, minimal coordinates, and surface joints to model the dynamics of the sliding person, stream-filament theory to model the water streamline, and friction and drag coefficients to model the interaction between water streamline and sliding objects.
The rest of this paper is organized as follows. Section 2 describes the proposed computational framework: In Sect. 2.1 the previously introduced concept of a surface joint is extended to model tube-like surfaces with constant cross-sections; in Sect. 2.2 a circular cross-section surface-joint is used to compute an approximation of the steady-state waterstreamline using Saint-Venant stream-filament theory; in Sect. 2.3 a simplified three-contactpoint kinematical model of a sliding person is introduced; in Sect. 2.4, friction and drag coefficients are used to model the interaction between the varying water streamline and sliding objects; and in Sect. 2.5, all these elements are put together to model the dynamics of a sliding person. Section 3 uses the computational tools introduced in Sect. 2 to model an actual water-slide ride and compares the simulation results with measurements and with a CFD simulation in OpenFOAM. Finally, Sect. 4 discusses the results.

Description of the method
In the rest of this paper, we use the following notation. A coordinate frame K i is composed of an origin point O i and three perpendicular directions x i , y i , and z i with z i = x i × y i . A matrix i R j represents the rotation matrix transforming the coordinates of frame K j into the coordinates of frame K i . A vector k i b j represents the physical property b of target frame K j , measured with respect to frame K i and expressed in the coordinates of frame K k . If the upper-left index is omitted, then the expression is valid for any coordinate decomposition.

Object/water slide interactions modeled as a kinetostatic transmission element
A multibody system can be regarded as a concatenation of kinetostatic transmission elements mapping motion from a set q in of input state objects to a set q out of output state objects at position, velocity and acceleration level [5]. The state objects comprise spatial reference frames and scalar variables. The motion transmission behavior described by the element comprises four suboperations position: q out = F ( q in ) , velocity:q out = J Fq in , acceleration:q out = J Fq in +J Fq in , force: where J F = ∂F/∂q in represents the Jacobian of the transmission element. This object-oriented abstraction was introduced by [5] to model the dynamics of general multibody systems. In the case of closed kinematical chains, the closure conditions are modeled by geometric measurements, denoted in this paper with the symbol g. Concatenated kinetostatic transmission elements are used to compute the equations of motion of the mechanical system in minimal form where vector q collects the minimal generalized coordinates,( ·) denotes the time derivative of (·), M(q) is the generalized mass matrix, b(q,q) denotes the generalized Coriolis, gyroscopic, and centrifugal forces, and Q(q) are the generalized applied forces. The surface joint, described in detail in [6] and depicted in its abstract behavior in Fig. 1, is a kinetostatic transmission element mapping the motion of an input frame K 1 and two surface parameters u and v to the motion of a Darboux frame K D moving on a surface defined in frame K 1 , at position, velocity, and acceleration level. The generalized coordinates u and v can be freely chosen and may be either translational (displacement) or rotational (angle) quantities. The natural choice for axis z D is the normal to the surface where r u and r v are the partial derivatives of the surface, and σ is a nonzero constant parameter, which can switch the orientation of the axis z D . The axes x D and y D can be freely selected as long as they lie in a tangent plane and are perpendicular to each other. The axis y D is chosen to be parallel to a linear combination of the partial derivatives r u and r v , and the axis x D completes the system to a right-handed coordinate frame, i.e., where α u and α v are constant parameters. This yields a well-defined Darboux frame K D even if one of the partial derivatives r u or r v vanishes. The choice of axes is different than that made in [6], but it suits better the intuition for this application, as explained later. The bivariate function r(u, v) and its partial derivatives with respect to the surface parameters u and v describe the differential geometry of the surface and the kinetostatics of the surface joint (motion and force transmission) as described in [6]. For general water-slide surfaces, the bivariate function r(u, v) can be computed using, for instance, surface fitting with bivariate B-splines as done in [6]. For tube-like water-slide surfaces with constant crosssection, which are the most common water-slide geometries, r(u, v) can be parameterized by defining a centerline 1 1 r T (u) = c(u) and a corresponding traveling tangent frame K T as a function of a first parameter u, as well as a cross-section T T r D (v) = s(v) as a function of a second parameter v, in coordinates of the tangent frame K T , as shown in Fig. 2.
There are some options for the parameterization of the tangent frame K T . In this work, we choose the Darboux frame parameterization described in [7], which yields a relative orientation with axis 1 x T being tangent to the centerline c(u), i.e., the axis 1 y T perpendicular to a horizon vector field h(u) as and the axis 1 z T completing the right-handed coordinate frame with To ensure that 1 R T (u) never becomes singular, the horizon vector function h(u) can be defined so that the horizon vector is never parallel to the tangent vector at u. The resulting parametric surface r(u, v) is given by Since s(v) = [0, s y (v), s z (v)] T , this can be rewritten as with the first-and second-order partial derivatives in which the dependencies on the surface parameters u and v have been left out for brevity. In this work, both the centerline c(u) and the horizon vector field h(u) are defined using the B-Spline curve-fitting routine CURFIT [8], and the cross-section s(v) is defined as a circle. The coordinate u is thus a "downstream" coordinate along the centerline of the slide, and the coordinate v is the cross-section angle defining the location of the origin of frame K D along the cross-section circle. The constants in Eq. (5) are chosen as α u = 0 and α v = 1, so that the axis x D is perfectly normal to the circular cross-section for any values of u and v. The constant σ in Eq. (4) is chosen such that the axis z D points toward the center of the circular cross-section. As shown in Fig. 2, this yields an intuitive and relatively small set of control-points defining the differential geometry of the water slide and hence the kinetostatics of a reference frame sliding on it. As shown later, dependencies of the water-channel cross-section are marginal (Saint-Venant terms) so that for future applications, general crosssections with a simple Bernoulli streamline theory can be used.

Computation of the stationary water-streamline for one inlet
In classical (nonstationary) Saint-Venant theory [9], one considers a one-dimensional shallow water flow along a tilted slope and takes into account a continuously varying water height h through the scalar equation Here w is the downstream velocity of a water particle, ζ is its coordinate along the water filament, and g is the gravitational acceleration. The term gS 0 represents the change of potential energy of the water particle due to the so-called "bottom slope" S 0 (= sin(α 0 ), where α 0 is the slope angle in "downward" direction). All other effects (friction) are subsumed in the term gS f , in which S f is termed the "friction slope" after rescaling the friction effects by premultiplication with g. This equation can be generalized to a general multibody equation for a particle of virtual mass m w = 1 kg moving along a curve on a surface (see Fig. 3) by the following steps.
1) As the streamline is stationary, we have 2) As the velocity w of the water particle is constant at any streamline coordinate ζ , it is a unique function of ζ , and thus from the chain rule it follows thaṫ which is the second term in Eq. (13). 3) A "natural" coordinate frame K w is introduced at the contact point of the point mass with the surface, with axis x w pointing downstream along the streamline (one-dimensional direction of Eq. (13)), axis z w pointing along the normal to the surface, and axis y w completing a right-handed coordinate frame. Moreover, a Darboux frame K M is introduced at the contact point of the cross-section according to Sect. 2.1, with axis x M pointing downstream along the tube centerline. 4) For a water particle column normal to the surface along the streamline, the only external forces acting on it are the weight m w g, the friction force −F R x w acting in opposite direction of the streamline velocity, and the normal force Nz w . 5) In streamline direction x w , the 3D equation must render Eq. (13) for the case of a straight slope. 6) In normal direction z w , the 3D equation must render the normal force N equal to the weight projection and the centripetal components due to the water velocity along the curved streamline. 7) In transversal direction an equilibrium between the centrifugal forces along the "cross" direction y w and the corresponding projection of the weight must hold. Putting all these steps together leads to the standard multibody equation with the following correspondence of the terms in Eq. (13) after premultiplication with m w : In addition, the term Nz w accounts for the normal constraint force needed for dynamic equilibrium. Note that for a straight slope, we obtain N/m w = g · z w , so that N/m w is in the classical Saint-Venant equation equal to the gravity component normal to the surface. Thus, for a curved streamline, N/m w is equal to the gravity term of Saint-Venant plus the effect of centrifugal forces, such as in the "g-force" effect known, for example, from roller-coaster design and motion simulation.
In the following, we further discuss the five terms.

1) Friction force
The friction force is velocity-squared dependent, as the Reynolds number for this application amounts to approximately Re ≈ 87500 for the water streamline (streamline velocity w ≈ 3.5 m/s, kinematic viscosity ν ≈ 1.0 × 10 −6 kg/(ms), hydraulic radius R ≈ 0.025 m). It is assumed here to consist of three parts: a constant term, a curvaturedependent term, and a term depending on the partial derivative of the curvature with respect to the downstream coordinate u. As depicted in Fig. 4, these terms account for the effects of constant friction, centrifugal forces, which can be found as pressure loss in bends in hydrodynamics textbooks [10], and "splash-effects" caused by water hitting a segment of high curvature (small bending radius) coming from a segment of smaller curvature (large bending radius), respectively. Hence where the curvature χ of the centerline and its partial derivative χ u with respect to the downstream coordinate u are given by with The Heaviside function H(χ u ) yields 1 for positive argument or zero else. The friction parameters f 1 , f 2 , and f 3 can be identified via numerical optimization (see Sect. 3). This work considers only one single inlet and a circular cross-section, but the described approach may be extended to several inlets and general cross-sections.

2) Acceleration term
The acceleration of the particle can be computed in terms of the streamline direction aṡ whereψ andθ are the components of the angular velocity of the streamline direction x w = w w along y w and z w , respectively. Equivalently, this can be written in terms of the surface coordinates aṡ w = r uü + r vv +â withâ = r uuu 2 + 2r uvuv + r vvv 2 .
Projecting Eqs. (16) and (21) in directions of y w and x w , respectively, yields which leads to the ordinary differential equations for the surface coordinates u and v, where the longitudinal streamline accelerationẇ still has to be determined.

3) Normal force
Multiplying Eq. (16) by z w leads to the normal force

4) Saint-Venant correction factor -continuity equation
The most involved term accounts for pressure change due to the change of water height. Saint-Venant's equation assumes that the channel bottom slope is small and the streamline only has small curvatures, which allows for the assumption of a hydrostatic pressure distribution [9]. Multiplied by the virtual mass m w , it yields a correction factor κ proportional to the corresponding virtual weight. For the case of a curved streamline, as in a water-slide, it is more realistic to define κ to be proportional to the "felt weight", which is equal to the normal force N . The Saint-Venant correction factor becomes The factor ∂h ∂ζ is calculated using the continuity equation Q Volume = Aw = const. Due to the angle β between the streamline velocity w and the tangent to the tube centerline (axis x M ), the water-channel cross-section is equal to an elliptical cross-section with semiminor axis r (radius of the tube) and semimajor axis a, as shown in Fig. 5. Since the height of the water channel h in the water slide is small, the ellipse is approximated by a circle of radius R, which corresponds to the radius of curvature of the covertex of the ellipse.
For the semimajor axis, we obtain and for the radius of curvature of the covertex, These can only become singular if the waterflow is perpendicular to the tube centerline tangent, which can be excluded for a water slide. Using the continuity equation ∂A ∂h b ∂h ∂ζ the partial derivative of the water height h with respect to the water position ζ can be computed as

5) Differential equation in streamline direction
Projecting Eq. (16) in direction of x w , inserting Eqs. (26) and (31), and collecting terms forẇ on both sides yield for the particle acceleration in direction of the streamlinė Inserting Eq. (32) into Eq. (24) yields a set of two ordinary differential equations of second order, which, together with the ordinary differential equation for the rate of change of the water heightḣ can be solved numerically for general initial conditions q 0 = [u 0 , v 0 ] T andq 0 = [u 0 ,v 0 ] T as well as an initial height h 0 at the water inlet. The solution yields the trajectory of frame K w describing the water streamline.
To compute the boundaries of the water channel, which are needed to compute the interaction of the water with the sliding person, two further coordinate systems K wL and K wR are first inserted, whose origins are shifted with respect to frame K w by the water height h in the z-direction and half the water width b/2 in positive and negative y-directions, respectively (see Fig. 6). Then two additional Darboux frames K D,wL and K D,wR with the respective positions [u wL , v wL ] T and [u wR , v wR ] T on the water slide are introduced, and the constraint conditions are solved. This positions the Darboux frames at a minimal distance to the two frames K wL and K wR (Fig. 6), yielding the sought water boundaries.

Modeling of the sliding person
The sliding person, assumed to slide with crossed legs and closed arms as shown in Fig. 7, is modeled as a rigid-body open kinematical chain consisting of a pelvis (frame K P ), a lower body connected to the pelvis by means of a revolute joint allowing for the flexion ϕ L of the legs, and an upper body connected to the pelvis by means of a revolute joint allowing for the forward tilt ϕ U of the trunk (see Fig. 7). The open kinematical chain is constrained to slide on the water-slide surface with three contact points at the pelvis, the heels (frame K L ), and the back (frame K U ), respectively. To impose these constraints, the pose of the pelvis frame K P is prescribed by its position [u P , v P ] T on the water-slide surface and a rotation ϕ P about the corresponding water-slide surface-normal z P . Two additional Darboux frames K D L and K D U with the respective positions [u L , v L ] T and [u U , v U ] T on the water-slide surface are introduced to formulate the closure conditions which yield a closed kinematical chain with two loops and three degrees of freedom as shown in Fig. 7. The latter are chosen as q = [u P , v P , ϕ P ] T of the pelvis frame K P on the water-slide surface.
In this work, we consider a person sliding with crossed legs and closed arms, but the presented approach can be extended to other sliding styles and to tire slides.

Friction and drag forces acting on sliding objects
To simulate the friction effects acting on a sliding object, the friction force is introduced, where μ D is a friction coefficient, N D is the normal reaction force acting on the corresponding Darboux frame K D , and w D is its velocity.
To approximately compute the forces resulting from the interaction between a sliding object and water channel, a sphere is fitted locally to the geometry of the object and then rigidly attached to the Darboux frame K D , as shown in Fig. 8, so that the drag force can be computed as Here C D is the drag coefficient, computed according to [11], ρ w is the density of water, (w D − w) is the relative velocity of the origin of the Darboux frame K D with respect to the water, and A c is the projection of the spherical cap under water onto a plane perpendicular to the relative velocity vector. The area A c is computed by finding the spherical cap under water and then projecting it according to [12]. The force is "turned on" when the Darboux frame K D enters the water channel and "turned off" when it leaves it. The exact points in time when this happens can be computed during integration if an integrator with root finding functionality is used. The spherical cap under water can be obtained in closed-form by computing the points S 1 and S 2 defining the intersection between water channel and sphere in the yz-plane of frame K M (see Fig. 8). Hereby the angle β between the tangent to the water streamline and the water-slide centerline is assumed to be very small such that frames K w and K M are almost parallel to each other, and the computed streamline velocity w and water height h may be used directly in the computations. Expressing the position of each intersection point S i in coordinates of frame K M as the vector sum and projecting it on z M with M r S i · z M = h yield the scalar equation Here γ represents the angle between the two tangents to the cross-section s(v) at the positions of K M and K D (see Fig. 8). Equation (39) delivers the closed solution for the two angles ϕ 1 and ϕ 2 defining the intersection points S 1 and S 2 .

Dynamics of a water-slide rider
The computational tools developed in Sect.  (2) with generalized coordinates q = [u P , v P , ϕ P ] T , where u P and v P describe the position of the contact point of the rider's pelvis on the water-slide surface, and ϕ P describes the rotational degree of freedom of the corresponding Darboux frame about the surface normal, as shown in Fig. 7 and Fig. 9. The corresponding vector of generalized forces Q(q) = Q G (q) + Q F (q) + Q D (q) contains the projection Q G (q) of the gravitational forces, the projection Q F (q) of the friction forces F F , and the projection Q D (q) of the drag forces F D (see Eq. (38)) acting on each of the sliding Darboux frames K D P , K D L , and K D U . The corresponding projecting Jacobians can be computed from the kinetostatics of the closed kinematical chain and are not shown here for brevity. In this work, the basis of the mass distribution, radii of gyration, and relative length of the body segments are taken from [13] and [14]. The radius of the spheres used to compute the drag forces are 0.5 m for the pelvis and upper-body and 0.05 m for the lower-body.

Parameter identification for the water-streamline model
The water streamline was measured using video cameras Canon EOS 2000D with a resolution of 640 × 480 (25 fps) and a Parrot ANAFI drone with a resolution of 4096 × 2160 (24 fps). To this end, marks were taped by hand on the water slide surface at 28 different cross-sections along the water-slide centerline and at every = 0.1 m along the crosssection profile. To compare the simulated results to the measurements, the middle of the water channel was determined at every measured cross-section i (see Fig. 10) and then parameterized by the angle v M i = M i /r. The velocity of the water streamline was found by observing table tennis balls carried along by the water (Fig. 11). At five straight sections, the velocity of the water w M i in the middle of the section was estimated by dividing the length of the section by the time it takes the table tennis ball to pass through it. Figure 12 shows the cross-sections at which position was measured for water-channel (light grey) and sliding-person (dark grey) parameter identification, as well as the water-slide sections at which the travel velocity was estimated (red hatched).

Method
A parameter identification was carried out for the stationary water-streamline model presented in Sect. 2.2 using numerical optimization. The corresponding nonlinear optimization problem was formulated as minimize with the friction parameters f 1 , f 2 , and f 3 as optimization parameters. In the cost function, v C i is the angle parameterizing the location of the middle of the computed water channel at the measured cross-section i. Similarly, w C i is the velocity of the computed water channel at the middle of the measured section, and w 1 is the corresponding weighting factor. The constant factor k 1 = 1 rad s/m is a rough approximation for the ratio of the expected range of the computed angles (here approximately 1 rad) to the expected range of the computed velocities (here approximately 1 m/s) and hence renders velocity costs in rad 2 . Overall, the optimization used n 1 = 27 sampling points for the cross-section angles and m 1 = 5 sampling points for the water velocity, and was carried out subject to lower and upper bounds l 1 and u 1 . The optimization problem was solved using the nonlinear optimization routine nag_opt_nlin_lsq (e04unc) [17], a sequential-quadratic programming least-squares solver, with the tolerances listed in Table 1. The integration of the equations of motion was performed with the routine LSODAR, which is part of ODEPACK [18], a Livermore solver for ordinary differential equations with automatic method switching for stiff and nonstiff problems with tolerances relTol = 1.0e-10, absTol = 1.0e-9, and output time-step 0.005 s.
To find initial guesses of the three friction parameters, the following steps were carried out.
1a) Estimate the friction parameter f 1 from a point mass sliding down the entire water slide with the velocity measured at the last straight water-slide section, without velocity changes, under the solely effects of gravity and squared-velocity proportional friction. For the given water-slide segment, this yields the energy equation where h is the height difference, u is the length of the water slide, and w M 5 is the water velocity measured in section 5 (see Fig. 12). This resulted in f 1 = 0.0689. 1b) Optimize the friction parameter f 1 with a simplified version of Eq. (43) considering only the water velocity (w 1 = 1) at the last straight water-slide section (for the given water-slide segment, section 5) with both f 2 and f 3 set to zero, using the value of f 1 estimated in (1a) as an initial guess. This resulted in f 1 = 0.0695467.    estimated in (1b) and (2), respectively, and using f 3 = 0 as an initial guess. This resulted in f 3 = 0.690685.
These initial guesses were then used to determine by optimization the parameter identification with different weighting coefficients w 1 , as shown in Table 2. Finally, a weighting coefficient w 1 = 0.8 was chosen, as here the difference between position cost and velocity cost is the smallest, with the velocity cost smaller than its position counterpart. The results are shown in Table 3 and Fig. 13. Figure 13 shows the position of the middle of the water channel and the velocity of the water-channel cross-section. Once the parameters are identified, a simulation run of the 22.3-s-long point mass motion requires a CPU-time of 1.0 s.

Discussion
To get an impression of the value of the friction coefficient f (see Eq. (17)), a sectional optimization of the water slide was performed. To this end, the water slide was divided in 50 equally sized sections. At each transition and at the beginning and end of the water slide, the friction coefficient f was optimized, which led to 51 optimization parameters. At the Fig. 13 Results of the position and velocity of the water channel compared to the measurement. The position is computed by subtracting π/2 from the transversal coordinate v C and multiplying the result by the waterslide radius r. A position of 0 m indicates that the water channel is on lowest position in the cross-section (transversal coordinate v C is π/2). In the case of a negative position, the water channel is on the right side of the water slide (in downstream direction) locations between the transitions, the friction coefficient was linearly interpolated, so no large jumps occurred. Figure 14 shows the friction coefficient f on the one hand from the sectional optimization and, on the other hand, from the optimization of the three friction parameters described above. Furthermore, the positive partial derivative of the curvature of the centerline is shown whereby it can be seen that the peaks of the partial derivative of the curvature of the centerline are almost similar to the peaks of the results of the sectional optimization. It must be noted that this study has limitations, as it is not possible to optimize global behavior by local cost functions alone, but it is useful to unveil some basic tendencies. Fig. 15 Influence of the number of measured cross-section to F 1 and the resulting friction parameters f 1 , f 2 , and f 3 . The measured cross-sections were counted from the end in group 1 (dark grey) and from the beginning in group 2 (light grey) of the water slide To assess whether the number of used measured cross-sections is sufficient and thus produces meaningful results, an optimization with a lower number of measured cross-sections was carried out while keeping all velocity measures. This was done using optimization problem (43) with w 1 = 0.8 and n 1 from 1 to 26. Two groups were formed: group 1 counting from the end and group 2 counting from the beginning. To compare the results, the determined friction parameters were used to compute F 1 for all measured cross sections. The results are shown in Fig. 15. We can recognize that by using 10 measured cross-sections or more a reasonable result is received.
To get an impression of how the results are affected by small changes of the three friction parameters, a sensitivity analysis was carried out. The results are shown in Table 4. Increasing the friction parameters leads to a better position determination but to a worse velocity determination and vice versa when decreasing the friction parameters. The change of each friction parameter produces a considerable change of results, from which it follows that it is important to consider all of them.
As shown in Eq. (32), the Saint-Venant correction factor κ = N∂h/∂ζ leads to two correction terms for the calculation of the tangential acceleration. The first correction term N(Rα − b)Ṙ/(m w bw) takes into account the changes in the radius of the water-channel cross-section, and the second correction term N A/(m w bw 2 ) takes into account the shallow water. To determine the influence of both correction terms, the location of the middle of the water channel in the water slide was simulated with both correction terms, without radius correction, and without any correction terms. Figure 16 shows the results of all three simulations. No remarkable differences can be seen, which indicates that the consideration of the Saint-Venant correction is not necessary and the water streamline can be found by a simple Bernoulli streamline model with friction losses.

Comparison between the computed water streamline and a CFD-simulation
The water channel was also simulated with OpenFOAM [16] using the turbulent Reynolds Averaged Simulation (RAS) method with the two-equation k-omega Shear Stress Transport (SST) model and the interFoam solver. The PIMPLE algorithm was used with one  . 16 Location of the middle of the water channel using both Saint-Venant correction terms (black), Saint-Venant correction term without radius correction (dashed dark grey), and no Saint-Venant correction (light grey). The position is computed by subtracting π/2 from the transversal coordinate v C and multiplying the result by the water-slide radius r. A position of 0 m in the cross-section indicates that the water channel is on lowest position in the cross-section (transversal coordinate v C is π/2). In the case of a negative position, the water channel is on the right side of the water slide (in downstream direction)

Fig. 17
Location v OF i of the middle of the water channel as computed from the OpenFOAM simulation at cross-section i nonorthogonal corrector and a maximum of 25 outer correctors, where the outer corrector residual control for the velocity and pressure was set to relTol = 0.1 and tolerance = 1e-4 and to relTol = 0 and tolerance = 1e-4 for the final step. To compute the water-channel boundaries at the different cross-sections, the solution was sliced into thin water slices with a thickness of 0.05 m. Each slice was assigned to a cross-section, and the water-channel boundaries were computed using only the cells contained in the slice filled with at least 20% water, by projecting the cell centers in the yz-plane of the corresponding frame K T , as shown in Fig. 17. For each cross-section, the location v OF i of the middle of the water channel was computed as v OF i = π − 0.5 acos  In total, three simulations in OpenFOAM were performed on a processor Intel(R) Core(TM) i7-9700 CPU @ 3.00 GHz using eight cores in parallel. The differences between the three simulations (OF1, OF2, OF3) and the corresponding CPU-times are shown in Table 5. All three simulations were stopped after 22.5 s simulation time. Figure 18 shows the results of the OpenFOAM-simulations and the surface-joint approach compared to the measurements (black crosses), parameterized in terms of the middle of the water channel in the water-slide cross-sections. We can see that the simulation using the surfacejoint approach shows better results compared to the measurements using the simulation in OpenFOAM.
To simulate in OpenFOAM with an acceptable amount of time, a rough mesh was used. First "blockMesh" was performed with cells of length and width 0.05 m and height 0.04 m. To get the contour of the water slide, "snappyHexMesh" was executed with two refinement steps for the inlet and one or two refinement steps for the water slide. A finer mesh and setting the courant numbers to one could increase the quality of the simulation but would require an even larger amount of CPU time.

Parameter identification for the sliding-person model
A water-slide ride was measured using the same video cameras and marks described in Sect. 3.1, but only 26 of the 28 measured cross-sections could be correctly processed. At five straight sections the travel velocity of the sliding person w M SP ,i in the middle of the section was estimated by dividing the length of the section by the time it takes the sliding person to pass through it. As the position of the sliding person was recorded with a drone on most measured cross-sections, several rides were evaluated to get the searched positions. The person was instructed to start always in the same way and to keep the arms and legs crossed during the ride. The footage of an additional fixed camera showed no major differences for the different rides. The initial condition of the position was determined on the basis of the measured cross-section 0, and the initial condition of the velocity was identified using a straight section located directly in front of the measured cross-section 0. To compare the simulated results to the measurements, the positions M P i , M L i , and M U i of the pelvis, lower and upper body on each measured cross-section i, respectively, were computed from the video material, as shown in Fig. 19 exemplarily for M L i , and then divided by the radius of the water slide r to obtain the angles v M P i , v M L i , and v M U i of the three body parts.

Method
A parameter identification was carried out for the multibody model presented in Sect. 2.5 using numerical optimization. The corresponding nonlinear optimization problem was formulated as minimize with the friction coefficient μ D (see Eq. (37)), set equal for all three contact points, as optimization parameter. In the cost function, v C P i , v C L i , and v C U i are the computed cross-section angles of the contact points of the passenger pelvis and lower and upper bodies at crosssection i, respectively. Similarly, w C SP ,i is the computed travel velocity of the sliding person, and w 2 is the corresponding weighting factor. The constant factor k 2 = 1 rad s/m is a rough approximation for the ratio of the expected range of the computed angles (here approximately 2 rad) to the expected range of the computed velocities (here approximately 2 m/s) and hence renders velocity costs in rad 2 . Overall, the optimization uses n 2 = 25 sampling points for cross-section angles and m 2 = 5 sampling points for travel velocities, and the optimization was carried out subject to lower and upper bounds l 2 and u 2 for the friction coefficient. The optimization problem was solved with the nonlinear optimization routine nag_opt_nlin_lsq (e04unc) [17], a sequential-quadratic programming least-squares solver, with the tolerances listed in Table 6. The constraint equations were solved at position level with the routine HYBRJ, which is part of MINPACK [19], a Newton-solver that uses a modified Powell hybrid method, with tolerances set to xtol = 1.0e-8, ftol = 0, and factor = 100. The integration of the equations of motion was performed with the routine LSODAR, which is part of ODEPACK [18], a Livermore solver for ordinary differential equations with automatic method switching for stiff and nonstiff problems and root finding functionality, with tolerances relTol = 1.0e-10, absTol = 1.0e-9, and output time-step 0.005 s. The root finding functionality was used to find the exact moment in time at which the Darboux frames enter and leave the water. The results of the parameter identification with different weighting coefficients w 2 are shown in Table 7.

Fig. 20
Position of pelvis, lower body, and upper body, and velocity of the sliding person with and without drag forces in the u-direction in comparison with the measurement (black crosses). The position is computed by subtracting π/2 from the transversal coordinates v C P , v C L , and v C U , respectively, and multiplying the result by the water-slide radius r. A position of 0 m in the cross-section indicates that the corresponding Darboux frame is on the lowest position in the cross-section (transversal coordinates v C P , v C L , and v C U are π/2, respectively). In the case of a negative position, the Darboux frame is on the right side of the water slide (in downstream direction) Finally, a weighting coefficient w 2 = 0.6 was chosen, because here the difference between position cost and velocity cost is the smallest, with the velocity cost smaller than its position counterpart. The results are shown in Table 8 and Fig. 20. In addition, Fig. 20 shows the influence of the drag force (see Sect. 2.5 and Eq. (38)) on the resulting sliding-person trajectory. Once the parameters are identified, a simulation run of the 15.8-s-long ride requires a CPU time of 10.98 s. We can see that the computed motion matches the measurements fairly well, considering that the computation is performed faster than real time.  . 21 Influence of the number of measured cross-sections to F 2 and the resulting friction parameter μ D .
The measured cross-sections were counted from the end in group 1 (dark grey) and from the beginning in group 2 (light grey) of the water slide

Discussion
To get an impression of how the results are affected by small changes of the friction parameter, a sensitivity analysis was carried out. The results are shown in Table 9. Increasing the friction parameter leads to a better position determination but to a worse velocity determination and vice versa when decreasing the friction parameter. The change of the friction parameter produces a change of the results with a gain factor of 6 to 10. Thus friction plays an important role in simulation prediction quality.
To assess whether the number of used measured cross-sections is sufficient and thus produces meaningful results, an optimization with a lower number of measured cross-sections was carried out while keeping all velocity measurements. Two groups were formed: group 1 counting from the end and group 2 counting from the beginning. To compare the results, the determined friction parameter was used to compute F 2 for all measured cross sections. The results are shown in Fig. 21. We can see that eight measured cross-sections or more lead to reasonable results.
To get an impression of the influence of the body mass, the body length, the height of the center of mass above the water slide, and the radii of the spheres used to calculate the drag forces acting on the sliding person, additional simulations with variations of these Fig. 22 Influence of the mass and body length of the trajectory of the sliding person. The position is computed by subtracting π/2 from the transversal coordinate v C P and multiplying the result by the water-slide radius r. A position of 0 m in the cross-section indicates that the corresponding Darboux frame is on the lowest position in the cross-section (transversal coordinate v C P is π/2). In the case of a negative position, the Darboux frame is on the right side of the water slide (in downstream direction) Fig. 23 Influence of the height of the center of mass and radii, used to calculate the drag forces acting on the sliding person. The position is computed by subtracting π/2 from the transversal coordinate v C P and multiplying the result by the water-slide radius r. A position of 0 m in the cross-section indicates that the corresponding Darboux frame is on the lowest position in the cross-section (transversal coordinate v C P is π/2). In the case of a negative position, the Darboux frame is on the right side of the water slide (in downstream direction) parameters were carried out. Figure 22 and Fig. 23 show the resulting pelvis trajectories. The results of the original simulation are shown with a black solid line. We can see that the most sensitive parameters are the heights of the centers of mass with respect to the surface. We Fig. 24 Results of the position and velocity of the water channel compared to the measurement using only static friction, static and curvature-depending friction, and all three friction parameters. The position is computed by subtracting π/2 from the transversal coordinate v C and multiplying the result by the water-slide radius r. A position of 0 m in the cross-section indicates that the water channel is on the lowest position in the cross-section (transversal coordinate v C is π/2). In the case of a negative position, the water channel is on the right side of the water slide (in downstream direction) can also infer that though the drag forces play an important role and need to be considered, the values used for the sphere radii are not critical.
To show the influence of each of the friction parameters on the prediction of the water channel, the water channel was additionally computed using only static friction (with f 2 = 0 and f 3 = 0) as well as static and curvature-depending friction (with f 3 = 0). The results of these computations are shown in Fig. 24 and Table 10 together with the results obtained with the full friction model and the measurements. We can see that the full friction model yields the best results. Additionally, the different water-channel predictions were used to model the dynamics of the sliding person. The results of these additional simulations are shown in Fig. 25 and Table 11. We can see that an accurate water-channel computation plays an important role in the simulation of the dynamics of the sliding person and should not be neglected.

Fig. 25
Pelvis position and velocity with only static friction, static and curvature-depending friction, and all three friction parameters in comparison with the measurement (black crosses). The position is computed by subtracting π/2 from the transversal coordinate v C P and multiplying the result by the water-slide radius r. A position of 0 m in the cross-section indicates that the corresponding Darboux frame is on the lowest position in the cross-section (transversal coordinate v C P is π/2). In the case of a negative position, the Darboux frame is on the right side of the water slide (in downstream direction)

Conclusions
The computational framework presented in this paper allows for a very fast and sufficiently accurate computation of the dynamics of water-slide rides, which could be used to approach the fun/safety tradeoff that water slide design is constantly challenged with in affordable computational times. The proposed framework uses curve fitting with B-splines to model the differential geometry of the water-slide surface with circular cross-section, minimal coordinates and surface joints to model the dynamics of a kinematical chain sliding on the water-slide surface, stream-filament theory, and Saint-Venant equation (which has a very small effect on the results and thus can be neglected compared to the Bernoulli streamline model with friction losses in future work) to model the water-streamline, and friction and drag coefficients to model the interaction between water channel and sliding objects. The paper shows that sufficiently accurate and very fast computations are possible for real-time simulation and layout optimization purposes.
Future work includes extending the framework to handle water slides with more than one inlet, general water-slide cross-sections, and riders with different sliding techniques and sliding aids like tires. Furthermore, different water slides should be investigated to assess how robust the identified parameters are.