Fast Optimization of Tuned Mass-Spring Arrangement for Reduction of Acoustic Power Radiation of Submerged Ribbed Cylinder

In this paper, acoustic power radiation of a submerged finite length ribbed cylinder subject to a harmonic point load is minimized by a new fast scheme. For this purpose, two arrangements of non-uniformly distributed sequential point masses and mass springs attached on stiffening ribs of the cylinder are used to optimally reduce the acoustic power radiation. A fully coupled analysis is here carried out based on Finite/Boundary element (FEM/BEM) model. Instead of direct BEM formulation, two beneficial procedures have been proposed for computing BEM matrices in each frequency line. In order to fast solving of equations, the Krylov vectors (produced via Ritz or Arnoldi iterative procedures) and structural mode shapes (modal truncation approach) have been used and validated before performing optimization. As a result, the best strategy for evaluation of response and cost function is using Taylor series expansion for computing BEM matrices and applying Krylov vectors for order reduction. The results show good agreement with previous studies and experiments. The optimization results show noticeable reductions in the acoustic power radiation. In point mass optimization, the most of additional masses has been placed in regions which are near to the excitation point whereas for the absorber design, they are put in the places in opposite side of the excitation point.


INTRODUCTION
Nowadays, optimization of acoustic power radiated from submerged structures is important in many applications.Analytical modeling methods are applicable to few actual systems so, numerical methods is an essential part of simulations especially in optimization processes.In low frequency range, three numerical methods dominate for acoustic part simulation: Finite element Method (FEM/IEM), Boundary Element Method (BEM) and Rayleigh integral method.
In vibro-acoustic optimization context, initial studies were accomplished by Olhoff (1974) in 1974.In his work, the thickness function of a simply supported rectangular plate was determined such that the fundamental natural frequency of transverse vibrations attains an optimum value.
Knowing the potentials of elastic and inertial modifications on structural noise reduction, several researches have been made for vibro-acoustic simulation and optimization of structures by use of mass/springs.Nagaya and Li (1997) have optimized structural noise radiation from a rectangular plate by appropriate placement of dynamic absorbers.Constans et al. (1998) have shown that significant sound power reduction can be ensued by optimal placement and sizing of small point masses in the semi-cylindrical shell structure using Simulated Annealing Algorithm.They deduced that the optimal small point masses can alter the critical mode shapes to quieter modes of vibration and causes to sound power reduction.Similar work had been also done by Cheng and Wang for a fluid-loaded beam (Cheng & Wang (1998)).
Ratle and Berry (1998), used a Genetic Algorithm for minimizing emitted sound pressure level of vibrating plate carrying point masses.It has been shown that the use of Genetic Algorithm can help finding original solutions that would not have been found by any intuitive means; because no a priori bracketing of the optimal solution is needed.Also they pointed that choice of cost function and frequency band are important for fast and effective optimization.Qiqiang and Zhijian (2012) have optimized a stiffened cylindrical shell with FEM and BEM in 2012.They minimized the structural response with employment of power and modal radiation efficiency as objective functions.Point mass placement on rings, has been led to a predictable conclusion that embedding point mass on regions close to excitation is more efficient for structural response reduction.
In 2016, Kumar et al. (2016) studied the effect of attached discrete patches/point masses on sound radiation of a plate by means of FEM for structural part and Rayleigh integral for acoustic part.They have investigated the effect of position of discrete mass and thickness distribution of discrete patches.They deduced that both of them have significant effect on sound radiation of plate but under the constraint of fixed additional mass, point masses can alter natural frequencies more effectively by adjusting their positions.Michielsen et al. (2016) formulated a linear quadratic regulator based optimization problem in order to minimize the broad-band low-frequency domain vibration and acoustic response of a baffled simply supported plate by means of multiple optimally tuned mass-spring-damper (TMD) system.Their results indicate that TMD have great potential to reduce the broadband low frequency response of vibro-acoustic systems.Also it can be concluded that there are fundamental differences between the optimal TMD if one minimizes the kinetic energy or the far-field radiated sound power.
A survey of methods, applications and various features of structural acoustic optimization for passive noise control can be found in a review paper (Ranjbar et al. (2010)).Various gradient based optimization methods have been used in vibroacoustic context especially in narrow frequency bands and single mode optimizations.Sequential Linear Programming (SLP), the method of moving asymptotes (Tinnsten et al. (2002)), the method of feasible directions, Sequential Quadratic Programming which is found to be robust and suitable for finding local minima in several vibroacoustic cases, Zheng et al. (2006), level-set (Isakari et al. (2017)), Pattern (Direct) search (Lang et al. (1975)) and other gradient-based methods were successfully applied in the field of vibroacoustic.Most of gradient based methods are categorized as local search methods because they converge to the nearest local minima.Among stochastic search methods, Genetic Algorithm (GA), Particle Swarm Optimization (PSO) and Simulated Annealing (SA) have the most applications in vibro-acoustic optimization (Marburg, (2002)).
Another important issue is high computational cost of vibro-acoustic evaluation and optimization of coupled systems especially for large scale problems.That is why a few works have been reported in optimization of coupled systems.Christensen and Olhoff (1998) considered a desired directional pattern for noise emitting loudspeaker diaphragm and minimized summation of square of deviations (errors) in arbitrary positions by changing the diaphragm structural parameters.Minimization has been carried out in three distinct frequencies.Very thin structure imposes the consideration of two-way coupling in their formulation.A multi-criteria optimization approach has been used by Akl et al. (2002) to find optimal design of underwater shell structures.They minimized shell vibration, sound radiation, weight of the stiffening rings and the cost of the stiffened shell simultaneously.A structural-acoustic optimization approach has been presented by Shepherd and Hambric (2014) for minimizing the radiated power of structures with heavy fluid loading excited by complex forcing functions.The procedure has been demonstrated on a curved underwater panel excited by a point drive and by turbulent boundary layer flow.The objective function was a weighted sum of total sound power and panel mass.
From mathematical viewpoint, the expansion of structural response in terms of eigenvectors is the best way to compute the structural response.Computational efficiency can be increased if some assumptions apply.For frequencies less than the first eigenfrequency of the system and for light fluid, structural response can be calculated by expansion of static Ritz vectors (2002).Puri (2011) introduced dimension reduction techniques for fully coupled, interior structural-acoustic systems based on Krylov Subspaces.For the test cases investigated in his research, it is shown that using the reduced order modeling technique causes very significant reduction in simulation time, while maintaining the desired accuracy of the state variables, i.e. displacements and pressures (2011).
In contrast with internal acoustics and FEM/FEM models, a little work has been carried out for efficiently solving vibro-acoustic problems that involve exterior unbounded domains via FEM/BEM models.The most challenging Latin American Journal of Solids and Structures, 2019, 16(1), e144 3/16 problems are non-symmetric (and hence non-orthogonal eigenvectors), frequency dependent and huge matrices to be solved.Eigenvalue extraction at each frequency which is a time consuming procedure, makes it impossible in the case of frequency dependent mass matrix.
To the authors' best knowledge, the problem of "fully coupled vibroacoustic (e.g.acoustic power radiation) optimization of finite length submerged ribbed cylinders" has not been addressed in the literature yet.A submerged ribbed cylinder has been considered as a numerical test case (the model described by Zhou & Joseph (2005)) for coupled vibroacoustic FEM/BEM modeling.At first, two procedures have been proposed for computing BEM matrices in each frequency line.In order to solve governing equations, dimension reduction is inevitable, so the Krylov vectors (produced via Ritz or Arnoldi iterative procedures) and structural mode shapes (modal truncation approach) have been used and compared before performing optimization.The results showed good agreement with previous studies and experiments.Although Krylov reduction concept has been developed for simulating interior acoustic problems (cabin noise, etc.), using these basis vectors along optimization process of exterior vibroacoustic problems has not been reported in the literature.Finally, the optimum concentrated (point) mass and the dynamic absorber (TMD) arrangements were considered for the reduction of average radiated power using Genetic Algorithm.As expected, the majority of total point masses have been proposed by GA for the nearest regions to the excitation point but, this is not true for the case of dynamic absorbers (TMD).A noticeable point in the design of silent submerged vessels is that the most changes in this case proposed on the opposite side of the excitation point.

Coupled vibroacoustic equations
Dynamics of cylinder structure is modeled via FEM with reasonable accuracy due to diversity of complex geometry and boundary conditions.Dynamic equation of motion for a submerged structure in FEM has the following form in harmonic excitation (Rao, 2004): Where M st and K st are assembled mass and stiffness (of cylinder wall, caps and its stiffeners) matrices of whole structure respectively, η is structural (Hysteretic) damping coefficient, f ext is exciting vector including force and moment contributions and X s is global coordinate vector for whole structure including translational and rotational global DOF's in each node.P is elemental/nodal acoustic (dynamic) pressure vector of wetted area nodes and A is distribution matrix.A ik is force/moment excitation correspond to i th DOF due to unit pressure at k th element.Although FEM results are used for comparison, BEM has been used in this study for fluid-structure interaction simulation i.e. computing radiation impedance matrix (Z rad ) (Rao, 2004): Matrix L s relates nodal displacements of nodes to normal displacement of element centers by: Where X n is the element normal displacement vector.Displacement vector at each element center is [17]: w e is deflection of element "e" center in local coordinate system, d is a vector including deflection ∂w/∂x and ∂w/∂y (in local coordinate system) for four nodes of element "e" and also N is corresponding shape function in row matrix form Figure 1 shows element "e" and its parameters.

BEM in discretized form
Figure 2 shows a schematic view of exterior vibro-acoustic problem with various acoustic boundary conditions.
Where g is free space Green's function, dS r is the area of radiating element "r", n ̂ is local normal direction and N is number of shell elements on the interface area.Coefficient c has magnitude of 0, 1 or 0.5 for out of field, domain and interface points respectively.This procedure leads to a matrix equation between acoustic pressure and normal velocity vectors: Because source point is always on the boundary, c is selected to be 0.5 and H=0.5I.Compare to Eq. ( 2), the radiation impedance matrix for shell is: Substitution of Eq. ( 2) into Eq.( 1) results: Latin American Journal of Solids and Structures, 2019, 16(1), e144 5/16 Clearly, the total fluid effect emerges in third term which can be written as 2ρ f AG(ω)L s is known as frequency dependent "Added mass matrix".

Model updating due to adding point masses
If node r has additional mass equal to m r , only three elements of M st must be modified.With reference to FEM standard procedures, additional mass matrix can be constructed by kinetic energy expression:

Model updating due to dynamic absorbers
As shown in Figure 3, additional masses are connected to original system via S 1 S 2 . .S na springs.

GOAL FUNCTION
As mentioned in the introduction, although several goal functions have been proposed for optimization process, majority of researchers have chosen the average radiated power in specified frequency range for external vibroacoustic systems so it is desirable to define the goal function as: Where P(ω) is radiated power at a specific frequency: Latin American Journal of Solids and Structures, 2019, 16(1), e144 6/16 ℜ( ) and ( ) * denote the real part and conjugate of a complex number and the integral is evaluated over the entire radiation surface.Very wide range of frequency may result in a little enhancement in the goal function so it is necessary to find out the excitation frequency and the frequency range.

Model order reduction
Computing structural response from Eq. ( 8) is time consuming thus, use of reduction scheme is necessary because each evaluation of fitness (goal) function requires several calculations at some discrete frequencies.Therefore, two reduction strategies are used to model order reduction.

A. Dry mode shapes modal reduction
Normalized truncated dry mode shape matrix (  ) can be used as the first choice for projection: By pre-multiplying    into Eq.( 8), a new equation in reduced dimension space is formed:

B. Krylov reduction
In Krylov subspace based reduction method,   might be replaced by   (Krylov subspace-based projection matrix); so by pre-multiplying    into Eq.( 8): This procedure has shown its applicability in interior vibroacoustics (Puri & Morrey, 2011), however the method has not been addressed for exterior problems.In this case, the problem is different because system matrices are frequency dependent.A successful application of Krylov vectors for dimension reduction in exterior vibroacoustic analysis can be found in (Dadkhah et al. 2015).A ribbed simply supported rigidly baffled plate with flat side in contact with water (Figure 4) has been considered as a case study for dimension reduction using Krylov vectors by 5 and 25 vectors.Point FRF of cross point "A" is shown in Figure 4 (b) both by full analysis and dimension reduction using Krylov vectors.As can be seen, this technique has desirable accuracy at least in first four modes.This accuracy can be enhanced by selecting more Krylov vectors.Also it has been found that the shift frequency factor does not have significant effect on the accuracy of vibroacoustic analysis of submerged plate.More details can be found on (Dadkhah et al. 2015).In this work, Krylov subspace based Ritz vectors have been used for model order reduction of submerged ribbed cylinder.

Interpolation of radiation impedance matrix
Eq. ( 7) can be evaluated faster by employing an interpolation technique for radiation impedance matrix (  ) calculation.In fact, the components of this matrix vary smoothly with frequency, so it is possible to use frequency interpolation to decrease the computational efforts.In this work, the radiation impedance is calculated at so-called "master frequencies" and then interpolated at each frequency linearly (or by higher order polynomials).This technique is useful to reduce the computational time without major decrease in accuracy.

Interpolation of Modal (Reduced) Added Mass Matrix
Regarding to Eq. ( 15), the modal added mass matrix (3rd term in the LHS) is also frequency dependent, therefore, one can interpolate      ()    matrix instead of   , so the computation time can be saved further.
Obviously, this can be done only when reduction is applicable.The only drawback is that in evaluating goal function e.g.radiated power, the pressure quantity is needed which enforces the calculation of   in each frequency.This problem can be solved by the technique introduced in section 4.4.

Taylor Series Expansion of G
In BEM procedure, calculation of each element of G involves the integration of Green's function: Where   is Green's function between source (j element center) and receiver (point i) and   is the surface of element j.The term  − in Green's function leads to frequency dependency of G, therefore, all elements must be calculated at each frequency or saved before frequency loop to use during the simulation.Using Taylor series expansion in this study, time and cost of computation is saved further: () 4   24 − … = ∑ Latin American Journal of Solids and Structures, 2019, 16(1), e144 8/16 For small "kr" argument, only few terms are sufficient for reasonable accuracy.For example a good approximation up to kr=2ð can be achieved by only the first eight terms of the series.So: Where Hence, G can be rewritten as weighted summation of some frequency independent matrices which can be evaluated before simulation: In above relation, only scalars (  and   ) are frequency dependent.

Model description
A ribbed cylindrical shell with two end caps is considered as a numerical case study and verification according to Figure 5. Table 1 shows the structural and fluid properties.The excitation force is   () = 1 ⋅ cos()  which is inserted on point A in x direction (according to Figure 5. Fluid domain has been modeled separately by FEM acoustic elements and BEM.Total DOF of FEM/FEM model is 18699. Latin American Journal of Solids and Structures, 2019, 16(1), e144 9/16

FEM/FEM model
Calculated natural frequencies of submerged ribbed cylinder are shown in Figure 6.The maximum error between model used FEM/FEM model with those reported in (Zhou & Joseph, 2005) is 3.8% .The first six mode shapes of the submerged cylinder agree with (Zhou & Joseph, 2005).First two of mode shapes are shown in Figure 7.The results indicate that four terms are sufficient for acoustic impedance evaluation in 0 to 300 Hz frequency range in Taylor series expansion.Number of basis vectors i.e. number of columns of  is set to be 50.As indicated in Figure 8, FRF of point "A" shows a small error between full order FEM model and reduced BEM model by Ritz vectors.

Figure 8: Point FRF of A, FEM/BEM by Ritz vectors reduction and 106.4 Hz shift (___), FEM/BEM by dry modes reduction (____); Full FEM/FEM (----)
Latin American Journal of Solids and Structures, 2019, 16(1), e144 10/16 The pattern of acoustic pressure radiated by vibrating cylinder has been compared with experimental results (Zhou & Joseph, 2005) and plotted in Figure 9.This figure shows resonable accuracy in acoustic pressure calculation.The efficiency (precision and computational cost) of proposed methods has been studied by several tests before running optimization procedure.These comparisons have been done by a Core i7-2600, 3.4GHz CPU with 16 MB of RAM PC.

a)
Interpolation: After calculation of () in master frequencies, this matrix is calculated via interpolator polynomial (n=3) in other frequencies.Results show that much less computing time can be achieved by reduced added mass interpolation instead of () as indicated in Table 2 .

b)
Taylor Series Expansion: As can be seen from Table 2, using Taylor series expansion has been led to significant reduction in every evaluation of ().

c)
Full vs. Reduced System of Equations: Table 3 shows a brief comparison between calculation time in full and reduced (via dry modes and Ritz vectors) system of equations.As can be seen, significant reduction in calculation time can be achieved by Ritz vectors with reasonable accuracy.
As a result, the best strategy for response evaluation is applying Ritz vectors for order reduction and using Taylor series expansion for computing ().

Optimization Algorithm
Heuristic methods require a large number of function evaluations but they have outperformed classical optimization methods such as linear or even nonlinear programming in several cases especially when considering multimodal functions thus requiring an optimization algorithm able to sort the globally optimal solution out of a large number of locally optimal solutions (Ratle and Berry (1998)).Based on statistical methods, crossover and mutation operators play an anti-trapping role in a Genetic Algorithm (GA) and they work better than other similar optimization procedures.Therefore, GA was chosen as optimization algorithm in this study.Converging best cost function to a fixed value and also converging average of cost function in every generation to the best cost function is a measure of approaching global minimum of function.A number of successful GA applications have been reported in the literature

Optimal point mass placement
In this section, optimal point mass placement is obtained in such a way that the average radiated power,   (with  1 and  2 equal to 100 and 160Hz respectively) is minimized.These masses are placed only on rings and the harmonic force is exerted on point 30 (according to Figure 10).Based on design considerations, it is assumed that no more than 4% of total mass (40 kg) is permitted to add as point masses.Modal analysis of optimized system shows slight natural frequency and mode shape changes (maximum of -5.3Hz for 5 th mode) due to tight constraint of additional mass.Clearly, a little mass is proposed for point 30 (excited point) but three middle rings (which have the most amplitude in first two modes) have 66% of total point masses as highlighted in Table 4. On the other hand, 35% of total mass proposed to be placed on first column which includes the excitation point.However, adjusting mass distribution according to first or second mode amplitudes of points led to very weak effect in power reduction.

Optimal dynamic absorber placement
Similar to point mass in this case several scenarios have been set to obtain the best arrangment.Population size and number of generations that are the main parameters for GA have set to be 150 and 80 respectively and adaptive feasible option has been selected for mutation function (several mutation functions has been tested).Again according to Figure 12, converging average of fitness function in last generations toward the best fitness function is a good sign of global optimization.Table 5 shows the optimization results for placement of absorber on all points depicted in Figure 10.Again final value of cost function is almost fix in different initial populations.Most of mass changes proposed for the left column which is located on opposite side of excitation.Three mid rings have 48% of total added mass.The results show that use of mass-spring system (plus inherent structural damping of 1%) instead of point mass can decrease the cost function more significantly (20.7 dB reduction).
Modal analysis of optimized system shows slight natural frequency and mode shape changes of submerged ribbed cylinder (except second mode) similar to point mass optimization (maximum of +5.3Hz for 6 th mode) due to tight constraint of additional mass.
The radiated power in initial design, with optimized point mass arrangement and with optimized absorber arrangement is shown on Figure 13.As can be seen, optimized point mass arrangement does not have substantial effect on overall radiated power whereas optimized dynamic absorber (TMD) arrangement has reduced radiated Latin American Journal of Solids and Structures, 2019, 16(1), e144 14/16 power at 110, 130, 140, 150 and 160 Hz (which selected to evaluate Eq. ( 15)).Also there are several local peaks due to discrete absorber resonance peaks which are controlled by inherent damping mechanisms.The first natural frequency of structure has been decreased in both optimized situations.In addition, the second mode in Figure 8 (f=142.4Hz) has been vanished by properly tuning of absorber parameters and this is the main reason for obtaining good reduction.

CONCLUSIONS
In this study, vibro-acoustic optimization of a submerged shell structure was investigated.As shown by several studies, BEM imposes time and memory limitations due to frequency dependency of fully populated BEM matrices.Therefore, two procedures (interpolation techniques and a sufficiently precise Taylor series expansion) have been proposed for computing the BEM matrices in each frequency line.Table .2shows a drastic effect on calculation time especially for second procedure.Krylov vectors and structural mode shapes (modal truncation approach) have been used and compared before doing optimization.The orthogonality of   with respect to   and   is not sufficient for accurate results in dry mode reduction because fluid effects leads to non-diagonal modal mass matrix.As a result, the best strategy for evaluation of response and cost function is using Taylor series expansion for computing BEM matrices and applying Krylov vectors for order reduction.
In the case of point mass optimization, average radiated power reduction of 11 dB (which means 3.55 times smaller acoustic radiated power) has been achieved.Natural frequencies and mode shapes change slightly and three middle rings (which have the most amplitude in first two modes) have 66% of total point mass.35% of total mass proposed to be placed on first column which includes the excitation point.However, adjusting mass distribution according to first or second mode is not so beneficial.
In selecting tuned mass damper parameters the best position for minimizing sound power or structure kinetic energy in the case of single absorber is the point of excitation (as pointed by Michielsen & Arteaga,2016) but for multiple TMD's most of additional mass proposed for the left column which is located on opposite side of excitation (Contrary to point mass case).Three mid rings have 48% of total added mass.The second mode has no longer considerable participation on radiated power by properly tuning of TMD parameters and this is the main reason for achieving good reduction.20.7 dB reduction (which means 10.84 times smaller acoustic radiated power) is the best situation in dynamic absorber (TMD) attachment case.

Figure 1 :
Figure 1: Element "e" of shell and its local coordinate

Figure 4 a
Figure 4 a): A ribbed simply supported rigidly baffled plate in semi-infinite fluid domain.b): Point FRF of point A (Vertical excitation)

Figure 6 :
Figure 6: First six Natural frequencies of submerged ribbed cylinder

Figure 10 :
Figure 10: Additional point mass/absorber locations After 80 generations of 150 individuals for GA and adaptive feasible option for mutation function (several mutation functions has been tested), the best result is 11 dB reduction in cost function.The amount of reduction in cost function is relatively desirable and is almost fix value in different initial populations.As can be seen from Figure 11 converging average of fitness function in last generations toward the best fitness function is a good sign of global optimization.Optimal point mass values have been shown in Table4.

Figure 11 :
Figure 11: History of fitness function vs generations in point mass optimization

Figure 12 :
Figure 12: History of fitness function vs generations in dynamic absorber optimization

Table 1 :
Model properties

Table 2 :
Different techniques for calculation of fluid effect and their efficiency Latin American Journal of Solids and Structures, 2019, 16(1), e14411/16

Table 3 :
Calculation time in different reduction methods

Table 4 :
Details of point mass optimization, bounds, constraints and results