Integrated computational model for Sea Organ simulation

An integrated approach for the simulation of hydraulic and musical aspects of the Sea Organ in Zadar, Croatia, is presented in the paper. The computational model combines a real-life sea wave generation algorithm, a numerical model for water mass oscillation in organ pipes, and a musical sound generation model. The results show that, under the same wave conditions, internal water level oscillations differ between individual Sea Organ segments mainly because of different pipe geometries. Furthermore, analyses revealed that the resulting sound depends not only on pipe geometry, but also on the length and direction of waves.


Introduction
Over the last few years, the north-western waterfront of the Zadar promenade in Croatia has been transformed from a neglected periphery to a modern, attractive and well-frequented public space.Two urban installations -the Sea Organ and the Greeting to the Sun -have gained a worldwide acclaim for their architectural, visual and acoustical innovation.The response of general public and professional circles to these two installations has been extremely positive, which has resulted in several highly valuable awards and recognitions [1].However, the design and construction of the Sea Organ proved to be mostly experimental, due to the lack of adequate numerical models.A simplified computational model for simulating the sound of the Sea Organ, based on randomly generated sea waves, has already been presented [2].The simplified model [2] assumes a direct linear relationship between the sea surface elevations and water level oscillations within the pipes.This assumption, however, neglects the influence of the inertial and friction forces and cannot therefore take into account the differences between pipes of different geometries.In other words, only the influence of wave dynamics on the resulting sound was analysed, whereas the influence of the pipe diameter, length, and slope, was neglected.The purpose of this study is to develop and present an improved three-part computational model for simulating the Sea Organ, which includes an algorithm for generation of random sea waves, a one-dimensional numerical model for internal water mass oscillations, and the corresponding sound model.Such a model can then be used to simulate internal oscillations and the resulting sound generated by any sea organ-type structure, and thus help in the decision-making process and in the design of similar structures.The main contributions of this work are three-fold: an original one-dimensional numerical model for solving water mass oscillations inside organ pipes driven by seawaves is developed an integrated model, uniquely blending the numerical hydrodynamic model with a computational sound model, is presented hydraulic and musical aspects of the Sea Organ in Zadar are examined using the proposed model to detect and define the main effects of different pipe-system components and wave parameters on internal oscillations and the resulting sound.
The paper is organized as follows: The hydraulic and musical characteristics of the Sea Organ in Zadar are first described and illustrated; this is followed by description of an integrated model consisting of the computational model for spectral waves, the numerical model for water level oscillations, and the conceptual waveform model for sound generation; two sets of results are then shown and discussed, while main conclusions and future work are presented in the final part of the paper.

Hydraulic and musical characteristics of the Sea Organ
The Sea Organ was built in 2005 in the scope of reconstruction of the existing sea-wall at the Zadar promenade, which was converted into a descending stepped wall [1].This 75-m long instrument is divided into seven segments, and each segment is defined by a different number of steps and the corresponding crest height.A total of 35 pipes of various lengths and diameters are positioned underneath the stone steps [1].The horizontal distance between pipes in each segment is ~1.5 m and the distance between two adjacent segments is ~4.5 m (Figure 1).The hydraulic pipe system consists of three distinct parts (Figure 2): a horizontal entry pipe positioned below the sea surface an acoustical pipe with the labium positioned inside the resonating corridor, which is built under the promenade a sloped smaller-diameter pipe connecting the lower entry pipe and the upper acoustical pipe.
The entry and sloped pipes are made of polyurethane (PE), whereas the acoustical pipes are made of stainless steel [1].The sound generation process can be described as follows: the wind, waves, tides and passing boats initiate the movement of the sea surface in front of the organ seawall; the vertical movement of the sea surface initiates water oscillations in sloped pipes; the water movement then pushes the air into the acoustical pipe and out through the labium, which finally produces a sound of a predefined frequency.The sound is released through special openings into the surrounding area, where it can be heard and enjoyed by the audience.In the musical sense, each segment is tuned either to G-major chord (D-G-d-g-h) or C-major chord with an additional sext (C-Dc-e-a), as illustrated in Figure 1.All tones correspond to the frequencies ranging from 60 to 250 Hz [3].The choice of these

Methods
The proposed integrated model consists of three parts: a computational model for spectral waves a numerical model for water mass oscillations in organ pipes waveform model for sound generation.
The computer code for all three parts is written in MATLAB [4], without any external libraries other than matrix2midi [5] and writemidi [5] for generation of the final MIDI sound file.

Computational model for spectral waves
A computational model for simulating sea surface elevations is based on the spectral description of wind-generated waves.Here, the methodology for the phase-amplitude model and a two-dimensional frequency-direction spectrum is used, as proposed in [6].A harmonic linear wave propagating in the x-y horizontal plane, in the main direction θ relative to the negative x-axis, is initially considered [6]: where θ is the surface elevation at time t, a is the wave amplitude, ω is the angular frequency, k = 2π/L is the wave number, L is the wave length, and α is the wave phase.To account for randomness of real waves, and to describe their stochastic nature, an irregular surface elevation at a given x-y point in space is computed as the sum of a finite number of harmonic waves, defined by different wave amplitudes and phases [6]: where N is the finite number of frequencies (denoted by index i), and M is the finite number of wave directions (denoted by index j).Each individual harmonic wave in Eq. ( 2) has a unique amplitude a i,j , derived from a given wave density spectrum and the random phase , which is uniformly distributed [4].Note that the wave number k i is directly related to the wave frequency ω i , and therefore shares the same index i.The amplitude for each harmonic wave component is defined by a discrete wave spectrum as , where spectrum S(ω) is discretized by a finite number N of frequency increments .Several wave spectrums have been derived from field observations, such as the JONSWAP [7] or Pierson-Moskowitz [8].However, the Tabain spectrum [9], as derived from observations in the Adriatic Sea, was used in this study: (3) where g is the acceleration of gravity, H s is the significant wave height, and γ is the characteristic spectrum width parameter.The value γ = 1.63 is often recommended for the Adriatic Sea, while the function p(ω) is defined as follows [9]: (4) whereσ(ω) is the spectrum shape parameter depending on the frequency.It is defined according to the peak frequency ω p , as follows [9]: Note that the Tabain spectrum requires only the significant wave height to define the surface elevation.Although the Tabain spectrum has been selected here because it is the most realistic one for generation of the Adriatic Sea waves, other spectrums may also be easily implemented in this computational model.A two-dimensional wave density spectrum defines the surface elevation distribution over a range of frequencies and directions [6].Therefore, to include the variability of wave propagation directions, the one-dimensional Tabain spectrum defined by Eq. ( 3) is extended by an additional parameter D(θ), so that S(ω, θ) = S(ω)D(θ).The directional distribution D is defined as shown below [6]: Nino Krvavica, Igor Ružić, Nevenka Ožanić where c = Γ (s+1)/[Γ (s+1/2)2√π], Γ(s) is the Gamma function, s is the parameter that controls the directional span, and θθ 0 is the angle between any direction θ and the main direction θ 0 The wave amplitude that varies over the frequency and wave direction is defined as .To illustrate the influence of the directional span on the sea surface elevation, directional spectrums and the corresponding sea surfaces at an arbitrary moment are shown in Figure 3 for two different s parameters and H s = 0.6 m.The wave number may be defined by a dispersion relation , where d is the water depth.However, this equation is nonlinear and requires iterative procedures, such as bisection or Newton-Raphson.It would therefore be preferable to use an explicit equation, such as [10]: The reflection from the organ sea-wall is accounted for by a local increase in wave amplitude.Non-linear wave-wave interactions are considered negligible and are omitted at this stage of model development.Therefore, surface elevations at points near the sea-wall are defined as follows: (8) where the reflected wave amplitude is a r = K r a, K r is the reflection coefficient, and is the incident wave amplitude.The coefficient of reflection from a vertical structure is determined based on the crest height h c above the sea water level (SWL) [11] and the main angle of wave propagation θ 0 , as follows: (9) The wave-induced pressure at a given depth h is defined, according to the linear wave theory, as the sum of hydrostatic and dynamic pressures [6]: (10) The hydrostatic pressure depends only on the depth h.On the other hand, the dynamic pressure is time-dependent, and it can be calculated as the sum of dynamic pressures induced by each harmonic wave component [6]: The computational implementation of the random wave algorithm for generation of sea surface elevations is relatively simple.First, a spatial domain is defined by the length L X and Integrated computational model for Sea Organ simulation width L Y .Next, the spatial step Δl is defined, so that a matrix J x K can be constructed, where J = L X / Δl is the number of rows and K = L Y / Δl is the number of columns.Also, the temporal domain is defined by the total simulation time T and time step Δt.Therefore, we have T / Δt matrices of size J x K. Finally, each element of each matrix is computed according to Eq. ( 2), and the dynamic pressure is computed by Eq. ( 11) at predefined coordinates corresponding to each pipe entry location.

Governing equations
Wave-induced water oscillations in organ pipes are computed by a one-dimensional numerical model, which is based on the continuity equation and the time-dependent Bernoulli equation.
The continuity equation is written as [12]: where Q is the flow rate, V is the internal volume of water, and t is the time.Under the assumption of incompressible flow, the time-dependent Bernoulli equation, derived by integrating the Euler equation along a streamline between arbitrary points 1 and 2, is written as follows [12]: where z 1,2 is the water depth, p 1,2 is the static pressure, ρ is the fluid density, v 1,2 is the water velocity, and is the energy loss.The last term on the right-hand side of Eq. ( 13) is energy per unit mass required to change the velocity v in an infinitesimal part ds of a streamline.Equations ( 12) and ( 13) are applied to study water oscillations in the sea organ pipe (Figure 4).Under the assumption that the water level oscillates in the second pipe only, i.e., that the first pipe is always submerged, Eq. ( 12) can be rewritten as: (14) where l 2 is the length of the wetted part of the second pipe, and is the cross-sectional area of the second pipe.
Considering the streamline that runs along the centre of an organ pipe, the left-hand side of Eq. ( 13) is replaced by the total wave-induced pressure p wave at depth h corresponding to the centre of the pipe.The pressure in the pipe changes with the water level h p , and hence p 2 = phh p = pgl 2 sinα, where α is the angle of the second pipe with respect to the horizontal plane (Figure 4).The acoustical part of the sea organ is an open-pipe system with a small opening in front of the labium; therefore, the air pressure in the pipes is considered negligible.The change of the water velocity in each segment of the pipe system is assumed to be instantaneous.Hence, Eq. ( 13) is rewritten as follows: ( where v 1 and v 2 are the water velocities in the entry pipe and sloped pipe, respectively, is the length of the entry pipe, and is the energy head loss.Since the water velocity depends on the pipe diameter, the equation is rewritten in terms of the water flow rate Q = A 1 v 1 = A 2 v 2 , as follows: (16) where A 1 and A 2 are the cross-sectional areas of the entry pipe and sloped pipe, respectively.To ensure the correct sign for energy dissipation, the squared flow rate is written as the product of flow rate and its absolute value.Finally, the timedependent Bernoulli equation is rewritten as follows: (17) while the total loss coefficient β is defined as: (18) where ξ E is the loss coefficient at the pipe inlet, ξ A is the loss coefficient at the pipe elbow, ξ R is the loss coefficient due to pipe profile reduction, whereas the last two terms in Eq. ( 18) are the pipe friction losses defined by the Darcy-Weisbach equation Nino Krvavica, Igor Ružić, Nevenka Ožanić [12], where D 1 and D 2 are the pipe diameters and and are the friction coefficients, usually computed by the White-Colebrook equation [12].The governing equations are therefore defined as a system of two non-linear ordinary differential equations (ODEs) with two unknowns: the length l 2 and the flow rate Q.Unfortunately, analytical solutions for Eq. ( 14) and Eq. ( 17) are unavailable, and so some numerical method should be applied to approximate the solution.

Numerical scheme
A general approach for solving dynamic response of any system involves direct numerical integration of ODEs.This approach is based on satisfying governing equations at discrete points in time, with a known initial solution [13].There are many different numerical methods and all of them can be categorized as either explicit or implicit.In this work, the implicit trapezoidal rule [13] was used to numerically evaluate the governing system of ODEs.The proposed trapezoidal rule for solving ODE of the form is written as [13]: where y n+1 is an unknown value of function y(t) at the time n+1 = t n + Δt, y n is a known value at the time t n , f(t n , y n ) i f(t n+1 , y n+1 ), and y(t) are values for the time derivative of the function at the respective time t n i t n+1 , while Δt is the time step.The trapezoidal rule is the second-order accurate and A-stable numerical method [11].The proposed numerical method was applied to evaluate Eq. ( 14) and Eq. ( 17) as follows: As the resulting equations are non-linear, an iterative method should be used.A quasi-Newton method is preferred since analytical formulation for the Jacobian matrix of the governing system is non-trivial (mainly due to derivatives of the empirical friction equation), The Broyden method [14] was chosen here to solve the non-linear system of equations.In this iterative method, the Jacobian matrix is replaced by an approximation, which is then easily updated at each iterative step [14].Furthermore, the amount of computations at each step is reduced, and its convergence is superlinear [14].

Conceptual waveform model for sound generation
A simplified sound generation approach is considered in this subsection.It is assumed that each acoustical pipe can produce a single tone of a corresponding frequency, depending on its geometry and acoustical properties.However, the sound amplification depends on the velocity of air flowing out of the pipe.The second assumption is that the air velocity is directly proportional to the computed internal water level oscillations.The first step is to compute the water level oscillations at each pipe, as described in the previous subsection.Once h p is known, the rate of water level change is calculated as the following discrete derivation: Assuming that the sound is generated only when the air flows out of the pipe, the duration of the sound is defined as the time distance between two neighbouring zero-crossing points of Eq. ( 23).The sound wave amplitude is considered to be proportional to positive values between two zero-crossing points.Therefore, the sound wave (waveform) is generated using the following equation: where A is the sound wave amplitude proportional to the rate of change in water level, is the predefined frequency of the sound, and is the time between two neighbouring zero-crossing points.

Results and discussion
Two examples are presented to demonstrate performance of the proposed model.In the first example, internal water level oscillations are compared between different segments, but under the same wave conditions.The wave spectrum was generated based on the significant wave height H s = 0.5 m, with no directional spreading θ = 0°, and with the direction of propagation normal to the coast line θ 0 = 0°.The reflection coefficient was computed by means of Eq. ( 9), based on different crest heights for each segment.Internal water oscillations were computed using the proposed numerical model (Section 2.2).
The pipe geometry at each segment was estimated based on the given design plans (Figures 1 and 2); D 1 = 30 cm for the entry pipe, and D 2 = D 3 = 12 cm for the sloped and acoustic pipes; the entry pipe is positioned at level -0.5 m a.s.l., and the upper acoustic pipe is positioned at +1.6 m a.s.l.; pipe lengths Integrated computational model for Sea Organ simulation Since the pipe material is either PE or steel, friction losses were omitted, whereas minor hydraulic losses were computed using the coefficient selected from known tabled values or formulas from the literature [12].Results from the first example -external and internal water oscillations at different segments -are illustrated in Figure 5.
Internal oscillations generated by the same waves are slightly lower in pipes 1 and 6 when compared to other pipes, which can be explained by smaller reflection coefficient at the first two segments due to lower crest.More importantly, a different response of the internal water mass is evident between the segments.It seems that because of longer pipes (and larger initial volumes of internal oscillations are more damped in pipes 1 and 6 as compared to pipes 11 and 16, and especially when compared to pipes 21 and 26. The second example shows internal water oscillations and the resulting sound under two different wave conditions.The first scenario involves a realistic case of a moderate breeze from the NW (WNW) direction, which results in a significant wave height H s = 0.5 m, with directional spreading θ 0 = 20°, and the main direction of propagation θ 0 = 45°.The second scenario involves a fresh breeze from the SE direction, which results in a slightly higher significant wave height H s = 0.7 m, with the directional spreading θ = 20°, and the main direction of propagation θ 0 = -60°.Internal water oscillations were computed using the same input parameters as in the previous example.Results for the second example are shown in Figure 6 (for the NW scenario) and in Figure 7 (for the SE scenario).The importance of using a numerical model for simulating internal water oscillations is evident in the discrepancy between η and the corresponding h p , for both scenarios (Figures 6.a and 7.a).Furthermore, the influence of wave heights is noticeable when Figure 6.a and 7.a are compared; larger wave heights result in larger internal oscillations.Similarly, larger wave heights also result in louder tones (Figure 7.c compared to Figure 6.c).Pianoroll of the resulting MIDI file shows the influence of the wave direction on the sound melody produced by the organ; NW direction waves generate a melody with an ascending sequence of tones, whereas SE direction waves generate a melody with a descending sequence of tones.

Conclusion
An integrated three-part computational model for simulating hydraulic and musical aspects of the Sea Organ in Zadar, Croatia, is presented in this paper.Random waves were generated by an algorithm based on a directional wave spectrum.Water oscillations in pipes were simulated by a one-dimensional numerical model based on an implicit solver for the governing system of two ordinary differential equations.Finally, the sound and melody were generated by a simple waveform sound model.The results revealed several useful and interesting aspects of the Sea Organ.The geometry of the pipe system, mainly the pipe length and slope, controls the internal water oscillations and, consequently, the resulting intensity and duration of sound.The analysis showed that the first two segments, characterized by longer pipes and lower reflection coefficient, exhibit smaller internal water oscillations and weaker sound when compared to other segments.This finding is in agreement with the author's practical experience, according to which the first segment Integrated computational model for Sea Organ simulation produces almost no sound at all during mild wave conditions.In addition to a strong relationship between the wave height and sound amplification, the analyses also showed that the resulting melody is also influenced by the wave length and direction; northern waves generate a melody with an ascending sequence of tones, whereas southern waves generate a melody with a descending sequence of tones.Although some physical processes, such as the air pressure variations in pipes, non-linear wave transformations, and wavewave interactions, were assumed to be negligible, it should be noted that main processes were captured and simulated by the model quite realistically.Unfortunately, the comparison of the final sound produced by the model cannot be adequately presented in paper format.Nevertheless, the readers are invited to evaluate for themselves the similarities between the proposed model [15] and the sound recording from the Sea Organ in Zadar [16].Future plans include experimental analyses of the physical sea-organ model that will be constructed in the Hydraulic laboratory of the Faculty of Civil Engineering at the University of Rijeka.The focus of the experimental study will be placed on the influence of the labium -a narrow opening in the acoustical pipe -and its effect on pressure and water level oscillations in the pipe system.Based on these results, the proposed model will be thoroughly validated and additionally improved if needed.

Figure 1 .
Figure 1.Ground plan of the Sea Organ in Zadar[1]

Figure 2 .
Figure 2. Two typical cross-sections of the Sea Organ in Zadar [1]

Figure 3 .
Figure 3. Direction spectrum and generated sea surface elevation for H s = 0.6 m and: a) s=262, b) s=65

Figure 4 .
Figure 4. Typical longitudinal section of an organ pipe system

Figure 5 .Figure 6 .
Figure 5. Numerical results for external (sea-water level) and internal (water level in the pipe) oscillations at different segments

Figure 7 .
Figure 7. Results for SE scenario: a) external and internal water level oscillations for pipe 16, b) rate of water level change in pipe 16, c) generated sound (waveform) in pipe 16, d) generated piano-roll at segment 4