# II. Fields of Application

## A. General Formulation

This section describes the general formulation used in ATILA to model piezoelectric and magnetostrictive transducersradiating into a fluid, and provides a complete list of the possible types ofanalyses. Each type is then described in more detail in the followingsections.

### 1. Elastic, Piezoelectric and Magnetostrictive Structures

To model a radiating elastic, piezoelectric ormagnetostrictive structure using ATILA, thefinite-element mesh must contain the structure as well as part of the fluiddomain, as shown in the figure above. The mesh of the structure under study anda part of the space surrounding it are necessary when modeling amagnetostrictive structure. This is not necessary for piezoelectric structuresbecause common piezoelectric materials have a high relative permittivity. Thepiezoelectric materials are driven by electrodes (sets of electricallyconnected nodes) generating an electric field. The magnetostrictive materialsare driven by the magnetic field induced by sets of coils. Each coil set is suppliedwith an electric current. The total magnetic field Hcan be broken up as:

H = Hs + Hr

Hs is the magnetic field created by the source currents (the coils are notmodeled) in vacuum:

Hs = Hsp Ip

where Ip is the source current in thepth set of coils. Hsp is pre-determined for each coil using a classicalBiot-Savart numerical integration. Hr is the irrotational component of the magnetic field and canbe given as:

Hr = - grad f

where f is the reduced magnetic potential.

The unknowns are the nodal values of the displacement field U in the whole structure, the electrical potential F in the piezoelectric material, the reducedmagnetic potential f in themagnetostrictive material and in a part of the surrounding space, theexcitation currents I in the excitation coilsand the pressure P in the fluid. Theequations are those of elasticity in the structure, Poisson's equation in thepiezoelectric material, Maxwell's equations for the magnetostatic case (fluxconservation) in the magnetic domain and Helmholtz's equation in the fluid. Theelectromechanical coupling is limited to the piezoelectric and magnetostrictivedomains. There is no electrical interaction between these two domains.Kinematic and dynamic continuity conditions are prescribed on the interfacesbetween the fluid and the structure. An appropriate radiation condition isapplied to the surface boundary that surrounds the fluid domain. The naturalcondition prescribed on the external boundary of the magnetic mesh is a zeroflux condition for the reduced magnetic field. The natural conditionprescribed on the external boundary of a fluid mesh is a zero of the normalderivative of the pressure.

In matrix form, the complete set of equations of the problemcan be written:

(1)

where:

U: vector ofthe nodal values of the components of the displacement field

F: vectorof the nodal values of the electrical potential

f: vectorof the nodal values of the reduced magnetic potential

I: vector ofthe prescribed values of the excitation currents (one component for each coil)

P: vector of the nodalvalues of the pressure field

F: vector ofthe nodal values of the applied forces

q: vector ofthe nodal values of the electrical charges

f: vector ofthe nodal values of the reduced magnetic field flux across the magnetic domain boundary

fb: vectorof the values of the reduced magnetic field flux seen by the coils (one component for each coil)

y: vectorof the nodal values of the integrated normal derivative of the pressure on thesurface boundary S

[Kuu]: stiffness matrix

[KuF]: piezoelectricmatrix

[Kuf]: piezomagneticcoupling matrix

[KuI]: source - structurecoupling matrix

[KFF]: dielectric matrix

[KfI]: source- magnetization coupling matrix

[Kff]: magnetic(pseudo-) stiffness matrix

[KII]: inductance matrix invacuum

[M]: consistent mass matrix

[H]: fluid (pseudo-) stiffness matrix

[M1]: consistent (pseudo-)fluid mass matrix

[L]: coupling matrix at the fluid structureinterface (connectivity matrix)

[0]: zero matrix

w: angular frequency

r: fluid density

c: fluid sound speed

T: transpositionoperator

A large number of analyses can be done using this generalequation. ATILA can solve the following problems (quantities in parentheses are the results):

STATIC ANALYSES

- STA1: computation of the staticresponse of an elastic structure to a concentrated or distributed force(displacement field).
- STA2: computation of the staticresponse of a piezoelectric or magnetostrictive structure to a concentrated ordistributed force (displacement field and electrical potential).
- STA3: computation of the staticresponse of an elastic structure to a forced displacement (displacement field).
- STA4: computation of the staticresponse of a piezoelectric or magnetostrictive structure to a forceddisplacement or a forced electrical potential (displacement field andelectrical potential) or a forced magnetizing current (displacement field andreduced magnetic potential).
- STA5: static analysis of ahydroelastic system, i.e., a structure including a fluid domain (displacementand pressure fields).

  MODAL ANALYSES

- MOD1: modal analysis of an elasticstructure (computation of the eigenfrequencies and the normal modes).
- MOD2: modal analysis of a piezoelectricstructure (computation of the resonance and anti-resonance frequencies and thenormal modes for different electrical boundary conditions).
- MOD3: modal analysis of amagnetostrictive structure (computation of the resonance and anti-resonancefrequencies and the normal modes for different electrical boundary conditions).
- MOD4: modal analysis of a closed fluiddomain under zero pressure and/or zero flux boundary conditions (computation ofeigenfrequencies and normal modes).
- MOD5: modal analysis of a hydroelasticsystem, i.e., a structure including a fluid domain (computation ofeigenfrequencies and normal modes).
- MOD6: modal analyses of a one- or two-or three-dimensional periodic passive structure (for given wave vectors, computationof eigenfrequencies and normal modes)
- MOD7: modal analyses of a one- or two-or three-dimensional periodic active structure (for given wave vectors,computation of eigenfrequencies and normal modes)

  HARMONIC ANALYSES

- HAR1: harmonic analysis of an in-vacuumelastic structure (for a given frequency, computation of the displacement fieldin the structure).
- HAR2: harmonic analysis of anon-radiating piezoelectric or magnetostrictive structure (for a givenfrequency, computation of the displacement field, the electrical potential, thereduced magnetic potential and the electrical impedance of the structure).
- HAR3: computation of the pressure fieldradiated in an infinite fluid space by a vibrating surface where thedisplacement field is known at all points.
- HAR4: harmonic analysis of an elasticstructure driven by a force or a prescribed displacement, radiating in aninfinite fluid space (for a given frequency, computation of the pressure anddisplacement fields).
- HAR5: harmonic analysis of apiezoelectric or magnetostrictive structure electrically driven and radiatingin an infinite fluid space (for a given frequency, computation of the pressurefield, displacement field, electrical potential, reduced magnetic potential andthe electrical impedance of the structure).
- HAR6: harmonic analysis of an elasticor piezoelectric or magnetostrictive structure excited by an impinging wave(for a given frequency, computation of the pressure field, the displacementfield, the electrical potential and the reduced magnetic potential).
- HAR7: harmonic analysis of thediffraction of a plane-wave by a one- or two-dimensional periodic elasticstructure (for a given frequency, computation of the pressure field, thedisplacement field and the reflection and transmission coefficients).
- HAR8: harmonic analysis of a one- ortwo-dimensional periodic piezoelectric structure (for a given frequency,computation of the pressure field, the displacement field, the reflection andtransmission coefficients, the free-field voltage sensitivity or thetransmitting voltage response).
- HAR9: harmonic analysis of a radiatingpiezoelectric or magnetostrictive structure using a coupled FEM-BEM method (fora given frequency, computation of the displacement field and the electricalimpedance of the structure).

  TRANSIENT ANALYSES

- TRA1: transient analysis of anin-vacuum elastic structure driven by external forces (for a given time,computation of the displacement field in the structure).
- TRA2: transient analysis of anin-vacuum piezoelectric or magnetostrictive structure driven by external forces(for a given time, computation of the displacement field, the electricpotential and the reduced magnetic potential).
- TRA3: transient analysis of anin-vacuum piezoelectric or magnetostrictive structure electrically driven (fora given time, computation of the displacement field, the electric potential andthe reduced magnetic potential).
- TRA4: transient analysis of apiezoelectric or magnetostrictive structure electrically driven and radiatingin an infinite fluid space (for a given time, computation of the pressurefield, the displacement field, the electric potential and the reduced magneticpotential).

  The relationship between the excitation magnetic field andthe currents flowing through the sources are calculated by an external programcalled FLUXLABS developed by J.-C. Sabonnadière et al. (Laboratoire d'Electrotechnique de Grenoble, ENSIEG, BP 46,F-38402 Saint Martin d'Hères, France) which applies the Biot and SavartLaw to windings of different geometries.

  ### 2. Modeling of Periodic Structures with 1D or 2D Periodicity

  ATILA models structures having aone- or two-dimensional periodicity for radiation or scattering problems. Inthis case, the space is divided into 3 domains: domains I and III, which aresemi-infinite and homogeneous fluid domains, and domain II, which includes theperiodic structure and a portion of the near-field fluid (see figure below).

  Region II is modeled using finite-elements and the user mustsupply the mesh for a unit cell. In Regions I and III, the pressure field isexpanded on a plane-wave basis, homogeneous or evanescent, and is automaticallygenerated. The matrix formulation of this problem includes the continuityconditions between different way of representing the pressure field at the I/IIand II/III interfaces as well as Bloch-Flocquet conditions between successiveunit cells. The computation of transmission and reflection coefficientsassumes that only the specular components propagate to infinity.

  ### 3. Modeling of Periodic Structures with 1D, 2D or 3D Periodicity

  ATILA models structures having aone-, two- or three-dimensional periodicity and no losses for eigenvalueproblems. For a given wave vector, the code computes eigenfrequencies andnormal modes. The angular frequency w is a periodical function of the wavevector; thus the problem can be reduced to the first Brillouin zone. Dispersioncurves can be built varying the wave vector on the first Brillouin zone.

  ### 4. Modeling Internal Losses

  In ATILA, losses in active orpassive materials can be taken into account by using complex physical constantssuch as:

  s = s' + js"

  for a physical constant s. This is possible for nearly allthe elements (see Chapter IV). The user must ensure the coherence of the tensors (standard inequalities between real and imaginary parts. See for example: R. HOLLAND, "Representation of Dielectric, Elastic and Piezoelectric Losses by Complex Coefficients," IEEE, SU14,18-20 (1967)). Furthermore, the user must take into account the changes inthese parameters with frequency, if necessary, by updating the data file.

  In the following sections, each type of analysis isdescribed in detail. Part IV describes some test examples furnished with the ATILA code.

  ### 5. Transient Analysis Methods

  In this analysis, the matrix equation (1) becomes:

  where × and ×× denote the first and second time derivative, respectively. K¢ and K² denote the realand imaginary parts of K respectively. w0is the pulsation at which the materials losses are provided. The preceding system of equations may be rewritten as follows:

  (2)

  This differential equation is solved by an iterative method, taking a constant time step Dt; the successive values of X and Y are noted Xn and Yn and correspond to values calculated at the time n D t. Three methods are implemented in ATILA: the Central Difference Method, the Newmark Method and the Wilson-q Method.

  Both methods are based on the same algorithm to calculate Z:

  Z and the values of a1, a2, b1, b2, k, l and h depend on the chosen method and are explained below.

  #### a) The Central Difference Method

  This method consists in using a second order (parabolic) interpolation along the time axis and deriving the first and second orderderivative from it:

  Replacing these values in the equation (2) above written attime step n leads to the following parameter values:

  Z = Xn+1/Dt2,a1,= ½, a2, = b1,= b2, = k = l = 0and h = Dt.

  Assuming that the initial first derivative of X is zero, the algorithm is supplied with an initialvalue of X-1:

  This method is conditionally stable, i.e., knowing thehighest eigenfrequency fmax of the lossless system ofequations, the time step must satisfy:

  #### b) The Newmark Method

  This method consists in using a truncated Taylor seriesexpansion of Xn+1 and its firstderivative:

  (3)

  The parameters b and g may be arbitrarily chosen. Introducing these equations inequation (2) above written at time step n+1 leads tothe following parameter values:

  , a1,= g, a2, = b, b1,= (1-2b)Dt2/2, b2, = (1-g)Dt, k = l = 1 and h = Dt.

  Assuming that the initial first derivative of X is zero, the algorithm is supplied with an initialvalue of :

  This method will be unconditionally stable, provided thatthe following inequalities are satisfied:

  where fmax is the highesteigenfrequency of the lossless system of equations. Note that a value of g > ½ introduces a numerical damping, which is generallyavoided. Names are given to usual couples of (g, b): (½, 0) stands for the Central Explicit DifferenceMethod, (½, ¼) for the Average Acceleration Method, (½,) for the Linear Acceleration Method.

  #### c) The Wilson-q Method

  This method is very similar to the Newmark Method, exceptthat the equations (3) are first written for the time step n+q, q³1, instead of n+1 (Dt is replaced with qDt) and the parameters g and b are respectively set to ½ and . Once Xn+q and its first and second time derivatives are calculated as above, a linearinterpolation between Xnand Xn+q gives thesolution Xn+1:

  This leads to the following parameter values:

  , a1,= q/2, a2,= q2/6, b1,= q2Dt2/3, b2 = qDt/2,k = 1, l = q and h = qDt.

  The stability of this algorithm is more difficult todetermine, but is proven for values of q greater than1.366.

  Important note:

  Extending these algorithms to piezoelectricity ormagnetostriction may lead to inaccurate results, because of the quasi-staticnature of electrical degrees of freedom: the matrix [A] of equation (2) becomessingular, fmax tends to infinity. In the case of theCentral Difference Method, the implemented algorithm can deal with this,provided that the loss matrix [B] contains elastic losses only, i.e., and are both null. In the case of the Newmark Method, a value of gslightly greater than 0.5 will help in damping the spurious oscillationsrelated to electric transients.