Augmentation of DAA staggered – solution equations in underwater shock problems for singular structural mass matrices

The numerical solution of underwater shock fluid – structure interaction problems using boundary element/finite element techniques became tractable through the development of the family of Doubly Asymptotic Approximations (DAA). Practical implementation of the method has relied on the so-called ugmentationof the DAA equations. The fluid and structural systems are respectively coupled by the structural acceleration vector in the surface normal direction on the right hand side of the DAA equations, and the total pressure applied to the structural equations on its right hand side. By formally solving for the acceleration vector from the structural system and substituting it into its place in the DAA equations, the augmentation introduces a term involving the inverse of the structural mass matrix. However there exist at least two important classes of problems in which the structural mass matrix is singular. This paper develops a method to carry out the augmentation for such problems using a generalized inverse technique.


Introduction
The numerical solution of transient coupled field problems has been an important topic for many years.In particular, the fluid -structure interaction associated with underwater shock of submerged or semi -submerged structures has been the subject of many investigations and has led to the development of a number of methods of solution.
Numerical stability has been the most important topic of course, otherwise any answers are meaningless.Next in importance is accuracy, and both of these topics were addressed in [1] that also includes a number of pertinent references on the subject.Although that work was devoted to the coupling of a finite element model of a structure to a boundary element model of a fluid described by the Doubly Asymptotic Approximation (DAA) [2][3][4], the ideas presented there are pertinent to the numerical solution of coupled field problems in many disciplines.
One important conclusion drawn in [1] is that a so -called staggered solution procedure can be used in which distinct software modules alternately advance the solution in time for the fluid and structural systems and are interfaced through extrapolation of the coupling terms.However, this method is severely limited to small time step sizes with the equation systems in their original form.Fortunately, another conclusion made there is that a process called augmentation can be used wherein information from one system is introduced into the other system by simple substitution from the governing equations.The method found in [1] that is most suitable is called the displacement extrapolation technique, and in conjunction with appropriate linear multistep operators, achieves unconditional stability with regard to time step size for linear systems.Application of the method to nonlinear structural response has also demonstrated such stability.
A primary requirement of the method is that the structural mass matrix be nonsingular as well as remaining constant during the transient response analysis.However, there exist at least two classes of problems in which the mass matrix is singular.One is the case in which the finite element model includes offset beam elements coupling some, but not all, translations and rotations.Rotational degrees of freedom untouched by the beams will still have zeroes on the diagonal in a lumped mass approach.This also happens when substructuring/super -elements are used in part of the model.Rotational degrees of freedom untouched by these constructs will similarly possess zeroes on the diagonal.In these cases the structural mass matrix used in the augmented DAA equations is typically approximated by lumping mass terms to the diagonal, perhaps with some scaling to preserve wet surface energy.
The purpose of this paper is to extend the concept of augmentation for the surface interaction Doubly Asymptotic Approximations to solve such problems using a generalized inverse technique that can be shown to lead to solutions that exhibit least squares error [5].The Underwater Shock Analysis (USA) code [6] has been modified to incorporate this new technique to evaluate its usefulness and accuracy.Some example problems are first investigated in which the structural mass is, in fact, nonsingular to demonstrate the accuracy of the method.Then some of these cases are rerun in which an entry in the mass matrix is set to zero to demonstrate the robustness of the new augmentation.The results clearly show the remarkable accuracy and utility of this method.

Structural response equation
The discretized differential equation system for the dynamic response of a structure can be expressed in the general form where x is the structural displacement vector, M s , C s and K s are the symmetric structural mass, damping and stiffness operators, respectively.In the event that the system is linear these become constant matrices throughout the structural response.f is the external force vector, and a dot denotes a temporal derivative.
For excitation of a submerged structure by an acoustic wave, f is given by where p I and p S are nodal pressure vectors for the wet -surface fluid mesh pertaining to the (known) incident wave and the (unknown) scattered wave, respectively, A f is the diagonal area matrix associated with elements in the fluid mesh, and G is the transformation matrix that relates the structural and fluid nodal surface forces.More will be said about G in the next subsection.

DAA 1 equation
The first order Doubly Asymptotic Approximation may be written in terms of quantities evaluated only on the wet surface of the structure as where u S is the vector of scattered -wave fluid particle velocities normal to the structure's wet surface, ρ and c are the density and sound velocity of the fluid, respectively, and M f is the symmetric fluid mass matrix for the wet -surface fluid mesh.This matrix is produced by a boundary -element treatment of Laplace's equation for the irrotational flow generated in an infinite, inviscid, incompressible fluid by motions of the structure's wet surface; it is fully populated with non -zero matrix elements.When transformed into structural coordinates, the fluid mass matrix yields the added mass matrix, which, when combined with the structural mass matrix, yields the virtual mass matrix for motions of a structure submerged in an incompressible fluid [7].
For excitation by an incident acoustic wave, u S is related to structural response by the kinematic compatibility relation where the superscript "T " denotes matrix transposition.Equation ( 4) expresses the constraint that normal fluid particle velocity match normal structural velocity on the wet surface of the structure.The fact that the transformation matrix relating those velocities is G T follows from the invariance of virtual work with respect to either of the wet surface coordinate systems.Generally, G is a rectangular matrix whose height greatly exceeds its width, inasmuch as the number of structural degrees of freedom usually considerably exceeds the number of those in the fluid system.The more convenient computational form of the DAA 1 is obtained by substituting Eq. ( 4) into Eq.( 3) and then premultiplying the result by

Interaction equations and augmentation
The introduction of Eq. ( 2) into Eq.( 1) with Eq. ( 5) yields the interaction equations where D f 1 is the symmetric matrix multiplying p S in Eq. ( 5) and is given by As discussed in [1], numerically stable and accurate computation is possible through the application of a staggered solution procedure that may be effected as follows.The first of Eq. ( 6) is simply rearranged to give It is at this point that the new augmentation procedure deviates from the standard one.For comparisons between them, the standard method is first developed based upon the assumption that M s is nonsingular.Then ẍ can be determined from Eq. ( 8) and the term involving it in the second of Eq. ( 6) becomes By defining another symmetric matrix D s as then the second of Eq. ( 6) becomes Note that an additional term now appears on the left side of Eq. (11), D s p S .This process of injecting one of the coupled equations into the other is termed augmentation, hence the first of Eqs ( 6) and ( 11) are herein termed the augmented interaction equations.This form of the DAA 1 has several desirable computational features.As mentioned above, it results in a numerically stable computational system in combination with the structural dynamics equations.Physically this is related to the fact that Eq. (11) incorporates the complete virtual mass of the fluid and structural models, albeit in the sum of their respective inverses.In essence, a damping term has been added to the DAA 1 that favorably affects the roots of its characteristic equation and counteracts the delayed feedback of the radiation damping energy from the fluid system into the structural system.In addition, Eq. ( 11) is also self starting in the sense that only the known incident wave information is required at t = 0 since both ẋ and x are zero for a structure being excited from rest.
Although u I can be discontinuous at t = 0 so that uI is a delta function, a further modification to Eq. ( 11) is made to avoid this problem.The method is described in [6] and is not of interest in this paper.
Of course, the attributes of Eq. (11) can only be enjoyed assuming that the structural mass matrix is nonsingular.When such is not the case for some classes of problems, then the solution of Eq. ( 8) must be effected by other means.Here, use is made of generalized inversion techniques.

Generalized inversion
To motivate the development, it is noted that the structural coupling term G T ẍ in the DAA 1 equation in Eq. ( 5) is just the structural acceleration vector normal to the wet -surface of the structure, i.e., in the fluid coordinate system as To determine ẍ from this expression, the inverse of G T must be found.Inasmuch as G itself is rectangular, no inverse is available in the usual sense.However, the so -called left inverse may be used here that is defined as [5] Applying Eq. (13) to Eq. ( 12) yields where or This can now be substituted into Eq.( 8) to give By premultiplying this expression by G −1 L and using its definition from Eq. ( 13) there is obtained By noting that the combination of the various appearances of G in the first term on the right hand side of Eq. ( 19) results in the identity matrix, this becomes Inversion of the symmetric square matrix multiplying ẍf and premultiplying by that result and then by ρcA f produces the counterpart to Eq. ( 9) as A new symmetric matrix D sL can now be defined as the counterpart to Eq. (10) for a singular mass matrix as Finally, the augmented DAA equation in this case becomes Note that if M s is nonsingular and G is square and nonsingular, then Eq. ( 22) reduces to Eqs (10) and ( 23) reduces to Eq. (11).Equations ( 22) and ( 23) are the counterparts to Eqs (10) and (11) when the structural mass matrix is singular.Equation ( 23) does involve some additional matrix computation at every time step on the right hand side for the structural coupling vector as opposed to Eq. ( 11).However, once D sL has been determined in the preprocessing phases before the transient response analysis, the left hand side computations remain the same.
There is a profound difference in the two expressions Eqs (10) and (22).To compute D s , the entire structural mass matrix must be factored.If it is diagonal and zero masses occur for dry structure in the model, then these can be eliminated in the inversion as the rectangular G matrix only couples structural nodes on the wet surface.However, if the mass matrix is non -diagonal, then such a process might not be very simple even though it is generally sparse.As it is, the order of the mass matrix is the number of structural degrees of freedom.
On the other hand, D sL involves the factorization of the mass matrix once it has been transformed to the fluid system that only couples to the structure at the wet surface.Hence any singularities in mass for internal dry structure have no influence.As it is, even a zero mass on the wet surface can be handled as long as there are enough non -zero masses nearby to create non -zero entries on the diagonal of the central matrix in Eq. ( 22) that is to be inverted.Of course the factorization in Eq. ( 22) also involves a sparse matrix, but whose order is that of the number of fluid degrees of freedom.As noted earlier, this is typically much less than the number of degrees of freedom in the structural model.
Without going through any of the details, it is easily shown that the augmentation of the curved wave approximation or any higher order DAA formulations results in exactly the same modifications to the DAA 1 equations developed above: the replacement of D s with D sL , and the added operations for processing the structural coupling vector.

Infinite cylindrical shell
The first problem examined here is the infinite cylindrical shell excited by an incident plane step wave.The same half -model with axial plane strain conditions has been used that was described in [6].Both the DAA 1 and the second order modal form of the DAA (DAA 2m ) were studied with a nonsingular diagonal structural mass matrix for comparisons with the standard augmentation method.The DAA 2m modal constant was taken as 1/2, again in agreement with the results presented in [6].The computational model consists of 18 roughly square elements around the half circumference that define both the structure and the overlying fluid.
Figures 1 and 2 show the DAA 1 and DAA 2m structural velocity comparisons respectively at 0, 90, and 180 degrees around the shell and C and R denote circumferential and radial motions respectively.As can be seen, there is virtually no difference in the results except where the plots show a slight thickening of the traces.The next test consisted of the same computational model but in this case, the mass matrix was artificially made singular by setting a translational mass to zero at a point 140 degrees around the circumference from the first point of contact of the incident step wave.This was done only to compute D sL from Eq. ( 22) but was left as it was in the structural equations of motion.
Figures 3 and 4 are the counterparts to those above where the curves marked Old Augmentation are the same as in Figs 1 and 2 for a nonsingular mass matrix using the standard augmentation.It is noted that the responses are affected to some degree, particularly near the source of the singularity as might be expected.However, these results cannot be duplicated by the standard augmentation procedure.
It is interesting that the DAA 2m results in Fig. 4 are less affected by the singularity than are those of the DAA 1 computations shown in Fig. 3. Presumably this is due to the fact that the DAA 2m equations contain more terms on the right hand side than the DAA 1 , hence the structural coupling vector may have somewhat less influence.

Spherical shell
The next example is the spherical shell, again with an incident plane step wave excitation.The computational model used for this analysis is a quarter of a sphere whose surface mesh consists of quadrilaterals only, all of which are very nearly square.It contains 96 elements and 113 node points [6].In that reference, the DAA 2m was used with a modal constant of unity.Here the so-called Curvature Corrected second order Doubly Asymptotic Approximation (DAA 2c ) [8,9] is used.The formulations in these two references differ in some respects but become identical for bodies of constant curvature such as the spherical or infinite cylindrical shell.Both are available in USA for more general usage.
A comparison between the standard augmentation and the new one for a nonsingular diagonal mass matrix is presented in Fig. 5. Velocity responses are shown for 0, 45, 90 and 180 degrees and all are for radial motions.Again the results are virtually indistinguishable.