Scientific Research

An Academic Publisher

The Third-Order Viscoelastic Acoustic Model Enables an Ice-Detection System for a Smart Deicing of Wind-Turbine Blade Shells ()

**Author(s)**Leave a comment

KEYWORDS

Cite this paper

*Journal of Applied Mathematics and Physics*,

**4**, 1949-1976. doi: 10.4236/jamp.2016.410197.

1. Introduction

In cold weather, a layer of atmospheric ice (AI) is accreted on the outer surfaces of the blade shells (BSs) of a wind turbine. As is well known, this layer can cause unexpected downtime and increased maintenance cost of the turbine, thereby resulting in reduced efficiency of the energy production by the turbines.

Remark 1.1. All of the cold-climate areas can be found within the Köppen-Geiger climate classification system (http://en.wikipedia.org/wiki/Köppen_climate_classification). The climates associated with the freezing temperatures comprise various weather conditions, from the hot summer continental climate (see also http://en.wikipedia.org/wiki/Humid_continental_climate) with the minimum temperatures in the coldest months down to −10˚C to the tundra climate (see also http://en.wikipedia.org/wiki/Tundra) with the minimum temperatures in the coldest months down to −50˚C.

The coldest climate is the ice cap climate (see also http://en.wikipedia.org/wiki/Ice_cap_climate) with the winter temperatures below −50˚C in the Arctic (see also http://en.wikipedia.org/wiki/Polar_climate and http://en.wikipedia.org/wiki/Climate_of_the_Arctic). The minimum temperatures during a year can vary between −50˚C and −90˚C (http://en.wikipedia.org/wiki/Polar_climate and http://en.wikipedia.org/wiki/Climate_of_Antarctica).

Thus, the temperatures dealt with in connection with deicing of the BSs are in the interval from 0˚C down to the aforementioned values depending on a specific climate.

□

Note that all of the physical-quantity values are specified below in the SI without indication of the corresponding units, with the exception of the cases where the units are not the SI ones (such as the temperature values in Remark 1.1).

The accreted AI is generally of different qualities. It can be continuous or porous, more specifically, hexagonal crystalline or low-density amorphous, clear, frozen dense snow, soft or hard rime, glaze, etc. In terms of continuum mechanics, these qualities are distinguished with values of the material parameters of the AI (e.g., [1] , [2] ).

According to [2] , the AI (including porous and continuous ice) is regarded as porous medium where pores are filled with air. The corresponding (volumetric) density of the AI mass is expressed as follows

(1.1)

where is the AI porosity, i.e., the volume fraction occupied by air (),is the mass density of a continuous, non-porous ice, and is the mass density of air. All of these parameters present the equilibrium values. Note that

(1.2)

(1.3)

As follows from (1.1),

(1.4)

The mass densities of the seasonal and dense snows are below and above 300, respectively.

Nowadays, BSs are usually fabricated of composite materials. As is well known (e.g., [3] ), they are viscoelastic.

Viscoelasticity is also a property of AI. A viscosity of a medium presents the product of the corresponding elastic modulus and the relaxation time of the quasi-eq- uilibrium (QE) component of the Cauchy stress matrix to its equilibrium value. More specifically, if and are the bulk and shear moduli of an isotropic material, then parameters

(1.5)

are the volume (or compressional) and shear viscosities. The stress relaxation exists in any material medium, in gases, liquids, and solids, no matter if the medium is spatially non-homogeneous or spatially homogeneous. Also note that, for the sake of simplicity, the stress-relaxation time (SRT) is assumed to be a scalar rather than a matrix. Also, the same value of it is used in expressions (1.5) for both the viscosities. The physical picture for the viscosities of solids and the relevant theoretical-physics modeling approach are available (e.g., [4] , §34). Many aspects of theoretical analysis and ex- perimental study of viscoelastic materials are described in [5] and [6] .

The fact that ice is viscoelastic rather than purely elastic (i.e., inviscid) was re- cognized more than a century ago [7] - [10] . This line of research was continued in [11] . It appeared that solid ice has high values of viscosity. For instance, the experimental data in ( [11] , the last column of Table 3) show that the ice shear viscosity values are in the interval. The shear-viscosity values for glaciers obtained for a number of the Swiss glaciers are ( [10] , Table 1). This is confirmed with the fact that ( [12] , p. 305) values of the viscosity of polycrystalline ice vary in the range.

Basing on the pioneering work of J. C. McConnel published in 1891, the ex- perimental data of ( [10] , Table 2) show that the ice viscosity rapidly increases from to when the ice temperature decreases from −1.7˚C to −15.3˚C. These data can partly explain the spread of the numerical values in each of the aforementioned two intervals. “J. C. McConnel appears to have been the first to observe the exact conditions under which ice is capable of being deformed without fracture by stress. He showed that a crystal of ice can be sheared by very small stress in a direction at right angles to the optic axis, and that the rate of shear becomes greater as the stress is increased. His two main conclusions are 1) that the friction between the particles of ice along the shear planes becomes greater as the temperature falls, 2) that when the molecules of ice slide over each other the cube of the friction varies as the square of the velocity.” ( [10] , p. 253). Work [10] also notes that the shear viscosity of common ice at 0˚C is. Why the effective viscosities of glaciers are noticeably greater than the one of common ice remains unclear ( [10] , p. 259).

The following remark indicates a possible interval of the AI SRT.

Remark 1.2. The aforementioned experimental values for ice show that shear viscosity in (1.5) can be in the interval. A typical value of the shear modulus for this ice is (see [13] , Table 1). Applying these data to the first equality in (1.5), one obtains that SRT for continuous ice varies in the range, i.e., in five orders.

In the case of polycrystalline ice (e.g., [12] ), the range of the SRT values is even wider. Indeed, as already noted, values of the viscosity of polycrystalline ice vary from about to ( [12] , p. 305). Combining a typical value of the shear modulus for this ice (see [12] , Table 4.1) with the second equality in (1.5), one obtains that the indicated interval corresponds to the SRT values in the six-order range. □

According to Remark 1.2, the SRT values for continuous or polycrystalline ice can vary in five-six orders. The latter feature practically means that the corresponding viscoelastic acoustic model must be relevant at any value of the material SRT.

The AI layer accreted on it the outer surface of a BS is different at different locations on the surface and at different time points. More specifically, the AI landscape varies in both the space and time.

To prevent the losses mentioned at the beginning of this section, one usually applies deicing methods. The most common of them is heating. However, by now, all of the deicing techniques are based on a priori assumptions on the AI-layer parameters over the BS surface. In the case of heating, this often results in under-heating, i.e., reduced reliability of the deicing, or over-heating that damages the BS-material and, thus, reduce the cost-efficiency of the deicing.

The notion of a smart deicing system (SDS) presumes that the system is smart in the sense that is cost- and energy-efficient, safe to the BS material, and reliable. As noted above, the AI layer is space-time varying. The BS, which the AI landscape is accreted on, has a complex, curvilinear shape, which is generally space-time varying as well, due to the action of the operational load (e.g., [14] ). To be able to smartly remove the above spatiotemporal landscape of the AI, an SDS must provide the heat levels, which are individual to the landscape regions and specific time points. Consequently, these levels should be determined on the basis of the spatiotemporal landscape of the AI-layer material parameters including the layer thickness. To provide this landscape, each heating element (e.g., [15] ) should be accompanied by an appropriate sensor, which is located near, but outside the area of, the heating element and measures a mechanical variable that allows to identify the mentioned parameters of the AI layer. The heating elements can be implanted in each of the three blades of a wind-turbine rotor in such a way that they are distributed over the BS surface uniformly (say, one element per four square meters), whereas the sensors can be attached to the inner surface of the BS at their aforementioned locations. Thus, an SDS is supposed to consist of a sufficiently large array of the identical heating elements (e.g., square- or disk-shaped) and com- munication hardware that controls the array. Similarly, a related IDS is supposed to consist of the corresponding sensor array and communication hardware that controls this array. Both the controls can be implemented by a personal computer (PC).

To meet temporal variations of the AI landscape, both the SDS and IDS should operate in the real-time mode. Both the SDS and IDS have components implanted in, or attached to, the BSs. The BSs are rigidly attached to a wind-turbine rotor, which is in general rotating, and, as follows from the published review works, both the SDS and IDS must not have components located outside the body of a wind turbine. Therefore, both the systems must be located inside the rotating rotor but accessible to a remote operator at any time. For this reason, both the controls should be wireless. Moreover, the controlling PC must also identify the parameters of the spatiotemporal AI layer from the measurement data obtained by the sensors in the IDS, determine on the basis of them the heat levels individual to the heating elements in the SDS, and send the level data to the corresponding elements. The first of the three tasks should be implemented in the parameter-identification software package (PISP) installed on the PC.

The main questions to be answered in the IDS design are the following:

・ How can one design an IDS in order to enable it, firstly, to be relevant to an SDS for the BSs of a wind turbine, secondly, to be reliable, cost- and energy-efficient, and, thirdly, to supply the prospective PISP with the data measured by the ISD sensors?

・ What are the parameter-identification procedures that can be implemented in the PISP?

The purpose of the present work is answering these questions by development of design issues of an IDS relevant to an SDS for the wind-turbine BSs and procedures for identification of the AI-layer parameters on the basis of acoustics of viscoelastic solids. The work is underlain by the previous works of the authors [16] and [1] . Each of them considers a thin curvilinear deformable solid layer associated with the BS and develops a procedure for identification of the material parameters of the layer. Paper [16] deals with a one-layer shell separating its interior and exterior, for instance, the non-QE deformable solid shells of various engineering systems (e.g., the shells of driver or passenger compartments of ships, cars, busses, airplanes, and other vehicles) or BS of a wind turbine. Paper [1] deals with the two-layer system of the BS layer and the AI layer accreted on the BS-layer outer surface.

Both works [16] and [1] consider the stress in a layer system as the one caused by the operational load on the system. The works describe the QE component of the average normal stress (ANS) with a third-order partial differential equation (PDE) of acoustic of viscoelastic solids. The deviatoric (or shear) stress is neglected for compactness of the model. The equation is introduced in [16] , takes into account the stress-relaxation function (in the exponential approximation) in the integrand of the Boltzmann super- position integral, includes the stress-relaxation time (SRT) of the layer material, and is relevant at any value of the SRT. The latter property enables application to viscoelastic materials, in particular, AI (see Remark 1.2).

The parameter-identification procedures in papers [16] and [1] presume structural- health/operational-load monitoring (SH/OL-M) (e.g., [17] ) with the help of the sensors located on the layer-system inner surface, i.e., the surface not affected by the environ- ment. However, these papers consider a use of accelerometers as the sensors. This restricts application of the procedures to the QE systems only. In contrast to that, the present work regards a use of the sensors, which enable application of the parameter- identification procedures to general, non-QE layer systems. Also, papers [16] and [1] consider homogeneous conditions at the layer-system outer surface, i.e., the one af- fected by the environment, but do not include the operational-load-caused source term in the above equation. The present work overcomes these limitations as well.

The answer to the first question in the above bullet list is developed in terms of methods of SH/OL-M and with the emphasis on monitoring stress in BSs in Section 2. The answer to the second question in the mentioned list is developed in Appendixes A-C. Appendix A deals with the identification of the material parameters of a one- layer system. Appendix B shows how this identification approach is generalized for multilayer systems. Appendix C exemplifies this generalization in the case of the BS/AI-layer system, which comprises two layers. Section 3 summarizes the obtained results and presents the concluding remarks. The section "Notations", which also includes the list of the used abbreviations, completes the work.

2. Design Issues of an IDS, Which Is Relevant to an SDS for the BSs of a Wind Turbi

Wind turbines are driven by irregular wind under irregular weather conditions, more specifically, by the air flows at the outer surface of the BS/AI-layer. Practical techniques for measurement of the wind effect are still unknown, even if the outer surface of the BS is free from AI. The absence of the techniques deprives some of the most common acoustic methods. For example, structural dynamics (e.g., [18] ) and lookup tables (e.g., [19] ) cannot be applied because the data on the wind effect on the BS/AI-layer, which are necessary parts of the input data for these approaches, cannot be measured.

The only practical techniques are the ones, which do not use measurements of the above wind effect. These techniques comprise methods of SH/OL-M. They are based on the data measured by sensors and used in works [1] and [16] . Methods of SH/OL-M apply various sensors, such as accelerometers or strain gauges (e.g., [20] ).

In a deformable solid system, an acoustic signal, which is measured at a spatial point of the system and can be regarded in SH/OL-M, is the QE component of one or another mechanical variable at the point. Until now, the settings for measurements of the QE components of the accelerations or strains in the non-QE systems are unknown.

The settings for measurements of the QE component of the stresses in the non-QE systems are unknown either. However, the remark below proposes the corresponding idea.

Remark 2.1. If a non-QE deformable solid system includes a cavity filled by air, which remains at equilibrium and at the atmospheric pressure, then the QE component of the stress in the wall of this cavity can be measured by an appropriate pressure/ normal-stress sensor as the difference between the values of the normal stresses in the solid wall and air in the cavity. □

An example of the non-QE systems noted in Remark 2.1 is the BS of an operating wind turbine. Indeed, the BS is a hollow solid body where the cavity is isolated from the external air. Other examples are the non-QE closed deformable solid hells and interiors of driver or passenger compartments of ships, cars, busses, airplanes, and other vehicles.

The sensors that can provide the measurements indicated in Remark 2.1 are the pressure-sensing resistors also known as the force-sensing resistors (FSRs). They are used in SH/OL-M. For example, in the wireless array of FSRs of the Honeywell Tech- nology Center ( [21] , Sections 3.2.3 and 3.2.4), the FSRs are identical. Each of them is a small-area thin planar resistor, which includes the film sensitive to the pressure or, more precisely, normal-stress difference at the opposite planar surfaces of the film, or to the related force.

The difference of the normal stresses at the opposite planar surfaces of an FSR, , and the corresponding force, , measured by the device, are coupled with relation

(2.1)

where is the time and is the sensing area of the FSR. One of many examples of FSRs is the FlexiForce A301 device [22] .

Remark 2.2. The operating-temperature interval of the FlexiForce Standard Model A301 FSR is from up to. The thickness of the device is. The width and total length of the device is and (approximately), respectively. This length includes 6-millimeter pins.

The sensing region of this device is a disk with the diameter. The area of this region, i.e., parameter in (2.1), is. The version “445 N” can measure force in the range from 0 to. Taking into account that, the mentioned area, and relation (2.1), one obtains that the interval of the measured pressure is from 0 to (i.e., somewhat greater than 6 MPa). The force reading change per degree of temperature change is 0.36. The output signal of the above sensor is presumed to be processed by the MCP6004 low-power operational amplifier of Microchip Technology. The operating-temperature range of this circuit includes the industrial temperature range, which is from to.

According to the sensor data sheet [23] , the power consumed by one sensor together with the related electrical circuit, which includes one MCP6004 amplifier, does not exceed. □

The operating-temperature intervals noted in Remark 2.2 include a considerable part of the temperature values discussed in Remark 1.1. The interval of the measured pressure is also rather wide. Is it possible to measure the QE component of the normal stress in the BS/AI-layer system with FSRs? In order to answer this question, one can take a closer look at key aspects of the aerodynamics of the operation of a wind-turbine rotor.

The main aerodynamic effect, by which the operating rotor extracts the energy, is the pressure drop between the air domains directly in front of, and directly behind, the BS of a rotor. More specifically, the picture is the following.

If there is no wind, the BS as well as the air domains at the inner and outer surfaces of the BS are at equilibrium, and their pressures are the same and equal to the at- mospheric pressure.

If there is a wind and the wind-turbine rotor operates, then:

・ the air domain at the outer surface, which is directly in front of the BS, is at the above-atmospheric pressure,

・ the air domain at the outer surface, which is directly behind the BS, is at the below- atmospheric pressure, and

・ the BS is not at equilibrium,

whereas the pressure at the inner surface of the BS remains atmospheric. Con- sequently, in the course of the operation, there is always the pressure difference between the BS inner surface not affected by a wind and the BS outer surface affected by a wind. This difference manifests the presence of the stress distributed along the thickness of the BS.

This stress can be measured as follows. Assume that the working planar surface of an FSR is attached to the inner surface of the BS. Consequently, the opposite planar surface of the FSR contacts the equilibrium atmospheric-pressure air in the interior of the BS. Then, according to the well-known continuity of the QE component of the stress at, and normal to, the interface between two solids (e.g., [24] , (1.48)), the normal stress at the FSR/BS interface is equal to the QE component of the stress at the interface. Similarly, according to the well-known continuity of the QE component of the stress at, and normal to, the interface between a solid and air, the QE component of the stress at, and normal to, the opposite planar surface of the FSR, i.e., the FSR/air interface, is equal to zero (e.g., [24] , (1.49)). Since the FSR measures the normal-stress difference between its surfaces (see the text on (2.1)), it in fact measures the QE component of the mentioned stress in the BS at the FSR/BS interface.

The above part of the present section can be summarized as follows.

A sensor, which can measure the QE component of the normal stress in a non-QE closed deformable solid shell such as a BS of a wind turbine, is an FSR with its working planar surface attached to the BS inner surface and its opposite planar surface con- tacting the equilibrium air in the BS cavity. An array of the FSRs can be wirelessly connected to, and controlled by, a PC.

This picture is equally applicable to a more general case where the layer system measured by FSRs includes not only the BS layer but also the AI layer that can be accreted on the BS outer surface. The corresponding parameter-identification pro- cedures are developed in Appendixes A-C.

The above settings can be implemented in a hardware configuration that includes an electronic communication subsystem for the sensor array. The components of it are standard manufactured (or off-the-shelf) products, very small and light, reliable, high- speed, inexpensive, and providing flexible scaling. Importantly, the sensor array and communication subsystem does not presume parts placed outside the BS.

To be specific, the present section considers a wind-turbine rotor that has three identical blades. The blade length is assumed to be of 50 (i.e., similar to the one of a fiberglass-reinforced epoxy blade of the Siemens SWT-2.3-101 wind turbine [25] ). As one can estimate, the area of the BS layer of this blade is about 400. In order to sense the AI landscape discussed in Section 1, one should, near but outside the area of each heating element of an SDS (see Section 1), attach the FSRs, one per, say, four square meters. Thus, the area of 400 is covered with 100 sensors (see Row 1 of Table 1). As is noted above in this section, all of the sensors are attached to the inner surface of the BS. They can be connected by a flat cable of the Serial Peripheral Interface (SPI) bus (e.g., [26] ) or the Inter-Integrated Circuit (I^{2}C) bus (e.g., [27] ). More specifically, the con- figuration for a rotor blade is the following.

Table 1. Key components of the electronic equipment for a three-blade rotor of a wind turbine according to the present approach. The communication hardware comprises the components in Rows 2 and 4 - 6.

* if one orders at the quantity of 300.

First, the bus cable of the length of about 210 is attached to the inner surface of the BS (see Row 2 of Table 1) in the form of a meander.

Next, 100 sensors are rigidly attached (stuck) to the same surface along, and in a close proximity to, the cable at the distances of about 2 between any two neighboring sensors. Note that the cable and sensors are placed in such a way that each sensor corresponds to the surface area of about four square meters. Then each of the sensors is connected to the cable via an amplifying circuit (e.g., see [23] for the recommended circuit, which includes the amplifier) and the micro-controller (see Rows 3 and 4 of Table 1).

Finally, each of the three cables is carried out from a respective blade into the rotor hub and connected to a micro-controller, a GSM/GPRS communication module, and a rechargable long-life battery that are rigidly attached to the solid body of the hub (see Rows 5 and 6 of Table 1). The macro-controller can be of the same type as the one indicated in Row 4 of Table 1. The battery is recharged from the energy produced by the turbine.

The GSM/GPRS module wirelessly transmits the measured data to a remote PC, which is also equipped with a GSM/GPRS module, processes the data for identification of the material parameters of the AI that can be accreted on the BS, and controls the sensor array. The corresponding estimated cost for the described hardware does not exceed 31 kSEK per one wind turbine with three 50-meter blades (see the last row and column of Table 1). As follows from the last sentence in Remark 2.2, the power consumed by the above electronic system necessary for the operation of 300 sensors does not exceed 6.

The outcomes of the present section is the answer to the question in the first bullet in Section 1. The answer to the question in the second bullet is developed in Appendixes A-C.

3. Obtained Results and Concluding Remarks

Summing up the present work, one can note the following.

The present work is based on the third-order partial differential equation (PDE) of acoustics of viscoelastic solids for the quasi-equilibrium (QE) component of the average normal stress derived and used in the previous papers of the authors. This PDE includes the stress-relaxation time (SRT) for the material and is applicable at any value of the SRT.

The work specifies the notion of a smart deicing system (SDS) for blade shells (BSs) of a wind turbine. The stress in a BS is considered as the one caused by the operational load on the BS. The work developed key design issue of a prospective ice-detection system (IDS) able to supply an array of the heating elements of an SDS with the element-individual spatiotemporal data (see Section 2) and procedures for identifi- cation of the material parameters of atmospheric ice (AI) layer accreted on the outer surfaces of the BSs (see Appendixes A-C). Both the SDS and IDS flexibly allow for complex, curvilinear and space-time-varying shapes of BSs.

The proposed IDS presumes monitoring of the QE components of the normal stresses in BSs. The IDS is supposed to include an array of force-sensing resistors (FSRs) and communication hardware, as well as the parameter-identification software package (PISP), which provides the identification on the basis of the aforementioned PDE and the data measured by the FSRs. The IDS does not have hardware components located outside the outer surfaces of, or implanted in, the BSs. The FSR array and communication hardware are:

・ Reliable because they comprise standard manufactured (or off-the-shelf) products only, which can, moreover, operate at the temperatures between −40˚C and +40˚C;

・ Cost efficient because their estimated cost is 30 - 35 kSEK for a rotor with three 50-meter blades;

・ Energy efficient because their estimated power consumption is within 10 watts in the case of the above rotor.

The present work extends methods of structural-health/operational-load monitoring (SH/OL-M) with measurements of the op- erational-load-caused stress in closed solid shells and, if the prospective PISP is used, endows the methods with identification of material parameters of the shells. The identification algorithms that can underlie the PISP are computationally efficient and suitable for implementation in the real-time mode.

The identification model and algorithms can deal with not only the single-layer systems such as the BS layer without the AI layer (see Appendix A as well as Table A1 for the input data and the parameters that can be identified) or two-layer systems such as the BS with the AI layer accreted on it (see Appendix C as well as Table C1 for the input data and the parameters that can be identified) but also multi-layer systems (see Appendix B). The outcomes are applicable to not only the BSs of wind turbines but also the non-QE closed single- or multi-layer deformable solid shells of various engineering systems (e.g., the shells of driver or passenger compartments of ships, cars, busses, airplanes, and other vehicles). The proposed monitoring of the normal-stress QE component in the mentioned shells extends methods of SH/OL-M.

The outcomes of the present work complement and further develop the results of the previous works of the authors, more specifically, papers [1] , [16] , and [32] . The topic for the nearest research is a better adjustment of the settings for the FSR-based mea- surement of the normal-stress QE components in BSs and a calibration (e.g., see Remark A.2) of the parameter-identification model and algorithms, as well as the resulting improvement of the PISP.

Acknowledgements

The authors express their gratitude to the Swedish Energy Agency for a partial support of the present work via the project 37286-1 in the “Wind power in cold climates” program. The authors also thank Andrey Koptyug, Sports Tech Research Centre, Department of Quality Technology, Mechanical Engineering and Mathematics, Mid Sweden University, Östersund, Sweden, for drawing the attention to possible use of force- sensing resistors and suggestion on the hardware configuration described in Section 2.

Appendixes

This section comprises Appendixes A-C. They consider the stress in a layer system as the one caused by the operational load on the system and present necessary details on the models and methods for identification of the material parameters of the BS and AI layers.

A. Identification of the Material Parameters of the BS of an Operating Wind Turbine

The present section considers the case where the AI layer is not present at the outer surface of the BS. The QE component of ANS, , can be described with linear PDE

, (A.1)

which is the non-homogeneous generalization of partial differential equation (PDE) ( [16] , (2.11)). In (A.1), x, y, and z are the spatial coordinates and term is due to the operational load on the BS (in particular, associated with the position and wind- induced rotation of the distributed mass of the BS). Consequently, in general depends on x, y, z, and t. Note that, at equilibrium, need not be zero (e.g., [14] ).

The meanings of two other terms in (A.1) are the following: is the speed of the bulk acoustic waves (i.e., the ones that correspond to compressions/rarefactions in the shear-free case) in the BS and is the stress-relaxation time in the BS. Note that inequalities

(A.2)

(A.3)

hold due to the physical meanings of the two parameters. Speed is determined as follows (e.g., [16] , (2.7))

(A.4)

where is the bulk modulus and is the volumetric mass density. If both and are available, then the volume viscosity of the BS, , can be determined as follows (cp., the first equality in (1.5)).

There are other contributions of the operational load. For example,

(A.5)

(A.6)

Remark A.1. Parameters and are the material parameters of the BS to be identified. If an analytical or tabular dependence of the speed of the bulk acoustic waves in the BS on the volumetric mass density of it, say, is available, then is determined as the solution of equation

(A.7)

This enables one to evaluate from (A.4). Thus, and can also be available. □

The material parameters of the BS to be identified, and, are parameters of Equation (A.1). Usually, parameters of an equation can be identified if a solution of the equation is available. However, Equation (A.1) includes term, which is unknown. Moreover, the pressure and acceleration noted in features (A.5) and (A.6) are also unknown. Consequently, description (A.1), (A.5), (A.6) seems useless for the parameter identification.

However, before rejecting it, one can consider particular, more specific cases of this description in hope that a better specificity will change the role of the unknown terms in such a way that it will be possible to include the “successors” of them in the identifi- cation procedure. For the reason explained in Section 2, one can focus on the distri- bution of along the thickness of the BS, say, the x-axis. Note that

(A.8)

In order to pass to the x-version of Equation (A.1), we first rewrite it in the following equivalent form

(A.9)

where

, (A.10)

and then substitute equalities (A.8) into (A.9). This results in

, (A.11)

where

(A.12)

(A.13)

Without a loss of generality, one can assume that the x-axis is normal to the BS layer and point corresponds to the inner surface of the BS. Common relation for the QE component of the ANS (A.12) at the solid-air interface is (e.g., [24] , (1.49))

(A.14)

The inner air is located to the left from this interface. As noted in Section 2, it remains at the atmospheric pressure,.

Let be the thickness of the BS along the x-axis. Then point corresponds to the outer surface of the BS. One can also assume that the BS layer is spatially regular sufficiently in order to allow to regard the x-axis as normal to the outer BS surface at point .

For the sake of simplicity, we consider the air acceleration indicated in (A.6) to be identically zero. Then, the effects of the phenomena (A.5) and (A.6) at point can be formulated as follows

(A.15)

(A.16)

where the outer-air pressure along the x-axis at the BS/outer-air interface. The obtained description for the QE component of the ANS (A.12) comprises (A.11) and (A.13)-(A.16).

It is in one spatial coordinate and is formally simpler than description (A.1), (A.5), (A.6) in three spatial coordinates. Still, the new settings do not enable one to obtain solutions of equation (A.11) because term (see (A.13) and (A.10)) in it and term in (A.15) are unknown. However, data of a sensor can help to specify at least some terms in Equation (A.11) or its particular version.

Indeed, if one, at point, separates the BS from the inner air with a body of an FSR (see Section 2), then the interface at point is the one between two solids (rather than a solid and air) that motivates to change relation (A.14) to (e.g., see [24] , (1.48))

(A.17)

which expresses the continuity of the stresses at the mentioned interface. The passing from (A.14) to (A.17) make it possible to obtain a prototype for an ordinary differential equation (ODE), which can be used for the identification. Indeed, applying value to Equation (A.11), one obtains

(A.18)

that, after substitution of (A.17), takes the following form

(A.19)

where the derivatives on the left-hand side are known because values of function are available as the FSR data. Relation (A.19) is an ODE with the known solution. This ODE can be specified further if the spatial derivative on the right-hand side is expressed in terms of the solution.

The simplest way to this expression is based on an approximation for, as function of x, in the form of a polynomial in x with the t -dependent coefficients that can be determined from the information on (cp., the method of ( [16] , Section 4]). This information comprises three identities, (A.15)-(A.17). Consequently, the polynomial has three coefficients. We consider its following specific form

(A.20)

where the t-independent integer number

(A.21)

is the degree of the polynomial. Representation (A.20) provides the estimation for the aforementioned spatial derivative, namely

. (A.22)

Substitution of (A.20) into (A.15)-(A.17) results in

(A.23)

(A.24)

(A.25)

Combination of (A.22) and (A.24) transforms (A.19) into ODE

(A.26)

for function where

(A.27)

(A.28)

Note that, by virtue of (A.2) and (A.21), inequality

(A.29)

holds. We also note that, in view of the role, which parameter

(A.30)

plays in ODE (A.26) for value (A.17) of a solution of PDE (A.11), this parameter can be interpreted as the characteristic wave number for the mentioned value.

Equation (A.26) presents an ODE with:

・ Solution available as the measurement results provided by the FSR,

・ Unknown parameters and, and

・ Function, which is unknown due to the fact that term in (A.27) and term in (A.10), and, thus, term in (A.27) are unknown.

The latter feature indicates that function should, in one or another way, be included in the identification of parameters and or, by virtue of (A.28),. This can be implemented by the procedure, which is described below and starts with the following three auxiliary steps. The corresponding input data and parameters that can be identified are listed in Table A1.

One considers ODE (A.26) at any three successive time points, say, , , and measured by the FSR (see the upper half of Table A1)

(A.31)

where, , and are the finite-difference (FD) approximations for the first, second, and third time derivatives of function at the mentioned points calculated on the basis on of the data measured by the FSR (see the upper half of Table A1). We assume that the length of the time interval between the left and right time points, i.e., , is much smaller than the characteristic time of function in

Table A1. Identification method for one-layer systems: The input data and the parameters that can be identified.

this interval. This feature allows to assume that. Denoting this value with and substituting it into (A.31), one gets

(A.32)

Equations (A.32) present a system of bilinear equations with constant coefficients for unknown numbers, , and. The left-hand sides of equalities (A.32) present the linear part of this system. Solution of the system can be simplified if one expresses, , and on the left-hand sides in terms of the right-hand sides of the system. Assuming that matrix of the linear part is nonsingular, i.e.,

, (A.33)

one easily calculates the and coefficients in the resulting expressions

(A.34)

(A.35)

(A.36)

If the inequality in relation (A.33) is not valid, then the rank of the matrix in this relation is less than three. Consideration of this special case is not difficult. However, it is beyond the scope of the present work.

Solution of system (A.34)-(A.36) with respect to, , and is rather simple in the following three cases.

If and, then, , and are the identified values of the three parameters provided that inequalities (A.29) and (A.3) hold.

If and, then, , and are the identified values of the three parameters provided that inequalities (A.29) and (A.3) hold.

If and, then, , and are the identified values of the three parameters provided that inequalities (A.29) and (A.3) hold.

The case where and is more complex. In this case, we multiply (A.34) and (A.35) by and, respectively, and subtract the second of the resulting equalities from the first one. This results in relation or, equivalently,

(A.37)

Substitution of (A.37) into (A.35) transforms the latter into quadratic equation

(A.38)

If this equation has exactly one root, which meets condition (A.3), and the cor- responding value of indicated by (A.37) meets condition (A.29), then the de- termined values are the identified parameters. They determine the identified value of by means of (A.36).

Remark A.2. In each of the above cases, parameter is determined from (A.28) by means of the identified and parameter n (see (A.21)). This means that n must be a part of the input data and, thus, be available before the identification. It can be determined from the calibrating identification, i.e., the one at an already available (see (A.4)). In this case, n is calculated from (A.28) by means of the already known and identified.

Thus, the input data cannot be obtained without measured data for stress. □

Remark A.3. Importantly, the above identification of parameters and also provides identification of parameter (see the text above (A.32)). It is the corre- sponding value of unknown function (A.27), which includes unknown, wind-related terms and (A.10). Due to this, the identified source term in (A.32) pre- sents valuable information.

The above also shows that the proposed method is applicable in spite of the presence of the unknown terms in the model. This is an important practical advantage of the method. □

As soon as is identified, parameter is identified according to (A.28). The rest of the identification procedure is described in the last sentence of Remark A.1.

Remark A.4 (cp., the discussion in ( [16] , Section 4). The proposed method identifies the parameters in a time interval comprising three successive time points, thereby presuming that the parameters are independent of time in this interval. The method can be applied to each of the three-successive-point intervals that can be considered in the time-point sequence indicated as a prat of the input data for the method (see the upper half of Table A1). As a result, one obtains each of the identified parameters in the form of a piecewise constant function of the time. This function need is not single- valued because the three-point intervals are mutually intersecting and, at the intersections, the function can have two values. The time dependences of this type are rather irregular and need special techniques for smoothing or other processing. These techniques are beyond the scope of the present work. □

The next section outlines how the proposed approach can be generalized for multi- layer systems.

B. Generalization of the Parameter Identification from One-Layer Systems to Multi-Layer Systems

The previous section proposes a method for the identification of the material para- meters of the BS layer of an operating wind-turbine rotor. The present section explains how this approach is generalized for the identification of the parameters of multi-layer system.

One can consider the l-layer system where, the x-axis is directed along the normal to all of the layers, point, as before, corresponds to the inner surface of the system, and, , is the thickness of the j th layer. Then and, , are the coordinates of the boundaries of the layers. Note that is the outer surface of the layer system. Thus, the l layers correspond to l intervals: and,.

Equations (A.18) or (A.19) corresponds to the left boundary of the x-interval considered in the one-layer treatment. The counterparts of (A.19), i.e., the equations corresponding to the left boundaries of the above l intervals are

(B.1)

(B.2)

The boundary data at the layer surfaces are

(B.3)

(B.4)

(B.5)

(B.6)

(B.7)

Relation (B.3) is similar to (A.17) because it is set at the inner surface of the layer system. Relations (B.6) and (B.7) are similar to (A.15) and (A.16) because they are set at the outer surface of the layer system. Relations (B.4) and (B.5) express the continuity of the normal stress and normal acceleration at the interfaces between the layers.

Remark B.1. The number of the boundary equalities (B.3)-(B.7) is. □

Terms and, , are the boundary values included in (B.3) and (B.4). Relation (B.3) determines in terms of normal stress measured by the FSR. This enables one to rewrite equation (B.1) as

, (B.8)

which is similar to (A.19). One can also consider a possibility to express boundary values, , in terms of. Moreover, following the idea of the method presented in the previous section (see the text below (A.19)), one needs to express the spatial derivatives on the right-hand sides of (B.8) and (B.2) in terms of.

Proceeding in these directions, one can involve a polynomial similar to (A.20) for each layer. This results in polynomials

(B.9)

(B.10)

where n and are independent of j, and n is described as before (see the text on (A.21)). They in particular provide estimations (cp., (A.22))

, (B.11)

(B.12)

The number of the coefficients of polynomials (B.9) and (B.10) is, i.e., the same as the number of the boundary equalities (see Remark B.1). Consequently, these coefficients can be determined. The resulting expressions enable one to express, , , and spatial derivatives (B.12), (B.13) in terms of. This, in turn, allows obtaining the corresponding versions of ODE (A.26) and implementing the analysis analogous to the one described in the previous section.

The above generalization for multi-layer systems is exemplified with application to a two-layer system, which comprises the BS and AI layers, in the next section.

C. Identification of the Parameters of the AI Layer Accreted on the BS of an Operating Wind Turbine

Section A considers the case where the x -interval for the BS is. The present section deals with a more general case where the AI layer of the thickness h is accreted on the outer surface of the BS. Consequently, the x-interval for the AI layer is. Thus, the layer system under consideration comprises two layers.

We specify the two-layer version of equations (B.8), (B.2), boundary data (B.3)-(B.7), polynomials (B.9), (B.10), and terms (B.11), (B.12) in the following forms

, (C.1)

, (C.2)

, (C.3)

, (C.4)

, (C.5)

, (C.6)

, (C.7)

(C.8)

(C.9)

, (C.10)

(C.11)

where P is the QE component of the ANS in the AI and the coefficients without the subscript “s” corresponds to the AI. Note that, in view of (C.9),

(C.12)

One can determine the five coefficients of polynomials (C.8) and (C.9) from the five boundary relations (C.3)-(C.7). The results are:

(C.13)

(C.14)

(C.15)

(C.16)

(C.17)

where

(C.18)

Remark C.1. In the limit case as, i.e., in the limit where the AI layer is not present, relation also holds (see (C.18)), and values (C.13), (C.14), and (C.17) tend to (A.23)-(A.25), respectively, as must be. Moreover, in this limit case, tends to due to (C.9), (C.18) and (C.15). □

Application of (C.14) and (C.16) to (C.10) and (C.11), respectively, leads to

(C.19)

(C.20)

Substituting (C.19) and (C.20) into the right-hand sides of (C.1) and (C.2), re- spectively, one obtains

(C.21)

(C.22)

where is described with (A.28), parameters and in (C.21) are a part of the input data (see the upper half of Table C1), the role of n is the same as the one described in Remark A.2, and

(C.23)

(C.24)

(C.25)

Table C1. Identification method for two-layer systems: The input data and the parameters that can be identified.

(C.26)

In view of (C.23) and (C.18), the acceptable values of v are such that

(C.27)

Equalities (C.21) and (C.22) present two different ODEs for function. This function is a part of the input data (see the upper half of Table C1). The source terms (C.25) and (C.26) are unknown functions because they include unknown terms:, , and. Thus, ODEs (C.21), (C.22) and source terms (C.25), (C.26) are the two-layer-system counterparts of the one-layer-system (A.26) and source term (A.27).

The ODEs can be used for identification of the AI layer parameters q and provided that (and, thus (see (A.28)), , and are available (see the upper half of Table C1). This can be implemented with the help of equation (C.22). In order to enable that, the procedure should include identification of the corresponding source term. The value of u in (C.22) can be determined from (C.18) by means of identification of the two parameters, v and the corresponding source term in (C.21). Thus, the three parameters, q, , and the source term in (C.22), can be determined after the above two parameters are determined.

The corresponding procedure can follow the already familiar way (see the text on (A.31) and (A.32)). More specifically, one considers ODEs (C.21), (C.22) in the interval comprising any three successive time points, say, , , and measured by the FSR (see the upper half of Table C1) where any two of the three points, say, and, are used for ODE (C.21) to identify the related two parameters, and all of the three points are used for ODE (C.22) to identify the other three parameters. This leads to the following versions of the ODEs

(C.28)

(C.29)

where, , , and are the FD approximations for the first, second, and third time derivatives of function at the mentioned points calculated on the basis on of the data measured by the FSR (see the upper half of Table C1). We assume that the length of the time interval between the left and right time points, i.e., is much smaller than the characteristic times of functions and in this interval. This feature allows to assume that and . Denoting these values with and, respectively, and substituting them into (C.28) and (C.29), one gets

(C.30)

(C.31)

By solving system (C.30) of two linear algebraic equations, one determines v and where v should meet requirement (C.27). Then, u and h are evaluated according to (C.23) and (C.18),

(C.32)

(C.33)

If, then (see (C.33)) as well that means that the AI-layer thickness is zero, i.e., this layer is not present. In this case, there are no material parameters of the layer to be identified.

If, then (see (C.32) and (C.27)). In this case, one substitutes u into equation system (C.31) and identifies parameters, , and from it. Since system (C.31) is exactly of the same form as the one of system (A.32), the parameters can be identified with the same method, i.e., the one described in the text on (A.33)- (A.38). The first two of the parameters should meet conditions (cp., (A.29), (A.3))

(C.34)

(C.35)

As soon as q is available, parameter s is evaluated from (C.24) (where h is determined with (C.33)). Then (cp., Remark A.1), is calculated as the unique positive solution of equation

(C.36)

(see the upper half of Table C1 for function). The obtained value of provides the bulk modulus and porosity of the AI, and, by means of (cp., (A.4)) and (1.4), respectively.

The proposed identification method is a generalization of the one of Appendix A for two-layer systems. However, the discussion on the one-layer method, in particular, Remarks A.3 and A.4 (along with Remark A.2 already mention in the present section) are equally applicable to the present, two-layer case.

Notations

Abbreviations

AI―atmospheric ice

ANS―average normal stress

BS―blade shell

FD―finite difference

FSR―force-sensing resistor

IDS―ice-detection system

ODE―ordinary differential equation

PC―personal computer

PDE―partial differential equation

PISP―parameter-identification software package

QE―quasi-equilibrium

SDS―smart deicing system

SH/OL-M―structural-health/operational-load monitoring

SRT―stress-relaxation time

Roman Uppercase Letters

―force corresponding to (see (2.1))

―shear modulus of AI

―bulk modulus of AI

―bulk modulus of a BS

―counterpart of term in the case of the AI layer in the BS/AI-layer system

―term for a BS determined with (A.10)

―version of term in the case of the jth layer of a multi-layer system

―term in (A.1) due to the position and wind-induced rotation of the distributed

mass of a BS

―number of time points

―version of quantity in the case of the AI layer in the BS/AI-layer system

―version of quantity in the case of the jth layer in a multi-layer system

―QE component of the ANS in a BS

―term for the AI layer in the BS/AI-layer system determined with (C.26)

―term for the BS determined with (A.27) in Appendix A and with (C.25) in

Appendix C

―approximate value of function determined as is described in the text between (A.31) and (A.32)

―sensing area of an FSR

Roman Lowercase Letters

―version of coefficient in the case of the AI in the BS/AI-layer system

―version of coefficient in the case of the jth layer in a multi-layer system

―coefficient of polynomial (A.20) for the BS layer

―coefficient of polynomial (A.20) for the BS layer

―version of coefficient in the case of a multi-layer system

version of coefficient in the case of the AI layer in the BS/AI-layer system

―version of coefficient in the case of the jth layer in a multi-layer system

―coefficient of polynomial (A.20) for a BS

―thickness of the AI layer

―thickness of the BS layer

―integer index,

―integer index,

―number of the layers in a multi-layer system

―degree of polynomial (A.20) or polynomials (B.9), (B.10)

―parameter determined with (C.24)

―parameter determined with (A.28)

―speed of the bulk acoustic waves in AI

―version of parameter in the case of the jth layer of a multi-layer system

―speed of the bulk acoustic waves in a BS

―analytical or tabular dependence of the speed of the bulk acoustic waves in

AI on

―analytical or tabular dependence of the speed of the bulk acoustic waves in a BS on

―time

―successive time points, at which an FSR measures the normal stress at the inner surface of a BS

―parameter determined with (C.18)

―three scalar spatial coordinates

―thickness of the first layers in a multi-layer system (see the second paragraph of Section B)

Greek Uppercase Letters

―difference of the normal stresses at the opposite planar surfaces of an FSR; stress in a BS at the location of an FSR on the inner surface of the BS

―values of the stress measured by an FSR at time points

―FD approximations for the first, second, and third time derivatives of function at time point calculated on the basis on of values and

―outer-air pressure along the x-axis at the interface between a BS and the outer air or at the interface between the last layer in a multi-layer system (e.g., the AI layer in the BS/AI-layer) and the outer air (in the AI-layer case, this interface is at)

Greek Lowercase Letters

―coefficients described in the text on (A.34)-(A.36)

―coefficients described in the text on (A.34)-(A.36)

―porosity of AI

―volume (or compressional) viscosity of AI

―volume (or compressional) viscosity of a BS

―stress-relaxation time in AI

―version of parameter in the case of the jth layer of a multi-layer system

―stress-relaxation time in a BS

―parameter determined with (A.30)

―shear viscosity of AI

―parameter determined with (C.23)

―volumetric mass density of AI

―volumetric mass density of air

―volumetric mass density of a continuous, non-porous ice

―version of parameters in the case of the jth layer of a multi-layer system

―volumetric mass density of a BS

Submit or recommend next manuscript to SCIRP and we will provide best service for you:

Accepting pre-submission inquiries through Email, Facebook, LinkedIn, Twitter, etc.

A wide selection of journals (inclusive of 9 subjects, more than 200 journals)

Providing 24-hour high-quality service

User-friendly online submission system

Fair and swift peer-review system

Efficient typesetting and proofreading procedure

Display of the result of downloads and visits, as well as the number of cited articles

Maximum dissemination of your research work

Submit your manuscript at: http://papersubmission.scirp.org/

Or contact jamp@scirp.org

Conflicts of Interest

The authors declare no conflicts of interest.

[1] | Mamontov, E. and Berbyuk, V. (2015) Passive Acoustic Signal Sensing Approach to Detection of Ice on the Rotor Blades of Wind Turbines. 16th International of Workshop on Atmospheric Icing of Structures (IWAIS), Uppsala, 28 June-3 July 2015, 75/248-80/248. |

[2] |
Sommerfeld, R.A. (1982) A Review of Snow Acoustics. Reviews of Geophysics and Space Physics, 20, 62-66. http://dx.doi.org/10.1029/RG020i001p00062 |

[3] | Gibson, R.F. (2012) Principles of Composite Material Mechanics. CRC Press, Boca Raton. |

[4] | Landau, L.D. and Lifshitz, E.M. (1986) Theory of Elasticity. Pergamon Press, Oxford. |

[5] | Gross, B. (1953) Mathematical Structure of the Theories of Viscoelasticity. Hermann & Cie, Paris. |

[6] | Junisbekov, T.M., Kestelman, V.N. and Malinin, N.I. (2003) Stress Relaxation in Viscoelastic Materials. Science Publishers, Enfield. |

[7] |
Hess, H. (1902) Elasticität und innere Reibung des Eises [Elasticity and Viscosity of Ice]. Annalen der Physik, 313, 405-431. http://dx.doi.org/10.1002/andp.19023130612 |

[8] |
Weinberg, B. (1905) über die innere Rei-bung des Eises [Viscosity of Ice]. Annalen der Physik, 323, 81-91. http://dx.doi.org/10.1002/andp.19053231103 |

[9] |
Weinberg, B. (1907) über die innere Reibung des Eises. II. Annalen der Physik, 327, 321-332. http://dx.doi.org/10.1002/andp.19073270208 |

[10] |
Deeley, R.M. (1908) The Viscosity of Ice. Proceedings of the Royal Society of London Series A, 81, 250-259. http://dx.doi.org/10.1098/rspa.1908.0077 |

[11] | Kobeko, P.P., Shishkin, N.I., Marei, F.I. and Ivanova, N.S. (1946) Plastic Deformation and Viscosity of Ice. Journal of Technical Physics, 16, 263-272. |

[12] | Hobbs, P.V. (1974) Ice Physics. Clarendon Press, Oxford. |

[13] |
Gao, H. and Rose, J.L. (2009) Ice Detection and Classification on an Aircraft Wing with Ultrasonic Shear Horizontal Guided Waves. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 56, 334-344. http://dx.doi.org/10.1109/TUFFC.2009.1042 |

[14] |
White, J.R., Adams, D.E. and Rumsey, M.A. (2009) Operational Load Estimation of a Smart Wind Turbine Rotor Blade. Proceedings of SPIE, 7295, 72952D/1-72952D/12. http://dx.doi.org/10.1117/12.815802 |

[15] |
Smith, B.A. and Maheri, A. (2012) Optimisation of the Efficiency of Carbon Fibre Heating Elements Implanted in Wind Turbine Blades. 2nd International Symposium on Environment-Friendly Energies and Applications (EFEA), Newcastle upon Tyne, 25-27 June 2012, 410-414. http://dx.doi.org/10.1109/efea.2012.6294045 |

[16] | Mamontov, E. and Berbyuk, V. (2015) Identification of Material Parameters of Thin Curvilinear Viscoelastic Solid Layers in Ships and Ocean Structures by Sensing the Bulk Acoustic Signals. 6th International Conference on Computational Methods in Marine Engineering, Rome, 15-17 June 2015, 502-513. |

[17] | https://en.wikipedia.org/wiki/Structural_health_monitoring https://en.wikipedia.org/wiki/Operational_loads_monitoring |

[18] | https://en.wikipedia.org/wiki/Structural_dynamics |

[19] | https://en.wikipedia.org/wiki/Lookup_table |

[20] |
Schroeder, K., Ecke, W., Apitz, J., Lembke, E. and Lenschow, G. (2006) A Fibre Bragg Grating Sensor System Monitors Operational Load in a Wind Turbine Rotor Blade. Measurement Science & Technology, 17, 1167-1172. http://dx.doi.org/10.1088/0957-0233/17/5/S39 |

[21] |
Schoess, J.N. (2001) Conduc-tive Polymer Sensor Arrays: A New Approach for Structural Health Monitoring. Proceedings of SPIE, 4335, 9-19. http://dx.doi.org/10.1117/12.434164 |

[22] | https://www.tekscan.com/products-solutions/force-sensors/a301 |

[23] | “The FLX-A301-A.pdf” File Can Be Downloaded from the Web-Site [22]. |

[24] | Pollard, H.F. (1077) Sound Waves in Solids. Pion, Lon-don. |

[25] | https://en.wikipedia.org/wiki/Wind_turbine_design#Blade_materials |

[26] | https://en.wikipe-dia.org/wiki/Serial_Peripheral_Interface_Bus |

[27] | https://en.wikipedia.org/wiki/I%C2%B2C |

[28] |
CA Mätsystem AB (2016) Täby, Sweden. www.camatsystem.com |

[29] | http://se.farnell.com/ |

[30] | https://www.elfa.se |

[31] | http://simcomm2m.com/En/module/?type=19 |

[32] |
Mamontov, E. and Berbyuk, V. (2014) A Scalar Acoustic Equation for gases, Liquids, and Solids, Including Viscoelastic Media. Journal of Applied Mathematics and Physics, 2, 960-970. http://dx.doi.org/10.4236/jamp.2014.210109 |

Copyright © 2020 by authors and Scientific Research Publishing Inc.

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.