Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 41, 1, pp. 119-135, Warsaw 2003 DIFFRACTION LOADS ON MULTIPLE VERTICAL CIRCULAR CYLINDERS BY LEAST-SQUARES TREFFTZ-TYPE FINITE ELEMENTS Małgorzata Stojek Department of Computational Methods in Metallurgy, University of Mining and Metallurgy e-mail: ghstojek@agh.edu.pl Thepaper dealswith the application of the so-called”frameless”Trefftz- type finite elements to the calculation of wave induced loads on the of- fshore, seabed, cylindrical structures. Themethod is based on the use of a suitably truncated T-complete set of Trefftz functions over individual subdomains linked bymeans of a least square procedure. The Neumann homogeneous boundary conditions and the Sommerfeld radiation con- dition are readily incorporated in the trial functions of special purpose finite elements. Key words: Helmholz equation, Trefftz, FEM 1. Introduction Linearized diffraction theory for vertical cylinders extending from the im- permeable seabed to the free surface in water of constant depth leads to the Helmholtz equation in the exterior domain accompanied by the Sommerfeld radiation condition at infinity and the Neumann boundary condition on the wetted cylinder surface. Traditionally, such a problem is reformulated as a bo- undary integral equation over the surface of thebody (Colton andKress, 1992) giving rise later on to an appropriate variational formulation and a boundary element (BE) approximation. The BE discretization can be done directly on the wetted body or some auxiliary smooth surrounding surface. In the latter casewehave the coupling of finite andboundaryelementswhich, althoughmo- re expensive, avoids problems with the integration of singular kernels arising fromBE approximations on non-smooth boundaries (corners, edges). 120 M.Stojek In both casesmentioned above, the application of theBEapproaches turns out expensive formediumand largewave numbers (Harari andHughes, 1992). An interesting alternative seems to be, therefore, the finite element method with two different approaches to an infinite domain. The first approach in- troduces an artificial boundary to make the computational domain finite (for overview see Givoli et al., 1997), the second one is based on the use of infi- nite elements (Bettes, 1992; Gerdes, 1998; Harari et al., 1998). However, the common use of the ”traditional” polynomial shape functions results in the so- called pollution error (Babuška and Sauter, 1997) related to the deterioration of stability of theHelmholtz variational form at largewave numbers. To avoid this, several non-standardFE formulations have beenproposed in the literatu- re. In thesemethods, stabilization is either attempted directly bymodification of the differential operator or indirectly, via improvement of approximability by the incorporation of particular solutions into the trial space of FEM (for overview see Ihlenburg, 1998; Farhat et al., 2001). The present paper shows the application of least-squares Trefftz-type ele- ments (Stojek, 1998) to the practical calculations ofwave induced loads on the offshore, seabed, cylindrical structures. Such T-elements belong to the non- conforming, knowledge-based finite elements andhave been developed directly for the scattering problems and long waves. Their main idea consists in the use of suitably truncated T-complete (Herrera and Sabina, 1978) solutions to the governing differential equation. The prescribed boundary conditions and the interelement continuity are enforced by means of a least-squares proce- dure. The Neumann homogeneous boundary conditions and the Sommerfeld radiation condition are readily incorporated in the trial functions of special purpose finite elements. The paper is organized as follows. The next section shortly recalls the theoretical formulation of the least-squares T-type elements. Section 3 sum- marizes the linear diffraction theory and the wave induced loads on vertical cylinders. Section 4 presents the numerical results for the interference effects among 2, 3 and 4 circular cylinders. Finally, the concluding remarks are given in Section 5. 2. Review of least-squares Trefftz-type finite elements for the Helmholtz equation We consider the following complex valued boundary problem in 2D (Stojek, 1998) Diffraction loads on multiple vertical circular cylinders... 121 ∇2ϕ+k2ϕ=0 inΩ ∂ϕ ∂n =0 on the wetted surface ΓB lim r→∞ √ kr (∂ϕs ∂r − ikϕs ) =0 at infinity Γ∞ (2.1) where k is a wavenumber, ϕ and ϕs denote the two-dimensional unknown total and scattered velocity potentials, respectively, ϕi =ϕ−ϕs is the given incident potential, Ω stands for the unbounded exterior domain occupied by the fluid, ΓB ∪Γ∞ = ∂Ω, and the boundary condition at infinity is the Som- merfeld radiation condition. We recall that the above physical problem may be formulated both in terms of the total or scattered potentials. In our case, the former is preferred since it leads to the homogeneous boundary conditions on the wetted surface ΓB. Let us divide the whole domain Ω into n subdomains Ωk provided with respective local coordinate systems (rk,θk). Over each Ωk we represent an approximate solution ϕk as ϕk =ϕ i+ϕsk =ϕ i+ m ∑ j=1 Nkjckj =ϕ i+Nkck (2.2) or ϕk = m ∑ j=1 Nkjckj =Nkck (2.3) where Nk is a truncated complete set of local solutions to Helmholtz’s equ- ation in Ωk (so called T-functions) and ck denotes a vector of undetermi- ned complex valued coefficients. The corresponding outward normal velocity vk =∇ϕk ·n on ∂Ωk is vk = v i+Tkck (2.4) or vk =Tkck (2.5) The Neumann boundary conditions on ΓB and the continuity in potential ϕ and normal derivative v on all subdomain interfaces ΓI are enforced by minimizing the following least-squares functional I(ϕ) = I(c)=w2 ∫ ΓB |v|2 dΓ + ∫ ΓI |ϕ+ I −ϕ− I |2 dΓ +w2 ∫ ΓI |v+ I +v− I |2 dΓ (2.6) where superscripts +,− in the integrals along ΓI designate the solutions from two respective neighbouringTrefftz fields and w is somepositiveweightwhich, 122 M.Stojek in general, may vary with the segment of matching. Note that the integral along Γ∞ in (2.6) is not included due to the assumption that the Sommerfeld radiation condition is fulfilled a priori. The vanishing variation of (2.6) δI = δc† ∂I ∂c = δc† (◦ R+Kc ) =0 (2.7) yields for the undetermined coefficients c, of the whole assembly of n subdo- mains, c⊤ = [c⊤1 ,c ⊤ 2 ,c ⊤ n ], the system of linear equations, Kc=− ◦ R, where K is Hermitian and positive-definite. The described approach makes use of the following different types of T-systems, defined in local polar coordinates (r,θ) in terms of the Bessel and Hankel functions, Jn,H (1) n , of the first kind 1. T-functions for the bounded subdomain (Herrera and Sabina, 1978) { J0(kr), Jn(kr)cosnθ, Jn(kr)sinnθ } (2.8) 2. T-functions for the unbounded subdomain (Herrera and Sabina, 1978) { H (1) 0 (kr), H (1) n (kr)cosnθ, H (1) n (kr)sinnθ } (2.9) 3. Special purpose functions for the doubly connected subdomain with a circular hole {( J0(kr)− J′0(kb) H′0(kb) H0(kr) ) , ( Jn(kr)− J′n(kb) H′n(kb) Hn(kr) ) cosnθ, (2.10) ( Jn(kr)− J′n(kb) H′n(kb) Hn(kr) ) sinnθ } obtained as the linear combinations of (2.8) and (2.9) fulfilling a priori Neumann homogeneous boundary conditions on the circle of radius b (Stojek, 1998). The T-sets (2.8) and (2.9) approximate the scattered potential in (2.2). They are complete regardless of the specific subdomain considered as long as the condition of being bounded or, alternatively, of being the exterior of a bounded subdomain is fulfilled (Herrera and Sabina, 1978). The T-system (2.10) has been developed directly for the total potential (2.3). Note that the Diffraction loads on multiple vertical circular cylinders... 123 assumption of the completeness of the T-systems is sufficient for our method to be convergent. The above systems are non-singular inside the given domain which yields important numerical advantages, i.e. all boundary integrals can be evaluated bymeans of standard quadrature rules. 3. Wave induced loads on vertical cylinders We consider an inviscid, incompressible fluid in a constant gravitational field. The space coordinates are denoted by xs = (x,y,z) and the gravitatio- nal acceleration g is in the negative z direction. Assuming that the flow is irrotational, the incident complex potential of the progressive linear wave in water of constant depth h0 is known and given by (Mei, 1992) φ i =− igH 2ω coshk(z+h0) coshkh0 eikx (3.1) kx= k1x+k2y= kxcosα+ky sinα where α is the angle of incidence, k stands for the wavenumber vector, ω= √ gk tanhkh0 is a frequency and H is the height of the wave. Let us limit ourselves to an isolated vertical cylinder extending from the seabed and protruding above the free surface. Thus, the incident wave is scat- tered horizontally and the three-dimensional diffraction problem governed by Laplace’s equation may be reduced to two dimensions. Finally, we pose the problem to be solved as follows. Given the incident potential ϕi(x,y)=−igH 2ω eik(xcosα+y sinα) (3.2) find the total potential ϕ(x,y) fulfilling two-dimensional problem (2.1). Once the two-dimensional potential ϕ is found the total velocity potential φ may be calculated from φ(x,y,z) =ϕ(x,y) coshk(z+h0) coshkh0 (3.3) where φ(x,y,z)|z=0 =ϕ(x,y). Then,by integration of thehydrodynamicpres- sure, p = iωρφ, over the wetted body surface we obtain the total horizontal force 124 M.Stojek F= 0 ∫ −h0 F z dz= 0 ∫ −h0 dz ∫ ΓB −pn dΓB (3.4) We investigate the interference effects among neighbouring vertical cylin- ders by means of the force ratio R. We define it as the force amplitude on one member of the system, |Fi|, compared to that predicted for an isolated cylinder, |F0|, under the corresponding wave conditions (the same kb) R= |Fi| |F0| (3.5) 4. Numerical studies In this section, the interference effects amongneighbouringvertical circular cylinders of the same radius b are presented. Since the wave induced loads dependnot only on the diffraction parameter kb, and the angle of incidenceα, but aswell are influenced by the relative geometry of cylinders, the force ratio (3.5) dependence R= f ( kb, L b ,α ) (4.1) is investigated, where L is the distance between cylinder centers. In what fol- lows, the different configurations of cylinders are examined and the results for the force, in the direction of the incident wave, are summarized in comparison with those generated by the program INTACT (the property of Technische Universität Hamburg-Harburg) based on the generalization of an exact alge- braic method for interactions among multiple bodies in water waves (Kage- moto and Yue, 1986). 4.1. Pair of cylinders Let us nowconsider the case of twoneighbouring cylinders (Fig.1) exposed to the action of the incident wave with the following parameters H =1 (height of the incident wave) k∈ (0.02,2) (range of the wavenumber) ∆k=0.02 (increment in the wavenumber) α=0◦ (angle of incidence) (4.2) Diffraction loads on multiple vertical circular cylinders... 125 Fig. 1. Scattering of the incident wave by two vertical, circular cylinders. Problem definition Fig. 2. Subdomain subdivision for two vertical, circular cylinders of Fig.1. Markers stand for the origins of local coordinate systems The radius of each cylinder b is equal to one. The T-fields subdivision, with m=41 T-functions per subdomain, is shown in Fig.2. Here, the use is made of three different T-sets, namely, • (2.8) in 2 bounded subdomains, 126 M.Stojek Fig. 3. Interference between the circular cylinders (Fig.1) shown as the force ratio (4.1) for the first cylinder versus the diffraction parameter kb. Mesh Fig.2, number of T-functions m=41, reference solution-program INTACT Fig. 4. Interference between the circular cylinders (Fig.1). Amplitude of the nondimensional force on the first cylinder versus the diffraction parameter kb. Mesh Fig.2, number of T-functionsm=41, reference solution-program INTACT Diffraction loads on multiple vertical circular cylinders... 127 • (2.9) in the unbounded subdomain, • (2.10) in 2 doubly connected subdomains with a circular hole. The results of numerical studies are summarized in Fig.3, where the force ratio (4.1) is plotted versus the diffraction parameter kb. A comparison of the numerical predictions and the experimental results is presented in Isaacson (1979). Finally, the amplitude of the nondimensional force on the first cylinder for the entire, chosen range of kb is shown in Fig.4. 4.2. Three cylinders in a row Letusnow focus our attention on three identical cylinders in a row (Fig.5). This time, the parameters for the incident wave are assumed as follows H =1 (height of incident wave) k∈ (0.02,2.5) (range of wavenumber) ∆k=0.02 (increment in wavenumber) α=0◦,30◦,60◦ (angle of incidence) (4.3) Fig. 5. Scattering of the incident wave by three vertical, circular cylinders in a row. Problem definition The radius of each cylinder b is again equal to one, but the distance be- tween two neighbouring cylinder centers is increased to L=4b. The subdivi- sion into one unbounded and 15 bounded T-fields is presented in Fig.6. Note that the use of the three doubly connected subdomains, each equipped with special purpose functions (2.10) for a circular hole, is essential for this sub- domain mesh. The number m of T-functions over one subdomain is taken to be m=41. 128 M.Stojek Fig. 6. Subdomain subdivision for three vertical, circular cylinders in a row (Fig.5). Markers stand for the origins of local coordinate systems Thenumerical results for the interference effects are presented for the three cylinders in the following way: • force ratio (4.1) versus the diffraction parameter kb for α=0◦ in Fig.7, • amplitude of the nondimensional force on the first cylinder versus the diffraction parameter kb for different α (Fig.8). One can observe that the first cylinder in a row is themost affected by the interference of the waves. Nevertheless, in the case of α = 0◦, the maximal force ratio R for the cylinder in the middle is very close in value to that for the first cylinder (Fig.7). 4.3. Four cylinders In this section, the interference effects for four cylinders (Fig.9) are in- vestigated. The respective T-fields subdivision, consisting of one unbounded subdomain with T-set (2.9) and four doubly connected subdomains provided with special purpose functions (2.10) for a circular hole, is given in Fig.10. The number m of T-functions over each subdomain is 41. The height of the Diffraction loads on multiple vertical circular cylinders... 129 Fig. 7. Interference among three circular cylinders in a row (Fig.5). Mesh Fig.6, number of T-functions m=41,markers indicate the solutions obtained by program INTACT Fig. 8. Interference among three circular cylinders in a row (Fig.5). Amplitude of the nondimensional force on the first cylinder versus the diffraction parameter kb. Mesh Fig.6, number of T-functions m=41,markers indicate the solutions obtained by program INTACT 130 M.Stojek Fig. 9. Scattering of the incident wave by four vertical, circular cylinders. Problem definition Fig. 10. Subdomain subdivision for four vertical, circular cylinders of Fig.9. Markers stand for the origins of local coordinate systems Diffraction loads on multiple vertical circular cylinders... 131 incidentwave H is equal to one, and the range of thewavenumber is increased to 5, namely k∈ (0.02,5). The increment in the wavenumber for the most of the range is taken to be ∆k = 0.02, only when a denser approximation is required, is the step reduced to 0.001. First, the incident wave propagating in the x Cartesian coordinate di- rection (α = 0◦) is examined. The results for force ratio (4.1) versus the diffraction parameter kb are displayed in Fig.11. Next, the angle of incidence α = 45◦ is considered and the force ratio R versus the diffraction parame- ter kb is plotted in Fig.12. Here, one can observe that for the first and third cylinder the force ratio R is much greater in the vicinity of kb = 2.76 than away from the ”resonance”. The details for the range k∈ (2.5,3) are given in Fig.13. Finally, the amplitude of the nondimensional force on the first cylinder versus the diffraction parameter kb for α=45◦ is shown in Fig.14. 5. Concluding remarks Theapplication of the so-called least-squares T-elements to the calculation of the diffraction loads on the multiple vertical cylinders with circular cross sections in waves has been presented. Such an approach has led to sufficiently accurate solutions for only few subdomain divisions due to the fact that the oscillatory nature of thewave is inherent in all kinds of Bessel’s functions, and as a consequence, in contrast to a polynomial approximation, the dimension of the T-subdomain need not be related in any way to the wavelength. We also recall that the T-type approach has the benefit of a priori satisfaction of the radiation condition at infinity, though, special methods developed for infinite domains are not needed. In addition, the use of the special purpose functions for the doubly con- nected subdomain has resulted in the quick convergence and good accuracy for amplitudes and phases of the calculated forces on the wetted surfaces. It should also be recalled that for applications in offshore engineering, the most important issue is the accuracy of the results on thewetted surface. Thus, one can suppose, that the best approximate solution should fulfill optimally the relevant boundary conditions. Indeed, the presented approach manifests such a property. The computational results showed the influence of the neighbouringbodies on the diffraction forces. The force amplification factor obtained for the array of cylinders and α=45◦ in the ”resonance range”was almost 5which showed the importance of the effects of wave trapping in the interior of the array. 132 M.Stojek Fig. 11. Interference among four circular cylinders (Fig.9). Mesh Fig.10, number of T-functions m=41, markers indicate the solutions obtained by program INTACT Fig. 12. Interference among four circular cylinders (Fig.9). Mesh Fig.10, number of T-functions m=41 Diffraction loads on multiple vertical circular cylinders... 133 Fig. 13. The ”resonance” range for the interference among four circular cylinders (Fig.9). Mesh Fig.10, number of T-functions m=41, markers indicate solutions obtained by program INTACT Fig. 14. Interference among four circular cylinders (Fig.9). Amplitude of the nondimensional force on the first cylinder versus the diffraction parameter kb. Mesh Fig.10, number of T-functions m=41, reference solution-program INTACT 134 M.Stojek References 1. Babuška I., Sauter S.A., 1997, Is the pollution effect of the FEMavoidable for the Helmholtz equation considering high wave numbers? SIAM J. Numer. Anal., 34, 2392-2423 2. Bettes P., 1992, Infinite Elements, PenshawPress, Sunderland 3. ColtonD., Kress R., 1992, Integral EquationMethods in Scattering Theory, Wiley, NewYork 4. Farhat C., Harari I., Franca LP., 2001, The discontinuous enrichment method,Comput. Methods Appl. Mech. Engng., 190, 6455-6479 5. Gerdes K., 1998, A summary of infinite element formulations for exterior Helmholtz problems,Comput. Methods Appl. Mech. Engng., 164, 95-105 6. Givoli D., Patlashenko I., Keller J.B., 1997, High-order boundary con- ditions and finite elements for infinite domains,Comput. Methods Appl. Mech. Engng., 143, 13-39 7. Harari I., Barbone P.E., Slavutin M., Shalom R., 1998, Boundary in- finite elements for the Helmholtz equation in exterior domains, Int. J. Numer. Meth. Engng, 41, 1105-1131 8. Harari I., Hughes T.J.R., 1992, A cost comparison of boundary element and finite elementmethods for problems of time-harmonic structural acoustics, Comput. Methods Appl. Mech. Engng., 97, 77-102 9. Herrera I., Sabina F.J., 1978, Connectivity as an alternative to boundary integral equations:Construction of bases,Proc. Nat. Acad. Sci. USA,75, 2059- 2063 10. Ihlenburg F., 1998,Finite Element Analysis of Acoustic Scattering, Springer 11. Isaacson M. de St. Q., 1979, Interference effects between large cylinders in waves, J. Petroleum Tech., 31, 4, 505-512 12. Kagemoto H., Yue D.K.P., 1986, Interaction among multiple three- dimensional bodies in water waves: an exact algebraicmethod, J. Fluid Mech., 166, 189-209 13. Mei C.C., 1992,The Applied Dynamics of Ocean SurfaceWaves,World Scien- tific, Singapore 14. Stojek M., 1998, Least-squares Trefftz-type elements for the Helmholtz equ- ation, Int. J. Numer. Meth. Engng, 41, 831-849 Diffraction loads on multiple vertical circular cylinders... 135 Obciążenia dyfrakcyjne fal wodnych na układy pionowych cylindrów o przekroju kołowym za pomocą elementów skończonych ”typu Trefftza” Streszczenie Artykuł prezentuje zastosowanie elementów skończonych ”typuTrefftza” do obli- czeń obciążeń dyfrakcyjnych fal wodnych na układy pionowych cylindrówo przekroju kołowym.Metoda jest oparta na aproksymacji rozwiązań w podobszarach za pomo- cą skończonej sumy szeregu funkcji ”T-zupełnych” (funkcji kształtu). Rozwiązanie globalne jest otrzymane przez ”zszycie” rozwiązań lokalnych i wymuszenie warun- ków brzegowych za pomocą metody najmniejszych kwadratów. Jednorodne warunki brzegowe Neumanna i warunek radiacyjny Sommerfelda są zawarte w specjalnych elementach skończonych. Manuscript received May 6, 2002; accepted for print October 8, 2002