MATHEMATICAL MODEL FOR RIVERBOAT DYNAMICS

Present work describes a simple dynamical model for riverboat motion based on the square drag law. Air and water interactions with the boat are determined from aerodynamic coefficients. CFX simulations were performed with fully developed turbulent flow to determine boat aerodynamic coefficients for an arbitrary angle of attack for the air and water portions separately. The effect of wave resistance is negligible compared to other forces. Boat movement analysis considers only two-dimensional motion, therefore only six aerodynamics coefficients are required. The proposed model is solved and used to determine the critical environmental parameters (wind and current) under which river navigation can be conducted safely. Boat simulator was tested in a single area on the Ljubljanica river and estimated critical wind velocity.


Introduction
Determining resistance and seakeeping for different boat types using a CFD approach is advancing rapidly and many recent studies acknowledge its importance [1,2].Prediction of resistance is the oldest application of CFD in ship hydrodynamics and its accuracy has improved significantly.Simulations show an average error of 3.3%  (displacement) for resistance in both high and low Fr (Froude number) scenarios in recent studies [3].The use of CFD tools is now readily available and puts prediction of resistance and seakeeping within reach for most users.In this study, Ansys CFX package served as the simulation platform for the calculation of force and moment coefficients.
The primary goal of this work is to establish a complete and simple mathematical model of boat dynamics that can be used in the analysis and estimation of riverboat critical environmental parameters for safe river navigation [2,[4][5][6][7].A critical point on the Ljubljanica river (Figure 6) was under investigation.Figure 1 illustrates a typical boat on the Ljubljanica river.Riverboats are flat bottomed boats with a small keel and draft and a large superstructure.The ratio of length to beam compared to sea boats is small.In this study, the wave resistance phenomena is not taken into account in the model, because it is assumed that predominately the calm water prevails and its magnitude is much smaller compared to the drag force.Wave resistance is only relevant in the xB direction (Figure 2) and does not influence dynamics in other directions as friction force (Drag) does.To keep model simple the effect of wave resistance force compared to other friction forces can be neglected.External forces such are wind and current are incorporated into the model by using the squared resistance law.CFD simulations served to incorporate the effect of the apparent direction of external fluid flow on boat direction into the model (Figure 2) and the boat aerodynamic coefficients were obtained for rotated boat geometry.For the upper portion of the boat (air), boat geometry was rotated by 180 o , and for the underwater portion of the boat, boat geometry was rotated by 90 o .Afterwards, a single spot on the Ljubljanica river called Špička was investigated in order to estimate critical wind velocity.

Model
The mathematical model of boat motion makes a few assumptions: (a) the boat is rigid body and oscillations are ignored, (b) the boat is symmetrical with respect to its xB axis (Figure 1), (c) boat motion is restricted to two-dimensional motion only, and roll and pitch oscillations are ignored.Restrictions are possible because of the specific boat geometry, as they are much wider than conventional boats.As mentioned in the introduction, the foundation of the model is the square drag law: where F is the force exerted on the body by external fluid flow, M is the moment exerted on the body by external fluid flow, C is an aerodynamic coefficient depending on the relative angle of fluid flow, S is boat frontal cross-sectional area, h is the moment arm from arbitrary  2).The drag force coefficient is an external force component in the direction of the B x axis in the boat coordinate system, the side force coefficient is an external force component in the direction the of B y axis in the boat coordinate system, and the yaw moment coefficient is a moment of external fluid flow acting around the zB axis, with its origin on the boat's centre of gravity (CG).In the simulations of fluid flow motion around the riverboat there is a distinct difference between two cases: the upper structure portions and the underwater structure portion (Figure 1).The upper structure portion is subject to air flow and the inlet simulation velocity was set to 20 m/s.The underwater part is subject to water flow and the inlet simulation velocity was set to 2 m/s.These are the typical maximal average velocities for air and water on the Ljubljanica river.
In order to elaborate on the dynamical model of boat motion, boat kinematic equations describe the motion of about the boat's centre of gravity CG:  boat course with respect to the inertial coordinate frame (x ,y).Velocities can be determined using the following set of Newton-Euler equations: where: ; , , external and propulsion forces in the xB direction in the boat frame.
; , , external and propulsion forces in the yB direction in the boat frame.
; , , , Fx external and propulsion moments in the boat frame, where P x is the distance to CG from the stern midpoint (propeller coordinate).0 J boat inertial moment.
Index A stands for air, W for water, R for rotation, and P for propulsion.All forces and moments are calculated from the apparent flow velocity and direction.Boat motion and environmental flows are merged into apparent flow variables (velocity and direction) called the external force and moment.There is distinction between two types of external effects: air and water.The rotational moment R M requires extra words to clarify its origin.It comes from the boat's rotation about its centre of gravity (zB axis in Figure 1), causing rotational resistance.There is another distinction between the two sides of the boat for the rotational moment.The forward portion (F) is the portion of the boat where the y component of the apparent velocity v ( v y ) has the opposite direction with respect to boat rotation.In this case, the total velocity acting on the side of the boat is the sum of v y and x  , where ω is boat angular velocity and x is the distance from its centre of gravity CG.The backward portion (B) is the opposite side from the F portion, where the apparent velocity acts in the same direction as boat rotation and the total velocity is the difference of v y and x  .In cases with high apparent velocity, the moment can also change signs.Thus, it is a need to distinguish between two cases for the back side.When the y component of apparent velocity is below the maximal tangential velocity ( /2 L


) at the extreme boat bow or stern, and above the maximal tangential velocity, the following set of equations defines the boat rotational moment:  9) where R S C is the side force coefficient at 90 R   .The set of Equations ( 2)-( 4) determines the complete boat motion in two dimensions.In addition to Equations ( 2)-( 4), there must be specified the value of external (apparent wind and current) and propulsion forces, and their moments, according to the apparent flow velocity and direction specified in Eq. ( 1), as explained at Figure 2. In order to determine the apparent velocity, vR, the next equation must be solved (explanation in Figure 2): and the apparent velocity magnitude of external fluid flow is with direction: The environmental fluid direction E  (explanation in Figure 2) is calculated from the true flow direction E  and boat direction  as: Adjustments are made for wind and current direction E  .Wind direction is defined as from where it blows, and current direction is defined as where it acts.The difference from 180 must be taken into account in Eq. (8).External forces and moments in Equations ( 3)-( 4) are calculated using Eq. ( 1), where the aerodynamic coefficients are calculated from Equations ( 9) and (10), with corresponding diagrams in Figures 5 and 6.Apparent velocity ( v R ) is the velocity defined in Eq. ( 6), and apparent angle ( R  ) is defined in Eq. (7).Diagrams in Figures 5 and 6 are obtained via CFD simulations in Ansys CFX system.The underlying turbulent model was set to standard k- Menter's model with tunnel boundary conditions (BC) (inlet, outlet, side, bottom and top).Inlet-outlet BC is standard, bottom BC was set to no-slip and top and sides BC were set to symmetry BC.
Results shown in Figures 5 and 6 illustrate air and water external forces and moments obtained via a fitting procedure [8].The air aerodynamic coefficients have the functional form defined in Eq. ( 9) and the water aerodynamic coefficients have the functional form defined in Eq. (10).The argument x in Equations ( 9) and (10) represents the apparent flow angle defined in Eq. (7), measured in radians.It is defined over the interval   ,   .In Figures 5 and 6, the CFD simulation results with plus symbols used as markers of the same colour as the coefficients they reference are displayed.The form of the fitting functions in Equations ( 9) and ( 10) is chosen such that R 2 , the coefficient for statistical determination, is close to 1.For the air coefficients, asymmetry in the results reflects the fact that the superstructure does not respect rotational symmetry.The underwater boat geometry is nearly rotationally symmetric, and the model functions for water aerodynamic coefficients were to reflect this.
A system of ODEs ( 2)-( 3) was integrated with an explicit Runge-Kutta scheme of orders two and four [8].The time increment was set to 0.01 t  .In time integration the instabilities are observed when environmental effects were strong and the time increments were too high.A time-adaptive scheme would be more appropriate for such a system.

Results
In Section 2, the mathematical model of riverboat dynamics is described with external reactions described by aerodynamic coefficients depending on the external flow direction.In order to calibrate the boat model Equations ( 9)-( 10), one must perform CFD simulations of fluid flow around the desired boat.In this case, the most common shape of a riverboat was constructed that can be found on the Ljubljanica river.An example boat with all its dimensions and a 3D view is presented in Figure 1.Its physical properties are described in Table 1.Fluid flow was calculated for two different boat portions: the underwater portion described by water fluid flow and the upper boat portions described by air fluid flow.In both cases, the boat geometry was rotated in increments of 10 .For the air flow, the angle interval ranged from 0 to 180 , and for the water flow, the angle interval ranged from 0 to 90 .In order to calculate the external forces from Eq. (1) on the cross sectional area, the density and moment arm must be known and can be obtained from Table 1, where all the parameters needed to initialize the boat simulation are presented.
The set of Equations ( 2)-(3) are solved using an explicit Runge-Kutta numerical scheme of the fourth order.Tests of the second order were performed but no obvious differences were observed.In order to have a stable integration procedure, a small time interval must be used when strong environmental flows are present.
Such a system does not require extensive computational resources and can be solved implicitly, or with an adaptive scheme, in order to achieve stability in numerical calculations across the wide range of external action magnitudes.Stability of the present ODE model is currently under investigation for both the explicit and implicit scheme, and will be incorporated into a future model.The algorithm was coded in the C++ programming language and can be obtained from the author upon request.The case study was performed with boat geometry similar to that found in Špička region of the Ljubljanica river (Figure 1).Maximal boat velocity is always set to 5 knots, and is reached within 7 seconds at full ahead power.All simulations are performed at maximal boat velocity.In the calculation, the yaw moment coefficient arm h must be also specified as noted in Eq. ( 1).CFD simulation results already account for the moment arm, providing the actual yaw moment as the result, and h is set to a unit value of 1m.
The first test of the boat simulation algorithm was a turning manoeuvre test, which was performed with varying rudder angles.The results are presented in Figure 7.The typical opposite trajectory motion was observed immediately following the rudder turn in the region [{0, 25}-{0, 35}].The same effect, but magnified, is seen in Figure 8 immediately after the coordinate system origin.It is seen that the turn radius extends from approximately 160m for a rudder angle of 10 °, to 25m for a rudder of angle 40 °.The second test was performed in a domain similar to the Špička region shown in Figure 6, in order to estimate critical wind velocity.The results are presented in Figure 8.On left side of the Figure 8, it can be seen the results for the following simulation scenarios: 12s full ahead, rudder angle 0 °, boat azimuth 45 °; 10s full ahead, rudder angle -30 °; 20s full ahead, rudder angle 0 °.The current velocity was set to 1m/s for all simulations and the wind ranged from 0-60 m/s.Without the intervention of a skipper, the boat could run ashore within a few meters depending on the wind speed.On the right side of Figure 8, the boat throttle was also full ahead, but after a time period of 20s the rudder angle was changed to avoid collision with the shore by as much as possible (skipper intervention).In the case with a wind speed of 60 m/s, the boat manoeuvre was pushed to its limits in order to avoid shore collision.Boat dimensions must be taken into account for this simulation.It can confidently state that the wind velocity limit is approximately 60 m/s in this scenario, but this is as a very strong wind that rarely occurs on the Ljubljanica river.In the event of such strong wind gusts, river traffic should be forbidden by the local maritime authority.

Conclusions
The squared resistance law used in boat motion modelling shows promising results for the prediction of critical environmental reactions for safe river navigation.The proposed model is simple, but contains all the relevant physics responsible for regulating the motion of a riverboat [2,[4][5][6].With the inclusion of wave resistance effects, similar to those shown in [9], the present boat model could be improved to account for higher velocities.The model can be easily extended for arbitrary boat shapes.The only changes required would be to update the parameters defined in Table 1 and in Equations ( 9) and (10).
Our riverboat model incorporates the asymmetric behaviour of the superstructure portion of the boat and the symmetrical behaviour of the underwater portion.Additionally, few complex simulations with very strong winds and currents were calculated.In some cases with strong wind, the boat became unable to turn any further into the wind due to a strong yaw air moment balancing against the propulsion moment.This example demonstrates the extreme conditions a boat may encounter at sea, or perhaps in certain rivers.The inability to turn any further into the wind severely compromised boat manoeuvrability.Integration of the ODE system was performed in an explicit manner with a constant time increment.This is not a comprehensive technique, and it should be replaced with an adaptive time integration scheme in order to achieve simulation stability across a wide range of environmental parameters.
A simulation in an area similar to the Špička region showed that the wind velocity limit is approximately 60 m/s for the type of boat most commonly used in that region.Such wind gusts do occur occasionally on the Ljubljanica river.The use of this kind of simulation tool makes it possible to analyse the environmental parameters required for safe river and sea boat navigation.

Fig. 1
Fig. 1 Riverboat 3D figure with corresponding dimensions and boat coordinate system (xB, yB, zB) centred on the centre of gravity (CG), designed in SolidWorks.

Fig. 2
Fig. 2 External forces with defined angles and relative direction in the boat coordinate system (xB, yB).

Fig. 3
Fig. 3 Diagram of wind vector (vA), current vector (vW), and propulsion vector (vP) with appropriate distances from the vector origin to the boat's centre of gravity (xA, xW, xP).

Fig. 4
Fig. 4 Boat air drag aerodynamic coefficients.CD, CS, and CY are drag, side force coefficient, and yaw moment coefficient. is the apparent wind direction.Cross markers show data obtained from CFD simulations and curves show plots of Eq. (9)

Fig. 5
Fig. 5 Boat water drag aerodynamic coefficients.CD, CS, and CY are drag, side force coefficient, and yaw moment coefficient. is the apparent current direction.Cross markers show data obtained from CFD simulations and curves show plots of Eq. (10)

Fig. 6
Fig. 6 Critical point on the Ljubljanica river -Špička.Intersection of the Grubar channel and Ljubljanica.

Fig. 8
Fig. 8 Wind and current effects on boat motion in a narrow river region.Current velocity was set to 1m/s and wind velocity was changed as shown in the plot legend.Arrows on the path indicate boat direction.

Table 1
Simulation and boat geometry parameters.
Moment arm (set to unit)