INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 46, 1997}2010 (1999) EXTENDING INTEGRAL EQUATION TIME DOMAIN ACOUSTIC SCATTERING ANALYSIS TO LARGER PROBLEMS V. SUCHIVORAPHANPONG, S. P. WALKER* AND M. J. BLUCK Mechanical Engineering Department, Imperial College of Science ¹echnology and Medicine, Exhibition Road, ¸ondon S=7 2BX, ;.K. SUMMARY A method is presented to accelerate the execution of integral equation time domain analyses of exterior acoustic scattering problems. Conventionally, these have costs which scale with the "fth power of the frequency of the excitation, and practical limits to such computations are reached when bodies approach perhaps &5}10 wavelengths long. The fast approach presented is based on exploiting the pulsed nature of the illumination to omit much nugatory calculation. There is an associated slight accuracy loss; this is investigated. The method has costs which can scale with frequency to the power as low as &3, such that, for example, costs on a 18)5 wavelength body are reduced by a factor of about 28, with this factor itself increasing with roughly the square of the body size. Associated with the reduction in operations is a reduction in the scaling of storage required, from the third to the second power of frequency. Examples of analysis of large scatterers are presented, extending to a &22 000 node &almond'. Copyright ( 1999 John Wiley & Sons, Ltd. KEY WORDS: time domain integral equations; acoustic scattering; wave propagation; transient scattering; acoustic cross-section 1. INTRODUCTION Large acoustic scattering problems arise in a number of engineering applications where it is desired to predict the acoustic "eld scattered by some body. For very large, smooth scatterers this can be tackled using ray tracing and similar techniques, whilst for small scatterers (up to say a handful of wavelengths long) full-"eld solutions are computationally practicable. However, these full-"eld solutions have computational costs, both in operations and storage, which scale very severely with problem size (a point to which we will return). There is a considerable range of problem, from perhaps a few to a few tens of wavelengths long, or where signi"cant features are small on an otherwise large body, where the "eld solutions are too expensive, but the ray-based techniques can be insu$ciently accurate. One approach for such analyses is the integral equation time domain technique [1}3]. In this paper we will describe a modi"cation to the usual integral equation time domain technique which * Correspondence to: S. P. Walker, Mechanical Engineering Department, Imperial College of Science Technology and Medicine, Exhibition Road, London SW7 2BX, U.K. E-mail: s.p.walker@ic.ac.uk CCC 0029-5981/99/361997}14$17.50 Copyright ( 1999 John Wiley & Sons, Ltd. Received 17 August 1998 Revised 30 March 1999 1998 V. SUCHIVORAPHANPONG, S. P. WALKER AND M. J. BLUCK reduces its cost and cost scaling sharply, and extends considerably the range of problem size tractable. Recently, an approach was presented [4, 5] which was able to reduce cost scaling by one power of frequency; the present one, based on work recently presented in electromagnetic scattering [6] applications, improves on this, reducing the cost scaling by approaching two powers. In the remainder of this introduction we will discuss the origin of the high computational costs of scattering calculations, "rstly by reference to the more widely employed integral equation frequency domain and "nite di!erence time domain approaches. In Section 2 we will outline the integral equation time domain method and its discretization, and demonstrate the origin of the cost scaling of its conventional form, as a precursor to discussing the fast method. In Section 3 we present modi"cations we have developed to reduce costs, and describe their practical incorporation into an integral equation time domain treatment. Section 4 will present results from its use. We seek to solve the scalar wave equation in the exterior region surrounding some hard scatterer (i.e. a homogeneous Neumann boundary condition), where some incident wave impinges upon the scatterer. One approach is to employ a frequency domain integral equation treatment. The di!erential equation is transformed to an integral equation on the surface of the body [7]: P) 2n/(r)"4n/*/#(r)# L /(r@) A B R ) n@ 1 #jk e+kR ds@ R2 R (1) Here the "eld at a point on the surface is given as a function of the (unknown) "elds at all other points on the surface. The surface must be discretized. Nodal spacing would depend on the highest frequency f of interest, with spacing typically &1/10 of the shortest relevant wavelength employed. With only the surface discretized, the number of nodes scales with f 2. In discretized form (1) becomes a dense matrix equation of size nodes by nodes, and thus size proportional to f 2. Direct solution then has a cost scaling with f 6. Storage of this dense matrix is also di$cult. As an example, a 10 wavelength diameter sphere with &10 nodes/wavelength, would have about 30 000 nodes, and an associated matrix occupying approaching 10 Gb. Recent work, especially in the computationally similar area of electromagnetic scattering, has applied &fast multipole' techniques to this matrix, with impressive results [8]. Alternatively, the di!erential equation can be tackled directly in the time domain (with either harmonic on pulsed time domain illumination, with frequency domain results extracted via Fourier transformation, if required). Using a "nite di!erence approach, we would discretize a region of the space surrounding the scatterer (truncating the region some way from the scatterer, and applying some approximate absorbing boundary condition). The volume discretization now results in the total number of nodes scaling with f 3. With the number of timesteps for which we must analyse similarly proportional to f we obtain a nominal cost scaling of f 4. In practice, as problems become bigger, more nodes per wavelength are required to maintain phase accuracy as disturbances propagate through a larger mesh [9]; with three spatial and one temporal dimension a!ected this raises the cost scaling to rather more than the fourth power of frequency. Storage is again also a di$culty, but less so than for the integral equation treatment. Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 46, 1997}2010 (1999) 1999 TIME DOMAIN SCATTERING ANALYSIS TO LARGER PROBLEMS 2. THE INTEGRAL EQUATION TIME DOMAIN APPROACH The time domain equivalent of (1) is 2np(r, t)"4np */# P ) R3 C p(r@, t*)# c Lt@ (r@, t*)D ds@ R)n (r, t)# R Lp (2) L 2.1. Numerical treatment A fuller description of this has been presented elsewhere [2, 10]; we will only summarize here. The surface of the scatterer is represented via curvilinear quadratic surface patches, of (say) eight nodes each, over which the geometry is approximated using polynomial functions S (m, g), a"1, . . . , 8 (3) a where m, g are parameterized intrinsic co-ordinates. The position vector to a point on such a patch is then 8 rm(m, g)" + S (m, g)r (4) a j(m,a) a/1 where j"j(m, a) are the node numbers of the local nodes on element m, and the position vector of each of these nodes is r . The pressure variation is represented similarly j(m,a) 8 (5) p (m, g; t)" + S (m, g)pm (t) b m a a/1 and the temporal variation of the "eld is modelled using one-dimensional quadratic Lagrangian basis functions with a nodal separation *t. Inserting these representations into (2) allows us to write 2np(r , k*t)"4np (r , k*t) i */# i M n@ (m, g) ) R 1 8 3 m m #+ + S (m, g) + ¹ (q(R , k))pml ab a b m R2 R m m a/1 m/1 m,g b/1 PP G M #+ m/1 C C PPm,g G m DH DH n@ (m, g) ) R 1 8 3 m + S (m, g) + ¹Q (q(R , k))pml ab a b m R2 c m a/1 b/1 DJ(m, g)Ddm dg DJ(m, g)Ddm dg (6) Integrals are evaluated using Gaussian quadrature. Some elements involve singularities, and are partitioned to integrate these accurately [11]. Doing this, we obtain 2np(r , k*t)"4np (r , k*t) i */# i 8 3 M N1 N' N' n@ (m , g ) ) R m q r m + + S (m , g )¹ (q)pml DJ(m , g )Du u #+ + + + ab a q r b q r q r R3 m a/1 b/1 m/1 p/1 q/1 r/1 M N1 N' N' n@ (m , g ) ) R 1 8 3 m q r m + + S (m , g )¹Q (q)pml DJ(m , g )Du u #+ + + + a q r b ab q r q r R2 c m a/1 b/1 m/1 p/1 q/1 r/1 (7) G G Copyright ( 1999 John Wiley & Sons, Ltd. C C DH DH Int. J. Numer. Meth. Engng. 46, 1997}2010 (1999) 2000 V. SUCHIVORAPHANPONG, S. P. WALKER AND M. J. BLUCK where N is the number of partitions (1 usually; 2 or 4 for singular integrals), N is the order of the 1 ' Gauss quadrature, u is the qth Gauss weight and s is the parametric retarded time t*. q Some (most) quantities p on the right-hand side are known; the associated retarded time was longer ago than one timestep. Such terms, matrix B below in (8), are multiplied out to form a vector say McN. Coe$cients multiplying unknown (nearby and present) values are gathered, and together with the free term form a sparse matrix equation for the "eld at time k: [A]Mp N"[B]Mp N"McN (8) k j,j:k The number of non-zero terms on any row of A is independent of f; for large problems the vast majority of the surface lies further away than c*t. The cost of repeated solution of this matrix equation is small. As discussed further below, the main work of the method is associated with the formation of the coe$cients in B (and A), and the multiplications by B to form McN. 2.2. Cost scaling in the conventional IE¹D Equation (2) and its discretized counterpart (7) show the "eld at any point on the surface, at some time, to be an integral of (some geometrical function of ) "eld values elsewhere on the surface at the associated retarded times. Again discretized with some frequency-dependant re"nement we have the number of surface nodes and patches proportional to f 2, such that the cost of evaluation of the "eld at any point and at any timestep is proportional to f 2. Doing it at each node thus generates a cost at each timestep proportional to f 4, and with again the number of timesteps proportional to f the overall cost scaling is seen to be with f 5. This cost scaling is implicit in the summations of equation (7). Storage needs scale as, but are rather greater than, those of the frequency domain integral equation approach. The coe$cients of B in (8) represent a dense matrix, of size a few times nodes by nodes. As in the frequency domain, for multi-wavelength problems this is far too large to store. The approach of necessity adopted is the recalculation of the matrix coe$cients as required. Cost scalings are unchanged, but a considerable cost penalty is naturally incurred. In exchange storage needs are vastly reduced. The next largest requirement is to store one transit's worth of the "eld history. With nodes scaling with f 2, and the timesteps in a transit with f, this scales with f 3. For the example sphere mentioned earlier, it comprises about 30 Mb. 3. THE FAST FORM OF THE IETD 3.1. Fast approximation The reduction in cost involves introduction of an approximation based upon the behaviour of high-frequency transient "elds. For a large body, illuminated with a short transient pulse, the pulse will tend to propagate over the body as a fairly narrow excited region, with perhaps a modest &wake', and with propagation of other similarly narrow active regions spawned by re#ections from geometrical discontinuities. There are two respects in which this can be exploited to save costs. (i) Over most of surface locations, most of the time, the "eld will be (close to) zero. If we can identify these locations and times we need not evaluate at them. If the excited region has Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 46, 1997}2010 (1999) TIME DOMAIN SCATTERING ANALYSIS TO LARGER PROBLEMS 2001 a dimension related to the pulse width (&1/f ), the fraction of the surface which is excited at any instant should decline with f. If so, this provides a reduction in cost scaling by one power of frequency. (ii) Where the "eld is to be found, it is found by integrating contributions from all the rest of the surface evaluated at the appropriate retarded times. By the same arguments as above, most of these integrations will be over regions in which the historical "elds were similarly (almost) zero. Similarly, the fraction of the surface exhibiting non-zero historical "elds should decline with frequency, providing a reduction in cost scaling by one further power of frequency. These savings will of course be an upper bound; the relevant regions are unlikely to be identi"ed perfectly, and there will be some &housekeeping' cost associated with their identi"cation. Note that one implication of the above is that interactions between most surface-point pairs are never actually calculated; only a fraction of the coe$cients of the full matrix is ever evaluated during the entire course of the analysis. The tiny fraction needed at each timestep is calculated as needed. Surface "eld history is again a major storage requirement, but this fast approach reduces this, and its scaling, very sharply. At any one location only those periods of history which are excited need of course be stored; a modest, frequency-independent number of values. Total history storage should then scale with the number of nodes; frequency squared. This is the same scaling indeed as merely the mesh geometry itself, and in practice mesh-related storage actually is the dominant requirement. 3.2. Computational implementation There are two aspects to achieve the cost saving; identifying where to calculate the "eld, and when doing so, identifying over which portion of the surface to integrate. We have investigated employing (rather complicated) searching algorithms to identify where to integrate. However, in practice a simple brute force approach is actually very economical, and accurate. For each location at which we are calculating the "eld, we simply test the retarded "eld magnitude in every element. This is both a very quick test, and is done for only the centre node of an element, so the total time spent on this testing is a tiny (O(1 per cent)) fraction of the total time. For an explicit treatment, a dynamically adaptive selection of locations at which to calculate would be simplest to implement: if at some time the "eld is found to be signi"cant at a point, it is calculated at neighbouring points till a halo of insigni"cant "eld values is discovered. In the implicit treatment it is not possible to "nd a "eld value in isolation from its neighbours, so this cannot be done. It is essentially always the case that the surface is active in the region where the incident wave impinges, and this is used as one criterion for the selection of where to calculate. Additionally, the progress of "elds as they move over the surface of the body is tracked. They propagate at the wave speed, and so are captured by including all locations within a distance &c *t' of those locations found to be active at the last step. These two approaches between them accommodate all convex geometries. They do not cope with propagation across cavities or between multiple scatterers; these are treated by identifying a small number of &watch nodes', where the "eld is calculated always, with these watch nodes becoming nuclei for the growth of calculation regions whenever a "eld is calculated there. Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 46, 1997}2010 (1999) 2002 V. SUCHIVORAPHANPONG, S. P. WALKER AND M. J. BLUCK Even using these fast methods, multi-wavelength scattering calculations are computationally large, and parallel computers need to be used. The fast algorithm is rather more complex to parallelize than conventional IETD. The nodes at which calculations are performed di!er at each timestep, and the elements over which integration must be performed varies from node to node, and from timestep to timestep. This makes load balancing rather awkward, but good balancing can still be achieved. (In part, by the deliberate random renumbering of nodes and elements, in contrast to "nite element treatments!). Most of the larger problems reported in this present paper were analysed on parallel machines; either a Cray T3D (512 processors) or a Fujitsu AP3000 (64 processors). 4. EXAMPLES OF PERFORMANCES 4.1. Accuracy There is naturally some loss of accuracy associated with neglect of calculation at and integration over some locations. Obviously, the savings achieved and the accuracy loss are linked, Figure 1. Surface pressure on front-centre of sphere, for indicated thresholding levels, with zoomed inset (diameter 1, pulse half-width 0)4, 1538 node mesh) Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 46, 1997}2010 (1999) TIME DOMAIN SCATTERING ANALYSIS TO LARGER PROBLEMS 2003 and what accuracy loss is tolerable depends on the application, but as the following examples will demonstrate, it is generally small at levels of approximation which still provide marked cost reductions. In Figure 1 we show the time dependant surface pressure on a unit diameter sphere, comprising 1538 nodes, illuminated with a Gaussian pulse of half-width 0)4. (In this, as in all cases, results are normalized: body dimensions and pulse widths are expressed in arbitrary units, pressures are expressed as a multiple of peak incident pressure, itself unity, and the wave speed is unity). This is done with a range of values of the threshold pressure below which calculation and integration are neglected, with results compared also to the analytical solution. As is seen, there is a steady decline in accuracy with threshold, but even with a threshold of 5 per cent of the incident wave magnitude the maximum error is &3 per cent of the peak "eld. A frequency domain surface pressure can be extracted from this time domain result, and this is used to obtain the bistatic acoustic cross section shown in Figure 2. The frequency extracted corresponds to the sphere being 1 wavelength in diameter. As before, little error is introduced by 2 the thresholding. Figures 3 and 4 show similar tests on an almond-shaped target [12]. There is no analytical result for comparison, but the modest degradation with threshold is again apparent. The (greatly) Figure 2. Acoustic cross-section of sphere, for indicated thresholding levels (wavelength 0)4 diameters, 1538 node mesh) Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 46, 1997}2010 (1999) 2004 V. SUCHIVORAPHANPONG, S. P. WALKER AND M. J. BLUCK Figure 3. Surface pressure on-axis, rear surface of almond, for conventional (0 per cent) and indicated thresholding levels, with zoomed inset (plane wave incident along axis on tip, pulse half-width 0)674, almond length 10 2450 node mesh) zoomed inset in Figure 3 shows decaying oscillatory behaviour after the unit magnitude incident pulse has passed; having a peak amplitude of only &0)003, this oscillation is lost in the thresholded analysis. One of very few large published time domain computations was recently presented by Ergin [3], showing the surface "eld history on a 10 wavelength long almond. The case studied is reproduced here as Figure 5. We have access only to the printed graph, but a few points measured from the graph are included in Figure 5. Whilst this is only a comparison between codes, rather than with experiment or analytical solution, the reasonable agreement does provide some comfort. 4.2. Cost savings and scaling We now turn to evaluate the cost saved in exchange for the introduction of errors such as those discussed. Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 46, 1997}2010 (1999) TIME DOMAIN SCATTERING ANALYSIS TO LARGER PROBLEMS 2005 Figure 4. Near "eld for the almond (location 0, 3)5, 0 relative to almond de"nition co-ordinates) for conventional (0 per cent) and indicated thresholds (plane wave incident along x-axis on tip, pulse half-width 0)674, almond length xx, 2450 node mesh) Figure 6 shows the computational cost of analysing propagation over a series of spheres, ranging in diameter up to nine wavelengths, and comprising up to 23 554 nodes. Timesteps, nodal separations, pulse widths and frequencies extracted are changed in concert between spheres, to maintain similar &nodes per wavelength' and so on. Details are summarized in Table I. Thresholding of 5 per cent is employed. Computational cost comprises overwhelmingly ('95 per cent) activities directly proportional to the number of surface patch integrations performed, allowing this to be used as an accurate measure over the wide range of machine types necessarily employed. The graph shows clearly the large savings being achieved, with the cost for the largest sphere reduced by a factor of about 25, and with that factor increasing sharply with body size. (Note that the larger spheres cannot be analysed using the conventional code, but accurate cost ratios are available merely by noting the fraction of the total work (integrations) which the fast method actually does). The logarithmic inset allows the scaling of cost with frequency to be inferred; depending on how it is extracted, it ranges between about 2)8 and 3)2; more or less the third power discussed above. The memory requirements of these same cases are shown in Figure 7. The elimination of the history storage scaling with frequency cubed has a marked e!ect, with the factor by which storage Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 46, 1997}2010 (1999) 2006 V. SUCHIVORAPHANPONG, S. P. WALKER AND M. J. BLUCK Figure 5. Comparison of the results of Ergin [3] (black circles) with a 1 per cent threshold case. Almond, 3)75 pulsewidths (&10 wavelengths) long is reduced being about 2 for the 9 wavelength case, with this factor increasing roughly linearly with frequency. Figure 8 shows a comparison of the conventional and fast approaches for a series of almonds, ranging in size up to 21 552 nodes, and 18)5 wavelengths. (Details are summarized in Table II). Similar behaviour to the spheres is seen, with costs reduced greatly; here by a maximum factor of &28 in the 18)5 wavelength case. Cost scaling is again to the power 2)9}3)2 of frequency. A "gure is not shown, but savings in storage are also similar to those of the spheres; for the 18)5 wavelength case storage needs are reduced by a factor of about 2. 5. DISCUSSION AND CONCLUSIONS The modi"cation to the usual IETD has been shown to reduce costs greatly, and to reduce the scaling with frequency of these costs by about two powers of frequency. On the biggest (18)5 wavelength) example shown here, the factor by which costs are reduced is &28. These reductions Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 46, 1997}2010 (1999) 2007 TIME DOMAIN SCATTERING ANALYSIS TO LARGER PROBLEMS Figure 6. Computational cost (expressed in terms of integrations required) versus body size for a series of spheres, conventional and 5 per cent thresholding Table I. Unit radius sphere mesh and illumination information. The pulsewidth parameter g gives the full-width at half-maximum of the incident Gaussian plane wave pulse Number of nodes 2402 4706 6146 9602 16 226 23 554 Bodysize (wavelengths) Timestep (and biggest nodal separation) Ave. nodal separation g 2)765 3)870 4)423 5)529 7)187 8)659 0)0879 0)0657 0)0546 0)0459 0)0385 0)0364 0)0723 0)0517 0)0452 0)0362 0)0278 0)0231 0)3516 0)2628 0)2190 0)1800 0)1529 0)1442 in operations are accompanied by reductions in the scaling of memory needed from the third to the second power of frequency. There is an associated slight loss in accuracy, but costs savings of the size mentioned are typically associated with maximum errors of only a few per cent of peak "eld magnitudes. Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 46, 1997}2010 (1999) 2008 V. SUCHIVORAPHANPONG, S. P. WALKER AND M. J. BLUCK Figure 7. Storage requirement versus body size for a series of spheres, conventional and 5 per cent thresholding Table II. NASA almond (9)936 units long) mesh and illumination information (timestep size"1)323]average nodal separation) Number of nodes 1202 2962 5538 10 866 21 552 Bodysize (wavelengths) 4)380 6)868 9)400 13)150 18)510 Copyright ( 1999 John Wiley & Sons, Ltd. Timestep size Biggest nodal separation Ave. nodal separation g 0)3000 0)1914 0)1399 0)1000 0)0710 0)3657 0)1814 0)1639 0)1019 0)0895 0)2270 0)1447 0)1057 0)0755 0)0537 1)2000 0)7655 0)5595 0)4000 0)2840 Int. J. Numer. Meth. Engng. 46, 1997}2010 (1999) TIME DOMAIN SCATTERING ANALYSIS TO LARGER PROBLEMS 2009 Figure 8. Computational cost (expressed in terms of integrations required) versus body size for a series of almonds, conventional and 1 per cent thresholding. Inset: logarithmic plot ACKNOWLEDGEMENTS It is a pleasure to acknowledge helpful discussions with Dr. S. J. Dodson, previously of this Department. REFERENCES 1. Bennett CL, Mieras H. Time domain integral equation solution for acoustic scattering from #uid targets. Journal of the Acoustical Society of America 1981; 69:1261}1265. 2. Bluck MJ, Walker SP. Analysis of three dimensional transient acoustic wave propagation using the boundary integral equation method. International Journal for Numerical Methods in Engineering 1996; 39:1419}1431. 3. Ergin AA, Shanker B, Michielssen E. Transient analysis of acoustic scattering using marching on in time with plane wave time domain algorithm. Presented at 14th Annual Review of Progress in Applied Computational Electromagnetics, Monterey, CA, vol. 2. Applied Computational Electromagnetics Society: Monterey, CA, 1988; 866}872. 4. Walker SP, Lee BH. Reduced-cost methods for large time domain integral equation scattering analyses. Communications in Numerical Methods in Engineering 1998; 14:751}761. 5. Lee BH. Improved integral equation methods for transient wave scattering. Ph.D. ¹hesis. Mechanical Engineering Department, Imperial College, University of London, 1996. 6. Dodson SJ, Walker SP, Bluck MJ. Costs and cost scalings in time domain integral equation analysis of electromagnetic scattering. IEEE Antennas and Propagation Magazine 1998; 40:12}21. Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 46, 1997}2010 (1999) 2010 V. SUCHIVORAPHANPONG, S. P. WALKER AND M. J. BLUCK 7. Morse PM, Feshbach H. Methods of ¹heoretical Physics (1st edn). McGraw-Hill: New York, 1953. 8. Song JM, Chew WC. Requirements scaling properties in large scale computing. Presented at IEEE Antennas and Propagation Society International Symposium, Atlanta, GA, vol. 3. IEEE Press: New York, 1998; 1518}1521. 9. Cangellaris AC, Lee R. On the accuracy of numerical wave simulation based on "nite methods. Journal of Electromagnetic =aves Applications 1992; 6:1635}1653. 10. Bluck, MJ. Integral equation methods for transient wave propagation. Ph.D. ¹hesis. Mechanical Engineering Department, Imperial College, University of London, 1993. 11. Lachat JC, Watson JO. E!ective numerical treatment of boundary integral equations: a formulation for threedimensional elastostatics. International Journal for Numerical Methods in Engineering 1976; 10:991}1005. 12. Woo AC, Wang HTG, Schuh MJ. Benchmark radar targets for the validation of computational electromagnetics programs. IEEE Antennas and Propagation Magazine 1993; 35:84}89. Copyright ( 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 46, 1997}2010 (1999)

1/--страниц