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.
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. H
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
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
MODAL ANALYSES
HARMONIC ANALYSES
TRANSIENT ANALYSES
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.
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.
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.
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.
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 nDt. 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.
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:
This method consists in using a truncated Taylor seriesexpansion of Xn+1 and its firstderivative:
(3)
The parameters b and
, a1,= g, a2, =
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 ) for the Linear Acceleration Method.
This method is very similar to the Newmark Method, exceptthat the equations (3) are first written for the time step n+. 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.
![]() |
![]() |
![]() |