Inference of thermal models for sensors

INFERENCE OF THERMAL MODELS FOR SENSORS A Master's Thesis Submitted to the Faculty of the Escola Tècnica d'Enginyeria de Telecomunicació de Barcelona Universitat Politècnica de Catalunya by Santiago Novio Vázquez In partial fullment of the requirements for the degree of MASTER IN ELECTRONIC ENGINEERING Autor: Santiago Novio Vázquez, Advisor: Manuel Domínguez Pumar. May 20, 2016 1 Abstract The presence of thermal diusivity in a spherical thermal anemometer gives it long-memory dependence. The identication problem and state realisation of this model is addressed by using diusive representation (DR). In order to do that, an ideal thermal sphere quadrupole model and its corresponding nite dierences model are proposed and simulated. From those models and the non-rational Cole-Cole transfer function, the identication problem is discussed. PRBS as an input signal is also discussed. The best choice in the number of poles and their position is found by analysing frequency response. For the 5 decades bandwidth used here, it has been found that 8 poles geometrically spaced by a scale factor of 5 is a good choice. Finally, 8-poles DR is veried as a good option to model a complex spherical thermal anemometer prototype which is being developed in Universitat Politica de Catalunya, by means of open loop and sigma-delta closed loop control simulations. 2 A Carme, a quien quiero mucho, por haber compartido nuestra vida durante el tiempo que he hecho este trabajo, por haber entendio que dedicase parte de ella a este trabajo, y por querer seguir compartiendola, A mis padres, por sus consejos, por estar siempre ahí aunque no estemos cerca, y por ver bien todo lo que hago por muy loco que sea. A mis hermanos, por ser como mis padres y formar parte de mis locuras. A todos ellos y al resto de gente que conozco, por aguantarme y porque ahora que tendré mas tiempo, van a tener que aguantarme más aún. 3 Acknowledgements I would like to thank Manuel Domínguez Pumar for putting his trust in me from the very beginning, for oering me his hand and constantly proposing me new challenges, for his patience to explain some of the concepts underlying this work, which were new for me and sometimes dicult to understand. I would also like to thank him for being like this, although due to personal circumstances, my work was intermittent and I did not get in touch with him for many days, and for giving me the freedom to set my own course in the work and giving me time to really understand what I was doing and to do things which were not so relevant. I would like to thank Teresa Atienza for sharing and contrasting the results of this work, which initially we did in parallel. I would also like to thank her for providing me with the data les of simulations in the laboratory, so that I could make the graphs included in the last part of this work and, therefore, validate the method with experimental results. I would like to thank Lukasz Kowalski for his predisposition to help at all times and for his information on the laboratory experiments. I would also like to thank him for designing a simulation of the ideal sphere in COMSOL, which in the end was not included in the work because that simulation had a computational cost which was too high, so it could not be done. I would like to thank all the sta in the Universitat Politècnica de Catalunya and all my mates for their welcome and support while this Master was being carried out. Thanks. 4 Revision History and Approval Record Written by: Reviewed and approved by: Date 19/05/2016 Date 20/05/2016 Name Santiago Novio Vázquez Name Manuel Domínguez Pumar Position Project Author Position Project Supervisor 5 Contents Abstract 2 Acknowledgements 4 Revision History and Approval Record 5 Contents 6 List of Figures 8 List of Tables 13 1 Introduction 14 1.1 Stateoftheart............................ 14 1.2 Methodology ............................. 14 2 Diusive Representation General Formulation 15 2.1 Finite-Dimensional Approximations . . . . . . . . . . . . . . . . . 16 2.2 Identication of the Diusive Symbol in Time Domain . . . . . . 18 3 Numerical Examples 20 3.1 Ideal Sphere Thermal Quadrupole Model . . . . . . . . . . . . . 20 3.1.1 Theoretical Diusive Symbol of the ISTQM thermal impedance 23 3.1.1.1 Simulation of ISTQM Theoretical Diusive Symbol ......................... 24 3.2 Diusive Representation of Non-Rational Cole-Cole Transfer Function .................................. 28 3.2.1 Choosing Input Signal . . . . . . . . . . . . . . . . . . . . 29 3.2.2 Choosing Poles and Bandwidth . . . . . . . . . . . . . . . 31 3.2.3 Other Contours . . . . . . . . . . . . . . . . . . . . . . . . 37 3.3 Diusive Representation of an Ideal Sphere Heat Transfer Model 38 3.3.1 Ideal Sphere Simulation by Using Finite Dierences Model 38 3.3.2 Diusive Representation Identication . . . . . . . . . . . 39 3.3.2.1 Frequency Response Analyse . . . . . . . . . . . 43 3.3.3 8-Poles Diusive Representation of Two ISFDMs with Different Physical Parameters . . . . . . . . . . . . . . . . . 47 3.3.3.1 Frequency Response . . . . . . . . . . . . . . . . 47 3.3.3.2 Open Loop Time Response . . . . . . . . . . . . 48 3.3.3.3 Closed Loop Time Response . . . . . . . . . . . 49 4 8-Poles Diusive Representation from Laboratory Data 53 4.1 8-Poles Diusive Symbol from Laboratory Data . . . . . . . . . . 53 4.2 Frequency Response . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.3 Open Loop Time Response . . . . . . . . . . . . . . . . . . . . . 56 4.4 Closed Loop Time Response . . . . . . . . . . . . . . . . . . . . . 56 5 Results 58 6 6 Conclusions and Future Development 58 Budget 59 Environment Impact 59 References 60 Appendix I: How to Find an Isolated Pole by Means of Diusive Representation 62 Glossary 68 7 List of Figures 1 Block diagram of discrete diusive representation Laplace transform................................... 17 2 Position of theoretical diusive symbol deltas (a) and frequency response (b)(c), of ISTQM of 10 mm outer radio and 6 mm inner radio under two Martian winds of velocities 0.3 m/s and 0.8 m/s (solid magenta and dashed green, respectively). As a reference, in dotted red line (b)(c) it is also plotted when the outer surface from the sphere is connected to a heat sink (T out =0). . . . . . . 25 3 (a) shows the rst 8 deltas of the theoretical diusive symbol. (b)(c) show thermal impedance frequency response of theoretical ISTQM, in green, and frequency response from 8-poles theoretical diusive representation, in dotted blue. . . . . . . . . . . . . . . . 26 4 ISTQM theoretical diusive symbol (a) and frequency response (b)(c) when inner and outer radio is re-sized. In green, reference sphere. In red, when outer radio is decreased one millimetre. In blue, when inner radio is decreased one millimetre. And in magenta, when both radius are decreased one millimetre. . . . . 27 5 ISTQM theoretical diusive symbol (a) and frequency response (b)(c) when physical parameters are decreased. In green, reference sphere. In red, when sphere conductivity is decreased. In blue, when specic heat is decreased. And in magenta, when mass density is decreased. . . . . . . . . . . . . . . . . . . . . . 27 6 Cole-Cole symulink block diagram . . . . . . . . . . . . . . . . . 28 7 10 poles diusive symbols of Cole-Cole function without noise. (a)(b) show time input and output signal details of a PRBS. (c) shows diusive symbols in logarithmic scale: theoretical diusive symbol in red, diusive symbol obtained using step input signal, in green, and using PRBS, in blue. (d)(e) show frequency responses of this diusive symbols. And (f) shows relative error between theoretical transfer function and transfer function obtained from 10-poles DR by using step and PRBS, without noise 30 8 10 poles diusive symbols of the Cole-Cole function with 30dB SNR Gaussian white noise. (a) shows time output signal details using as input: step signal (in green) and PRBS (in blue). (b) shows the approximate output, using Cole-Cole block diagram of g. 6 in red , using DR obtained with step input signal in red, and with PRBS in blue. (c) shows diusive symbols in logarithmic scale: theoretical diusive symbol in red, diusive symbol obtained using step input signal, in green, and using PRBS, in blue. (d)(e) show frequency responses of this diusive symbols. And (f) shows relative error between theoretical transfer function and transfer function obtained from 10-poles DR by using step and PRBS, with 30dB SNR Gaussian white noise . . . . . . . . . 31 8 9 Transfer function relative error from Cole-Cole diusive representation of 10 poles and adding and removing poles at low and high frequencies. (a)(d) adding and removing poles at high frequencies. (b)(e) adding poles at low frequencies. (c)(f) two poles at high frequencies are removed while dierent number of poles at low frequencies are added. . . . . . . . . . . . . . . . . . . . . . . 32 10 Transfer function relative error from Cole-Cole diusive representation with dierent number of poles while bandwidth is preserved. (a)(d) changing number of poles whereas bandwidth remains constant. (b)(e) adding poles to the number of poles used in (a)(d). (c)(f) eect of pseudoinverse truncation when there is a high number of poles. . . . . . . . . . . . . . . . . . . . . . . . 34 11 Noise eect on diusive representation using several number of poles and appropriate bandwidth. The eect of ltering is also shown when a high number of poles (3+20) is used. (c)(f) full rank and rank 11 are shown. . . . . . . . . . . . . . . . . . . . . 36 12 7 diusive symbols of 2+10-poles DR, geometrically displaced 7 times, recovering theoretical Cole-Cole diusive symbol shape. . . 37 13 (a) Interpolated diusive symbols module of 1+10-poles DR using dierent contours. (b) Transfer function relative error using 1+10-poles DR with dierent contours . . . . . . . . . . . . . . . 37 14 Scheme of an elementary cell . . . . . . . . . . . . . . . . . . . . 38 15 Frequency response of adapted theoretical sphere under two Martian winds of 0.3 m/s (magenta and blue) and 0.8 m/s (green and yellow). Using nite dierences model (solid lines) and ISTQM thermal impedance transfer function (dashed lines). Frequency response is also presented when there is a heat sink over the sphere surface, in dotted red line. . . . . . . . . . . . . . . . . . . 41 16 Theoretical diusive symbol deltas of adapted ISTQM under two dierent Martian wind velocities. In magenta, 0.3 m/s wind velocity, and in green, 0.8 m/s wind velocity. . . . . . . . . . . . 41 17 20 DS of 10-poles DR, 20 times displaced. In magenta from ISFDM of 0.3 m/s wind velocity and, in green, from ISFDM of 0.8 m/s wind velocity. . . . . . . . . . . . . . . . . . . . . . . 42 18 (a)(b): frequency response of diusive representation using TLADS (dashed lines) and nite dierences (solid lines). (c)(solid lines): relative error between DR transfer function using TLADS and ISFDM transfer function. (c)(dashed lines): relative error between DR transfer function using diusive symbol 1 positon to right from TLADS in g. 17 and ISFDM transfer funtion. . . . 43 19 Relative error transfer function if the chosen diusive symbol is 2 positions to the left from TLADS (in dashed lines) and 5 positions to left from TLADS (in dotted lines). . . . . . . . . . . . . . . . . 44 20 (a)(b): frequency response of diusive representation using TLADS (dashed lines) and nite dierences (solid lines). (c)(dashed lines): relative error between the transfer function using 10-poles DR and the ISFDM transfer function. (c)(dotted lines): relative error between the transfer function 8-poles DR and the ISFDM transfer function................................. 45 9 domain Ω− γ should be big enough to contain all the singularities. In practise, it is chosen a contour of type: γ(ξ) = |ξ|ei sign(ξ)( π 2+α), ξ ∈R (3) with α∈]0,π 2] being small enough. However, to obtain a good approximation of the model, the smaller α is, the more expensive numerically the discretisation in ξ is. Therefore, if the results are good, greater α can be addressed. In the case of thermal systems all the poles of H(p) are placed in R− , then the contour with the greatest α=π 2 can be applied. Thereby, the state-space realisation of thermal system is obtained of the form: ∂tψ(t, ξ) = −ξψ(t, ξ) + u(t), ξ ∈R+, ψ(0, ξ)=0 y(t) = R+∞ 0η(ξ)ψ(t, ξ)dξ (4) where ψ is a time-frequency representation of u called the diusive representation -(DR). And η(ξ) is the diusive symbol -(DS) of H(∂t) , which is a solution of the following integral equation directly obtained from Laplace transform (with respect to t): H(iω) = Z+∞ 0 η(ξ) iω +ξdξ (5) When H(t, ∂t) is time-variant, the diusive realisation is extended to [15]: ∂tψ(t, ξ) = −ξψ(t, ξ) + u(t), ξ ∈R+, ψ(0, ξ)=0 y(t) = R+∞ 0η(t, ξ)ψ(t, ξ)dξ (6) In eqs. (4)(6) initial condition ψ(0, ξ)=0 is supposed . From linearity, the contribution of the initial condition ψ(0, ξ) can be separately expressed by the exponential form: ψ0(t, ξ) = e−ξtψ(0, ξ) (7) thereby, the function ψ(t, ξ) + ψ0(t, ξ) is solution of eqs. (4)(6) with initial condition ψ(0, ξ)6= 0 . 2.1 Finite-Dimensional Approximations Finite-dimensional approximations of H(∂t) can be constructed by discretising the ξ variable, leading to an input-output approximation u7→ ey≃y=H(∂)u , of the form    ˙ ψk(t) = −ξkψk(t) + u(t), ψk(0) = 0, k = 1,··· , K ey(t) = K P k=1 ηkψk(t) (8) which can be reinterpreted by dening the atomic diusive symbol eη(ξ) := K P k=1 ηkδ(ξ−ξk) , where δ(ξ) is the classical Dirac distribution, of the form: ∂tψ(t, ξ) = −ξψ(t, ξ) + u(t), ξ ∈R+, ψ(0, ξ)=0 ey(t) = R+∞ 0eη(ξ)ψ(t, ξ)dξ (9) 16 ηk can be computed by using standard quadrature methods like linear interpolation, allowing to deal with numerical approximations in a unied framework. Using the classical interpolation function the discretised diusive symbol is of the form: ηk=η(ξk)Z+∞ 0 Λk(ξ)dξ =η(ξk)(ξn+1 −ξn−1) 2 (10) In the same way, the nite-dimensional approximation of the time-variant H(t, ∂t) is:    ˙ ψk(t) = −ξkψk(t) + u(t), ψk(0) = 0, k = 1,··· , K ey(t) = K P k=1 ηk(t)ψk(t) (11) If a time-varying K-dimensional approximate diusive symbol bη(t, ·) is obtained and the system under consideration is time-invariant, then bη(t, ·) converges (under a suitable hypothesis) to the exact diusive symbol η , in the following sense: bη(t, ξ)−→ t→+∞eη(ξ)−→ K→+∞η(ξ) (12) By applying laplace transform of (8), the block diagram shown in g. 1 is obtained. This block diagram is coded in Simulink R  and used in the sections below to compute the diusive representation output. In this picture, the ltering eect of ξ and the weight eect of η can be observed. Figure 1: Block diagram of discrete diusive representation Laplace transform. 17 2.2 Identication of the Diusive Symbol in Time Domain The identication problem is to obtain the diusive symbol η(ξ) of the unknown operator H(∂t) from experimental data. Numerical solutions impose a discretisation {ξk}k=1:K of the variable ξ . Therefore the operator approximation is given by (see Eq. 6): y≃ K X k=1 ηkψk(t) (13) The data available on the temporal mesh is also nite, tm=1,···,M ∈R+ , and then the time-continuous approximate model can be written under matrix form: ym≃Am,kηk (14) where, ym=y(tm) and Am,k =ψk(tm) . Given a known input, u(t) , and an output measure with error eym:= y(tm)+ δy(tm) , the identication problem is to yield a solution, bηk , for the following equation: eym≃Am,kηk, m k (15) which can be approximated by the nite dimensional least-square approach: min bηk M X m=1 |(Am,k bηk)m−eym|2 (16) the solution is given by the pseudoinverse of the matrix A , A+ : bηk=A+ k,meym= (A∗ k,mAm,k)−1A∗ k,meym (17) The identication problem is often ill-conditioned because the matrix A∗ k,mAm,k is closed to a non-invertible one. This problem is classically solved by adding judiciously a penalisation term ε [1,3] : bηk= (A∗ k,mAm,k +εIk,k)−1A∗ k,meym (18) Another way to address this computational problem is based on singular value decomposition (SVD) which provides improved understanding of the stabilisation problem while it gives an accurate solution. The singular value decomposition of matrix A can be written as A=UWV∗ (19) where the columns of U are the eigenvectors of AA∗ , the rows of V∗ are the eigenvectors of A∗A , and W is the diagonal matrix of 'singular values' of A . These singular values are the square roots of the eigenvalues of both AA∗ and A∗A , which means that U and V share the same eigenvalues. If the singular value decomposition of A is UWV∗ then the pseudo-inverse or Moore-Penrose inverse of A is: A+=U∗W−1V (20) 18 If A+ has more rows than columns, which is the case, and has full rank, it result in the least squares solution seen before in eq. (17): A+=U∗W−1V= (A∗A)−1A∗ (21) Then we can obtain the diusive symbols of the form: bηk=A+ m,k eym=U∗k,kW−1 k0m−k,k0m−k Vm,meym (22) The SVD of Ak,m returns k singular values, σk , in descending order along the diagonal of W , the other m−k singular values are 0 . The fraction number between the largest and the smaller singular value, κ(A) = σmax σmin is called 'the condition number of A ', which, by denition, is equal or larger than one. This number is indicative of the magnication of the error in eq. (22): kδbηk kbηk≤κ(A)keyk kAbηkkδeyk keyk (23) where keyk kAbηk=1 cos θ , being θ the angle between the vector ey and the range of A . The closer is cos θ to 1, the closer is approximate solution trajectory to data trajectory, and the condition number of the least squares problem becomes closer to the condition number of a nonsingular matrix, κ(A) . Conversely, if cos θ decreases towards zero the trajectory of the approximate solution moves away from the trajectory data, and the condition number of the least squares problem becomes arbitrarily large. We want to bound the error as much as possible in the determination of the diusive symbols trajectory, kbηk , to ensure that small changes in the input data do not aect our model. To do that, we should ensure that θ is as close as possible to 0 and, at the same time, that κ(A) is as small as possible. This is an important drawback in this identication method, since usually frequency bandwidth covers several decades, i.e. ξ1ξk , which makes κ(A) big from the beginning. Moreover, if discretisation points number, K , increases, θ decreases to 0 as it is desirable, but, since more linear dependence between trajectories appears when K increases, more singular values close to 0 also appear making that κ(A) increases, which is not desirable. Therefore, a balance in the choice of the number of discretisation points must be found to avoid an ill-conditioned problem. Singular values in W close to 0 (large values in W−1 ) make the noise and numerical errors on the transformation from ey -space to bη -space be increased. Singular values close to 0 can be made here equal to 0 , thereby we 'low-pass lter' the input vector, ey . Thus, we can lter the noise in the measurement, but in this case the angle θ will increase arbitrarily. Again, in this case, a judiciously tolerance should be chosen. In this work, we are going to use the Matlab R  function pinv() which is based on Moore-Penrose pseudoinverse and this function allows setting any singular value less than certain tolerance to zero. We are going to work only with full-rank matrices (without being limited by tolerance) since, as it is seen in sec. 3.2, ltering pseudoinverse means using more poles without any remarkable improvement. And moreover, only a small number of discretisation points are necessary to obtain good results, which is 19 also appropriate to reduce computational cost. 3 Numerical Examples The spherical thermal sensor which is being developed in the Electronic Engineering department is of complex nature and Diusive Representation seems a good candidate to obtain an accurate reduced order model of that sensor, as it was discussed in sec. 1. Some ideal non-rational heat transfer models close to this spherical thermal sensor are proposed in order to understand the real model behaviour and its approximation to the ideal one. These ideal models will allow us to know if the Diusive Representation is suitable to model the aforementioned spherical thermal sensor trought some numerical examples, showing its advantages and drawbacks. 3.1 Ideal Sphere Thermal Quadrupole Model An ideal sphere thermal quadrupole model (ISTQM) is proposed here in order to perform a numerical example of the diusive representation identication problem. This is done to have a similar model to the aforementioned spherical thermal anemometer prototype. Let us consider the problem of computing the heating temperature Tin(t) at the inner surface of a hollow sphere where a random heat ux Pin(t) is applied at this surface (this is the operation mode of the aforementioned spherical thermal anemometer prototype). Assuming heat conduction into the sphere is isotropic, the heating temperature T(r, t) and the ux passing through the sphere P(r, t) satisfy the following heat diusion equations in spherical coordinates: ∂T(r, t) ∂t =a∂2T(r, t) ∂r2+2 r×∂T(r, t) ∂r  (24) P(r, t) = −λ(4πr2)∂T(r, t) ∂r (25) with, •a=λ/ρc : thermal diusivity (m 2 s −1 ) •λ : thermal conductivity (Wm −1◦ C −1 ) •ρ : mass density (kg m −3 ) •c : specic heat capacity (J kg −1◦ C −1 ) A thermal quadrupole can be used in order to model heat diusion into the sphere [7,9]: Tin(s) Pin(s)=A B C D  Tout(s) Pout(s) (26) 20 Using eqs. (24) and (25) the terms A, B, C, D of the matrix characterising the heat transfer of the sphere of inner radio r1 and outer radio r2 (see Fig.) are:            A=r2 r1cosh φ−sinh φ δr1 B=1 4πr1r2 sinh φ δλ C= 4πr2λh1−r1 r2cosh φ+δr1−1 δr2sinh φi D=r1 r2cosh φ+sinh φ δr2 (27) where φ=δ(r2−r1) , and δ=ps/a . Considering that the sphere sensor is sensing a wind and all temperatures are referenced to room temperature, the outer heat ux Pout(t) can be related with the outer temperature Tout(t) through a wind heat transfer coecient hw , of the form: Pout(t) = hwATout(t) (28) this heat transfer coecient hw (W m −2◦ C −1 )) can be computed by using the Nusselt number Nu of a forced convention heat transfer of a uid over the sphere [21]: hw=Nuλair dsp (29) where λair is the air thermal conductivity, dsp = 2r2 the outer diameter of the sphere, and Nu is: Nu = 2 + (0.4Re1 2+ 0.06Re2 3)Pr0.4µ∞ µs1 4 (30) where µs is the air viscosity at sphere surface temperature, µ∞ is the air viscosity at room temperature, and Re and Pr are the Reynolds and Prandtl numbers respectively: Re =ρairvairdsp µair (31) Pr =µaircair λair (32) where ρair is the air mass density, vair the air velocity, dsp = 2r2 the outer diameter of the sphere, µair the air viscosity which is temperature dependent and is approximated here, cair the air specic heat capacity, and λair the air thermal conductivity. Therefore, using eq. (28) in thermal quadrupole with the terms of (27), the transfer function of the system, its thermal impedance, can be computed: H(s) = Tin(s) Pin(s)=1 4πr2 r2 r1+1 δr1r2hw λ−1tanh φ λ1−r1 r2+hwr1+hλδr1−1 δr2+hw δitanh φ (33) 21 The thermal resistence Rths is obtained by considering the static gain: lim s→0 Tin(s) Pin(s)=Rths =1 4πλ 1 r1−1 r2+1 4πr2 2hw (34) At high frequencies, transfer function of the thermal impedance behaves like a non-integer integrator with order equal to 0.5: lim s→∞Hs(s) = a0.5 4πr2 1λs0.5 (35) In the case of sphere surface is a heat sink ( hw→ ∞ , Tout ≃Tair ), thermal impedance and thermal resistance are [9]: lim hw→∞H(s) = Hs(s) = 1 4πr1λ 1 1 + δr1tanh−1φ (36) lim hw→∞Rths =1 4πλ 1 r1−1 r2 (37) In this case, at relative high frequencies, the transfer function of the thermal impedance behaves as a Cole-Cole function whose order is equal to 0.5: lim s→∞Hs(s) = 1 4πr1λ 1 1 + r2 1 as0.5 (38) which, at high frequencies, behaves again like the non-integer integrator seen in eq. (35). The theoretical diusive symbol associated with Cole-Cole transfer function H(s)=1/(1 + (τ0s)α) is known [15]: η(ξ) = sin(απ)/π (τ0ξ)−α+ 2 cos(απ)+(τ0ξ)α (39) where α is the Cole-Cole order and τ0 is the characteristic relaxation time. Because of this, we are going to perform a numerical simulation based on nonrational Cole-Cole transfer function in sec. 3.2. Then, in sec. 3.3, numerical simulations of the ISTQM proposed here is performed by using nite dierences model. 22 3.1.1 Theoretical Diusive Symbol of the ISTQM thermal impedance The theoretical diusive symbol of the ISTQM thermal impedance is obtained here in order to verify diusive symbols computed in simulations. From the thermal impedance quadrupole model, eq. (33), and following the method explained in [18] we have that the theoretical diusive symbol can be computed by using: η(ω) = 1 2πi lim ν→0[H(−ω−iν)−H(−ω+iν)] = =1 8π2r2i r2 r1+√a −i√ωr1r2h λ−1tanh −i√ω(r2−r1) √a λ1−r1 r2+hr1+hλ−i√ωr1 √a−√a −i√ωr2+h√a −i√ωitanh −i√ω(r2−r1) √a −1 8π2r2i r2 r1+√a i√ωr1r2h λ−1tanh i√ω(r2−r1) √a λ1−r1 r2+hr1+hλi√ωr1 √a−√a i√ωr2+h√a i√ωitanh i√ω(r2−r1) √a (40) tanh ix =itan x , and tanh −ix =−itan x , then theoretical diusive symbol is: η(ω) = 1 2πi lim ν→0[H(−ω−iν)−H(−ω+iν)] = =1 8π2r2i ihr2 r1+√a √ωr1r2h λ−1tan √ω(r2−r1) √ai ihλ1−r1 r2+hr1+−√ωλr1 √a−λ√a √ωr2+h√a √ωtan √ω(r2−r1) √ai −1 8π2r2i −ihr2 r1+√a √ωr1r2h λ−1tan √ω(r2−r1) √ai −ihλ1−r1 r2+hr1+−√ωλr1 √a−λ√a √ωr2+h√a √ωtan √ω(r2−r1) √ai =X ω∈R δλ1−r1 r2+hr1+−√ωλr1 √a−λ√a √ωr2 +h√a √ωtan √ω(r2−r1) √a =X ω∈R δ√ω√a[λ(r2−r1) + hr1r2] ωλr1r2+a(λ−hr2)−tan √ω(r2−r1) √a (41) which is nonzero only in the intersection points of the following curves: √ω√a[λ(r2−r1) + hr1r2] ωλr1r2+a(λ−hr2)= tan √ω(r2−r1) √a (42) , which position an innite discrete diusive symbol. 23 3.1.1.1 Simulation of ISTQM Theoretical Diusive Symbol We are going to simulate the ISTQM under Martian atmosphere to compare results with laboratory data. Since Martian atmosphere is 95.32% carbon dioxide and average temperature is about minus 60 ◦ C, we are going to use the Holman table for CO 2 [12] and eq. 29 to compute approximate convective heat transfer coecients. Moreover, sphere prototype is made of silver, which has the following physical parameters:                rin = 3 ·10−3m rout = 5 ·10−3m λ= 429 Wm −1 K −1 ρ= 10.49 ·103 kg m −3 c= 0.24 ·103 J kg −1 K −1 a= 1.704 ·10−4 m 2 s −1 (43) Frequency responses of ISTQM thermal impedance, eq. (33), with parameters (43), under two Martian wind velocities of 0.3 m/s and 0.8 m/s are shown in gs. 2(b)(c). It can be seen that there is a concordance between intersection points of curves from eq. (42), plotted in g. 2(a), and frequency response of thermal impedance, gs. 2(b)(c). It can be seen that, rst isolated intersection points (rst deltas), from each wind velocity shown in g. 2(a), match the cuto frequencies in Bode plots, gs. 2(b)(c). This pole for 0.3 m/s wind velocity is labelled in gs. 2(a)(b). From this pole to reach the next pole, magnitude decays 20dB per decade, which corresponds to frequency response of a rst order pole. However, beyond second pole magnitude decays 10dB per decade, which corresponds to frequency response of a non-integer integrator with order equal to 0.5. Thermal impedance after second pole tends to behave as ISTQM thermal impedance with surface connected to a heat sink, eq. (36). This frequency response behaviour is plotted in red dotted line in gs. 2(b)(c). Notice that tangent curve from eq. (42) tan √ω(r2−r1) √a (44) plotted in cyan in g. 2(a) does not depend on heat transfer coecient (wind velocity) but only on physical sphere parameters. Also notice that the curve from eq. (42) √ω√a[λ(r2−r1) + hr1r2] ωλr1r2+a(λ−hr2) (45) plotted in solid magenta and dashed green in g. 2(a) does depend on heat transfer coecient (wind velocity) and resembles a Cole-Cole diusive symbol whose order is α= 0.5 , eq. (39). As we have seen before, ISTQM thermal impedance transfer function behaves as Cole-Cole transfer function at high frequencies. This behaviour can also be seen, at high frequencies, in g. 2(a) where singularities are close together recovering the shape of eq. (45). 24 It is remarkable that the rst isolated pole, with lower frequency, depends on heat transfer coecient (wind velocity), whereas the other poles, at high frequencies, almost only depend on sphere parameters. Moreover, at high frequencies, intersection points, which are close together, behave as Cole-Cole transfer function which is a continuous function. In gs. 2-(b)(c), the dotted red line shows frequency response when convective heat transfer coecient is innite, that is to say, sphere surface is connected to a heat sink (T out =0). In this case, frequency response only depends on sphere parameters. Therefore, at high frequencies, frequency response only depends on sphere physical parameters, at low frequencies, frequency response only depends on wind velocity, and at central frequencies, frequency response depends on both. 10-2 100102 Frequency(rad· s -1) 0 0.05 0.1 0.15 0.2 0.25 Intersection points (a) Diffusive symbol deltas 10-2 100102 Frequency (rad· s -1) -50 -40 -30 -20 -10 0 10 20 30 40 Magnitude [K/W] (dB) (b) Thermal impedance Ideal sphere wind 0.3 m/s Ideal sphere wind 0.8 m/s Ideal sphere conductivity inf 10-2 100102 Frequency (rad· s -1) -80 -70 -60 -50 -40 -30 -20 -10 0 Phase (degree) (c) Thermal impedance X: 0.01024 Y: 0.0155 X: 0.01042 Y: 36.41 Figure 2: Position of theoretical diusive symbol deltas (a) and frequency response (b)(c), of ISTQM of 10 mm outer radio and 6 mm inner radio under two Martian winds of velocities 0.3 m/s and 0.8 m/s (solid magenta and dashed green, respectively). As a reference, in dotted red line (b)(c) it is also plotted when the outer surface from the sphere is connected to a heat sink (T out =0). In order to verify that theoretical diusive symbol has been well found, frequency response of ISTQM is recovered from it. Fig. 2(a) shows at which frequencies are placed some of the innite poles of theoretical diusive symbol found in eq. (41) . Therefore, it is possible to use a subset of these poles, [ξ1,··· , ξj] , covering some bandwidth, to construct a nite diusive representation with the same frequency response inside this bandwidth. Since the frequency response at these poles H(ξj) is known by eq. (33), it is possible to dene an equation system of this nite diusive representation: η1 iξj+ξ1 +···+ηk iξj+ξk =H(ξj),(j=k) (46) This equation system is computed to nd 8-poles discrete diusive symbol from the ideal sphere under 0.3 m/s wind velocity. The result is shown in the following table : 25 not increase too much, but removing tree poles, the angle and relative error are considerable. However, when bandwidth is increased by adding poles at low frequencies (x+10) (gs. 9-(b)(e)) the condition number does not increase so quickly and the angle improves a little bit (see table 4). It happens until two poles are added, getting better matching with theoretical diusive symbols and transfer function. But, if more than two low frequency poles are added, the angle does not improve whereas the condition number is still increasing, and if more than four poles are added, the result worsens. Finally, taking into account that removing two poles at high frequencies transfer function relative error is similar, as it have been shown before. Figs. 9- (c)(f) show the case when these two poles at high frequencies are removed and poles at low frequencies are added. Contrasting relative error between gs. 9-(e) and (f), it is noticed that, it is possible to obtain similar or better results by using fewer poles and placing them correctly. 10-1 100101102 ξ (rad s-1) 100 105 η(ξ)k (a) Varying diffusive symbols at high frequencies Theoretical 10 10+1 10+2 10-1 10-2 10-3 10-2 10-1 100101102 ξ (rad s-1) 10-4 10-2 100 η(ξ)k (b) Adding diffusive symbols at low frequencies Theoretical 10 1+10 2+10 3+10 4+10 10-2 10-1 100101 ξ (rad s-1) 10-4 10-2 100 η(ξ)k (c) Diffusive symbols reducing bandwidth Theoretical 8 1+8 2+8 3+8 10-1 100101 frequency ω (rad s-1) 0 5 10 15 20 % (d) Relative error |Hη - Htheo| / |Htheo| 10 10+1 10+2 10-1 10-2 10-3 10-2 10-1 100101 frequency ω (rad s-1) 0 2 4 6 8 10 12 % (e) Relative error |Hη - Htheo| / |Htheo| 10 1+10 2+10 3+10 4+10 10-2 10-1 100101 frequency ω (rad s-1) 0 2 4 6 8 10 12 % (f) Relative error |Hη - Htheo| / |Htheo| 8 1+8 2+8 3+8 Figure 9: Transfer function relative error from Cole-Cole diusive representation of 10 poles and adding and removing poles at low and high frequencies. (a)(d) adding and removing poles at high frequencies. (b)(e) adding poles at low frequencies. (c)(f) two poles at high frequencies are removed while dierent number of poles at low frequencies are added. 32 DR poles number bandwidth (rad/s) θ (degree) κ(A) 10 [3.14 ·10−2,6.28 ·10+2] 0.2471 3.8068e+05 10 + 1 [3.14 ·10−2,1.89 ·10+3] 0.2471 2.1672e+07 10 + 2 [3.14 ·10−2,5.67 ·10+3] 0.2471 8.0364e+11 10 - 1 [3.14 ·10−2,2.09 ·10+2] 0.3290 5.8238e+04 10 - 2 [3.14 ·10−2,6.9·10+1] 0.3522 1.1867e+04 10 - 3 [3.14 ·10−2,2.3·10+1] 1.1809 3.7029e+03 1 + 10 [1.05 ·10−2,6.28 ·10+2] 0.0100 8.5754e+05 2 + 10 [3.48 ·10−3,6.28 ·10+2] 0.0059 1.4348e+06 3 + 10 [1.16 ·10−3,6.28 ·10+2] 0.0059 1.9703e+06 4 + 10 [3.85 ·10−4,6.28 ·10+2] 0.0059 2.4346e+06 (8 low bw) [3.14 ·10−2,6.28 ·10+1] 0.3483 1.0848e+04 1+(8 low bw) [1.06 ·10−2,6.28 ·10+1] 0.2470 2.4275e+04 2+(8 low bw) [3.58 ·10−3,6.28 ·10+1] 0.2469 4.0667e+04 3+(8 low bw) [1.20 ·10−3,6.28 ·10+1] 0.2469 5.5992e+04 Table 4: Projection angle and condition number from Cole-Cole diusive representation of 10 poles and adding and removing poles at low and high frequencies. In gs. 10-(a)(d) and table 5, it is shown what happens when DR poles number is modied inside the same bandwidth. On the one hand, increasing the number of poles, the angle decreases, and the relative error of transfer function also decreases in a well approximated bandwidth. It can be seen that with more than 5 poles approximation becomes quite good. However, from eight poles onwards improvement is not signicant. On the other hand, increasing the number of poles the condition number also increases making solution more sensitive to noise and numerical errors (see gs. 11-(c)(f)). In simulations it has been found that, the appropriate number of poles actually depends on the distance between them. Moreover, a geometric scale shows good behaviour to place those poles. Experimentally, and in this numerical example, it is found that using poles geometrically scaled the appropriate scale ratio is between 2 and 8. The best results are obtained with a geometric scale ratio between 3 and 6. Figs. 10-(b)(e) show that, when some extra poles are added at low frequencies, extending bandwidth to 6·10−3 , the angle decreases at the expense of increasing a little bit the condition number (see table 5). This improves relative error at small frequencies using more than 4 poles. Table 5 shows that, in the case of 3+20 poles and truncating pseudoinverse, the condition number decreases while the angle increases. The rank of pseudoinverse is reduced by truncation, ltering correlated trajectories. In g. 10-(c) it is shown that, approximate diusive symbols move away from theoretical diusive symbols whereas transfer function relative error (g. 10-(f) and g. 11-(f)) remains similar when pseudoinverse is truncated. And, if truncation is too big (3+20 rank 11), the relative error increases. 33 Therefore, choosing a small number of poles results similarly to using more poles and truncating pseudoinverse matrix are obtained. But, in this case, using an appropriate small number of poles, identication problem is more stable and has low computational cost. 10-2 10-1 100101102 ξ (rad s-1) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 η(ξk) (a) Cole-Cole diffusive symbols Theo 4 5 6 8 20 30 10-2 10-1 100101102 ξ (rad s-1) 0 0.05 0.1 0.15 η(ξk) (b) Cole-Cole diffusive symbols Theo 1+4 1+5 1+6 1+8 3+20 3+30 10-2 10-1 100101102 ξ (rad s-1) 0.05 0.1 0.15 0.2 0.25 η(ξk) (c) Cole-Cole diffusive symbols Theo 3+20 full rank 3+20 rank 22 3+20 rank 19 10-1 100101 frequency ω (rad s-1) 0 5 10 15 % (d) Relative error |Hη - Htheo| / |Htheo| 4 5 6 8 20 30 10-1 100101 frequency ω (rad s-1) 0 5 10 15 % (e) Relative error |Hη - Htheo| / |Htheo| 1+4 1+5 1+6 1+8 3+20 3+30 10-2 10-1 100101 frequency ω (rad s-1) 0 5 10 15 % (f) Relative error |Hη - Htheo| / |Htheo| 3+20 full rank 3+20 rank 22 3+20 rank 19 Figure 10: Transfer function relative error from Cole-Cole diusive representation with dierent number of poles while bandwidth is preserved. (a)(d) changing number of poles whereas bandwidth remains constant. (b)(e) adding poles to the number of poles used in (a)(d). (c)(f) eect of pseudoinverse truncation when there is a high number of poles. 34 DR poles number tolerance θ degree κ(A) 4 full rank 1.9134 3.0490e+04 5 full rank 0.7188 5.1534e+04 6 full rank 0.3900 8.1121e+04 8 full rank 0.2867 1.7664e+05 20 full rank 0.0703 2.8383e+07 30 full rank 0.0146 2.8330e+09 1+4 full rank 1.9120 1.1965e+05 1+5 full rank 0.6643 1.8298e+05 1+6 full rank 0.2701 2.5860e+05 1+8 full rank 0.0401 4.6326e+05 3+20 full rank 0.00062 8.7480e+07 3+30 full rank 0.0011 6.4168e+09 3+20 full rank 0.00062 8.7480e+07 3+20 rank 22 0.0030 6.9622e+06 3+20 rank 19 0.0839 1.5307e+05 Table 5: Projection angle and condition number from Cole-Cole diusive representation with dierent number of poles while bandwidth is preserved. And also, when pseudoinverse with a high number of poles is ltered (rank x). Fig. 11 shows inuence of the noise in system identication. It can be seen that using a small number of poles, means a small condition numbers (see table. 6), values of poles found do not change too much, gs. 11-(a)(b). However, when the condition number is high, the values of diusive symbol can change a lot, making identication problem unstable. This can be seen in g. 11-(c) where 3+20 full rank Cole-Cole diusive symbols are drawn with noise and without noise. Figs. 11(c)(f) also shows eect of ltering this 3+20 full rank pseudoinverse matrix to rank 11, the transfer function relative error is worse than in the case that there is not noise since the angle increases, see table. 6. However, there is not variation when noise is added since the condition number is lower, see table. 6. Result ltering the pseudoinverse to 11 is worse than the result obtained with 1+10 full rank pseudoinverse seen in gs. 11-(b)(e). Therefore, it seems better to use an appropriate number of poles to obtain a good result resilient to noise than using a lot of number of poles and then ltering. 35 100 ξ (rad s-1) 0.04 0.06 0.08 0.1 0.12 0.14 0.16 η(ξk) (a) Cole-Cole diffusive symbols Theor 1+8 no noise 1+8 30dB SNR 1+8 rbw no noise 1+8 rbw 30dB SNR 100105 ξ (rad s-1) 0 0.05 0.1 0.15 0.2 η(ξk) (b) Cole-Cole diffusive symbols Theor 10 no noise 10 30dB SNR 1+10 no noise 1+10 30dB SNR 2+10 no noise 2+10 30dB SNR 10-2 100102 ξ (rad s-1) 0 0.5 1 1.5 η(ξk) (c) 3+20 Cole-Cole diffusive symbols Theor full rank no noise full rank 30dB SNR rank 11 no noise rank 11 30dB SNR 10-1 100101 frequency ω (rad s-1) 1 2 3 4 5 % (d) Relative error |Hη - Htheo| / |Htheo| 1+8 no noise 1+8 30dB SNR 1+8 rbw no noise 1+8 rbw 30dB SNR 10-1 100101 frequency ω (rad s-1) 1 2 3 4 5 % (e) Relative error |Hη - Htheo| / |Htheo| 10 no noise 10 30dB SNR 1+10 no noise 1+10 30dB SNR 2+10 no noise 2+10 30dB SNR 10-1 100101 frequency ω (rad s-1) 0 1 2 3 4 5 6 % (f) Relative error |Hη - Htheo| / |Htheo| 3+20 full no noise 3+20 full 30dB SNR 3+20 rank 11 no noise 3+20 rank 11 30dB SNR Figure 11: Noise eect on diusive representation using several number of poles and appropriate bandwidth. The eect of ltering is also shown when a high number of poles (3+20) is used. (c)(f) full rank and rank 11 are shown. DR poles number noise (SNR) θ degree κ(A) 1+8 - 0.0401 4.6326e+05 1+8 (reduced bw) - 0.2470 2.4275e+04 10 - 0.2471 3.8068e+05 1+10 - 0.0100 8.5754e+05 2+10 - 0.0059 1.4348e+06 3+20 - 0.00062 8.7480e+07 3+20 rank 11 - 0.9302 1.2895e+04 1+8 30dB 3.6954 4.6423e+05 1+8 (reduced bw) 30dB 3.7033 2.4325e+04 10 30dB 3.7048 3.8148e+05 1+10 30dB 3.6953 8.5933e+05 2+10 30dB 3.6951 1.4378e+06 3+20 30dB 3.6946 8.7662e+07 3+20 rank 11 30dB 3.8118 1.2895e+04 Table 6: Noise eect on projection angle and condition number using several number of poles and appropriate bandwidth. The eect of ltering is also shown in the case of 3+20 poles In g. 12 it is shown that using diusive representation with a small condition number and a small projection angle (for example 2+10 in table 6), it is possible to recover the shape of theoretical DS. In order to do that, several discrete diusive symbols geometrically displaced are inferred, then those diffusive symbols are interpolated, eq. (10), and nally, all those interpolated DS are drawn together. 36 10-2 10-1 100101102103 Frequency (rad· s -1) -0.04 -0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 η(ξk) Interpolated diffusive symbols Theoretical Cole-Cole DS. 7 DSs, of 2+10 poles DR, displaced. Figure 12: 7 diusive symbols of 2+10-poles DR, geometrically displaced 7 times, recovering theoretical Cole-Cole diusive symbol shape. 3.2.3 Other Contours Fig. 13-(a) shows a diusive symbol of 11 poles DR using several contours eq. (3). It can be seen in g. 13-(b) that, there is not any signicant improvement in transfer function by using other dierent contours of −|ξ| . Moreover, calculus complexity is greater and relative error becomes greater as contour approaches to imaginary axes. In table 7, it is also shown that although condition number improves, projection angle between diusive representation and theoretical outputs gets worse. Therefore, in thermal systems, since poles are expected to be found on negative real axis, contour over negative real axis is used in this work simulations. 10-2 100102 Frequency ξ (rad s-1) 0.04 0.06 0.08 0.1 0.12 0.14 0.16 |η(ξk)| (a) Cole-Cole diffusive symbols Theor contour -|ξ| contour |ξ|e-i5π/6 contour |ξ|e-i3π/4 contour |ξ|e-i4π/6 10-2 10-1 100101 frequency ω (rad s-1) 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 % (b) Relative error |Hη - Htheo| / |Htheo| contour -|ξ| contour |ξ|e-i5π/6 contour |ξ|e-i3π/4 contour |ξ|e-i4π/6 Figure 13: (a) Interpolated diusive symbols module of 1+10-poles DR using dierent contours. (b) Transfer function relative error using 1+10-poles DR with dierent contours 37 DR poles number Contour θ degree κ(A) 1+10 -1 0.0100 8.5754e+05 1+10 e−5i/6 9.8682 6.3755e+05 1+10 e−3i/4 14.6439 4.4843e+05 1+10 e−4i/6 18.7653 2.8468e+05 Table 7: Diusive representation of 1+10 poles: projection angles and condition numbers using dierent contours 3.3 Diusive Representation of an Ideal Sphere Heat Transfer Model The laboratory data shown in sec. 4 are obtained by simulating two Martian wind velocities of 0.3 and 0.8 m/s over the spherical anemometer prototype presented in sec. 1 into an hypobaric chamber [14]. This laboratory data is sampled for 10 4 s, and sample time is 0.05s. In this numerical example, an PRBS input of the same duration and sample time is applied to one sensor under the same dierent wind velocities of 0.3 and 0.8 m/s. Therefore, response can be correctly approximated only in the frequency band: 2π tK ,2π 2max{∆tk}=6.28 ·10−4,6.28 ·10+1 (51) 3.3.1 Ideal Sphere Simulation by Using Finite Dierences Model In order to obtain time response of ISTQM, an ideal sphere numerical simulation is performed using nite dierences model (FD) [11]. For this model, the sphere is divided into I spherical slices of thickness ∆r . The radius of each slice is de- ned by ri= (I0+i)∆r where I0=r1/∆r and 0⩽i⩽I , the elementary surface is Si= 4π∆r2(I0+i)2 , the elementary thermal resistance Ri= 1/4πλ∆r(I0+i)2 , and the elementary thermal capacity Ci=4 3π∆r3[3(I0+i)2+1 4]ρc . Considering each slice as the elementary cell shown in g. 14, temperature is governed by: dTi dt =1 RiCi (Ti−1−2Ti+Ti+1) (52) Figure 14: Scheme of an elementary cell 38 Since the sphere is immersed into dierent air ows, outer resistance is replaced by corresponding convective thermal resistance, Rh , and connected to room temperature, Tr . Therefore, since inner temperature Tin(t) is observed while power ux is applied to the sphere inner surface Pin(t) , the ideal sphere state-space representation using nite dierences is: ˙xs=Asxs+BsPin(t) Tin(t) = Csxs (53) where xs= [T0, T1,··· , Ti,··· , TI−1, Tr] , Cs= [1,0,··· ,0] , Bs=6a ∆r2      1 4πλ∆r(3I2 0+1 4) 0 . . . 0      (54) and the non-null elements of matrix As are:                    As(i, i −1) = a ∆r2 (I0+i)2 (I0+i)2+1 12 2⩽i⩽I−1 As(i, i) = a ∆r2 (I0+i)2+(I0+i+1)2 (I0+i)2+1 12 2⩽i⩽I−1 As(i, i + 1) = a ∆r2 (I0+i+1)2 (I0+i)2+1 12 2⩽i⩽I−1 (55)                            As(1,1) = −a ∆r2 (I0+1)2 I2 0 2+I0 2+1 24 , As(1,2) = a ∆r2 (I0+1)2 I2 0 2+I0 2+1 24 As(I, I −1) = 2a ∆r2 (I0+I)2 (I0+I)2−(I0+I) 2+1 12 , As(I, I) = −2a ∆r2 (I0+I)2 (I0+I)2−(I0+I) 2+1 12 −2a ∆r34πλRh[(I0+I)2−(I0+I) 2+1 12 ], As(I, I + 1) = 2a ∆r34πλRh[(I0+I)2−(I0+I) 2+1 12 ], (56) 3.3.2 Diusive Representation Identication Martian atmosphere is 95.32% CO 2 and average temperature is about minus 60 ◦ C. Therefore, the Holman table for CO 2 at minus 53 ◦ C (220 ◦ K) [12] is used to compute the convective heat transfer coecients, eq. (29). Convective heat transfer coecients vary with temperature and pressure, but only an approximate value is needed for this work. Moreover, since the ideal sphere diers from the complex anemometer prototype, some physical parameters of the theoretical silver sphere (43), are changed to obtain frequency response similar to that in the laboratory (g. 33). For example, since the laboratory sphere prototype is split into two/three hemispheres, 39 the ideal sphere radius may be smaller. Then, the outer radio is reduced to 3.5 mm, and the inner radio is conned to 0.5 mm since laboratory sensor/heaters does not cover all the hemisphere inner surface. With this outer radio, convective heat transfer coecients, hw , are approximated to 20 W/m 2 K under 0.3 m/s wind velocity, and 32 W/m 2 K under 0.8 m/s wind velocity. Moreover, since there is a thermal barrier (thermal-conductance paste of thermal conductivity 3.4 Wm−1K−1 ) and several discontinuities between sensor/heater and inner surface of hemispheres, the ideal sphere thermal conductivity is divided by 150. Setting these new parameters,                rin = 0.5·10−3m rout = 3.5·10−3m λ= 2.86 Wm −1 K −1 ρ= 10.49 ·103 kg m −3 c= 0.24 ·103 J kg −1 K −1 a= 1.136 ·10−6 m 2 s −1 (57) , transfer function response obtained from thermal quadrupole model under two winds of 0.3 m/s and 0.8 m/s (g. 15, magenta-blue and green-yellow, respectively) are close to that obtained in the laboratory (g. 33). Using nite dierences model described in 3.3.1, it is possible to compute thermal time response from the ideal sphere, which is needed to obtain its diusive representation. This ideal sphere nite dierences model (ISFDM) is veried comparing its thermal impedance transfer function frequency response with the ISTQM transfer function frequency response obtained in eq. (33). Bode plots from nite dierences model (solid lines) and from thermal quadrupole model transfer function (dashed lines), eq. (33), are shown in g. 15. As it can be seen, frequency response using both models is quite similar. Two wind velocities of 0.3 and 0.8 m/s are simulated, and the sphere is divided into 500 slides to compute nite dierences model. The ISTQM frequency response when instead of wind there is a heat sink is also shown in dotted red. In this case, frequency response only depends on sphere parameters. 40 10-3 10-2 10-1 100101 Frequency (rad· s -1) 20 25 30 35 40 45 50 Magnitude [K/W] (dB) (b) Thermal impedance Finite differences wind 0.3 m/s Finite differences wind 0.8 m/s Quadrupole model wind 0.3 m/s Quadrupole model wind 0.8 m/s Quadrupole model h=Inf 10-3 10-2 10-1 100101 Frequency (rad· s -1) -50 -40 -30 -20 -10 0 Phase (degree) (c) Thermal impedance Figure 15: Frequency response of adapted theoretical sphere under two Martian winds of 0.3 m/s (magenta and blue) and 0.8 m/s (green and yellow). Using nite dierences model (solid lines) and ISTQM thermal impedance transfer function (dashed lines). Frequency response is also presented when there is a heat sink over the sphere surface, in dotted red line. Theoretical diusive symbol deltas from this adapted ISTQM under two dierent Martian wind velocities are computed and shown in g. 16 and table 8. As expected, the lower frequency pole for each wind velocity has a dierent value whereas the other poles have similar values for the two wind velocities. 10-2 10-1 100101102 Frequency (rad· s -1) 0.2 0.4 0.6 0.8 1 1.2 Intersection point (a) Diffusive symbols deltas Figure 16: Theoretical diusive symbol deltas of adapted ISTQM under two dierent Martian wind velocities. In magenta, 0.3 m/s wind velocity, and in green, 0.8 m/s wind velocity. 41 10-2 10-1 100101 Frequency(rad· s -1) 0 10 20 30 40 50 Magnitude [K/W] (dB) (a) SFDM DR thermal impedance 8 poles DR adapted SFDM, wind 0.3 m/s 8 poles DR adapted SFDM, wind 0.8 m/s 8 poles DR higher conductivity, wind 0.3 m/s 8 poles DR higher conductivity, wind 0.8 m/s 10-3 10-2 10-1 100101 Frequency(rad· s -1) -60 -50 -40 -30 -20 -10 0 Phase (degree) (b) SFDM DR thermal impedance 10-3 10-2 10-1 100 Frequency(rad· s -1) 0 1 2 3 4 5 % (c) Relative error |Hη - Htheo| / |Htheo| 8 poles DR adapted SFDM, wind 0.3 m/s 8 poles DR adapted SFDM, wind 0.8 m/s 8 poles DR higher conductivity, wind 0.3 m/s 8 poles DR higher conductivity, wind 0.8 m/s Figure 25: (a)(b) frequency response of 8-poles DR transfer function from two dierent ISFDM and under two wind velocities. In magenta: reference sphere under wind velocity of 0.3 m/s; In cyan: sphere with twice conductivity under wind velocity of 0.3 m/s. In green: reference sphere under wind velocity of 0.8 m/s; In yellow: sphere with twice conductivity under wind velocity of 0.8 m/s. (c), relative error between this frequency response and theoretical ISFDM frequency response. 3.3.3.2 Open Loop Time Response In addition to frequency response of transfer function obtained by means of DR, time response is also computed here. This is done by programing in Simulink R  the block diagram shown in g. 1. Fig. 26 shows that concordance between real outputs of ISFDM and computed outputs of 8-poles DR of this ISFDM is quite good as expected. 0 2000 4000 6000 8000 10000 Time [s] 0 2 4 6 8 10 12 ∆Temp [K] (a) Reference SFDM - PRBS time response Finite differences, 0.3 m/s wind Finite differences, 0.8 m/s wind 8 poles DR, 0.3 m/s wind 8 poles DR, 0.8 m/s wind 300 320 340 360 380 Time [s] 6 7 8 9 10 ∆Temp [K] (b) Reference SFDM - PRBS time response detail 0 2000 4000 6000 8000 10000 Time [s] 0 2 4 6 8 10 12 ∆Temp [K] (c) Higher conductivity SFDM - PRBS time response Finite differences, 0.3 m/s wind Finite differences, 0.8 m/s wind 8 poles DR, 0.3 m/s wind 8 poles DR, 0.8 m/s wind 300 320 340 360 380 Time [s] 6 7 8 9 ∆Temp [K] (d) Higher conductivity SFDM - PRBS time response detail Figure 26: Comparison between real outputs of ISFDM and computed outputs by means of 8-poles DR of this ISFDM. 48 3.3.3.3 Closed Loop Time Response Laboratory anemometer operates in constant temperature anemometer mode (CTA), i.e. a heater-sensor in contact with air ow is heated at constant temperature above the environment air temperature by means of control electronics. Thus, the power going into the heater is related to the velocity of the air. In order to do this, a sigma-delta topology called thermal sigma-delta modulation has been employed with excellent results [2,8,16]. In this approach, the closed feedback loop includes the heater-sensor as part of the sigma-delta converter topology (see g. 27). Figure 27: Thermal sigma-delta closed loop control Thermal domain in g. 27 is replaced by ISFDM to perform closed loop simulation using nite dierences, and by 8-poles DR of this ISFDM to perform closed loop simulation using diusive representation. In this way, the results using these two approaches can be compared. In g. 28 it is shown how this thermal domain box is replaced by the diusive representation model (g. 1). Initial conditions, eq. (7), are also included as impulse response. 49 Figure 28: Thermal sigma-delta closed loop control replacing thermal sphere by diusive representation. Inside the dashed box initial conditions are included as impulse response Closed loop simulation is performed keeping ∆ Temp = 8K, using 1kHz clock frequency. The power injected to the heater (equivalent control) is computed by averaging 50 consecutive sampling periods in the case of ISFDM and 5000 in the case of diusive representation. In g. 29, it is shown that closed loop time response matches quite well in the two cases. It is also shown that, equivalent control, once control surface is reached, changes more quickly in the case of reference sphere (b) than in the case of the higher conductivity sphere (d), although control surface ( ∆ Temp = 8K) is reached rst in the case of the reference sphere (a) than in the case of the higher conductivity sphere (c). The transitions of g. 30 are shown in more detail in g. 29. These simulations starts with initial conditions to 0. 50 0 50 100 150 200 250 300 Time [s] 0 2 4 6 8 10 ∆Temp [K] (a) Reference SFDM closed loop response SFDM, wind 0.3 m/s 8 poles DR, wind 0.3 m/s 50 100 150 200 250 300 Time [s] 0.02 0.025 0.03 0.035 0.04 Equivalent control [mW] (b) Reference SFDM closed loop response SFDM, wind 0.3 m/s 8 poles DR, wind 0.3 m/s 0 50 100 150 200 250 300 Time [s] 0 2 4 6 8 10 ∆Temp [K] (c) Higher conductivity - SFDM closed loop response SFDM, wind 0.3 m/s 8 poles DR, wind 0.3 m/s 50 100 150 200 250 300 Time [s] 0.02 0.025 0.03 0.035 0.04 Equivalent control [mW] (d) Higher conductivity - SFDM closed loop response SFDM, wind 0.3 m/s 8 poles DR, wind 0.3 m/s Figure 29: Closed loop response with initial conditions to 0. (a)(b): ∆ - temperature evolution under 0.3 m/s wind velocity, from reference sphere (a) and from the higher conductivity sphere in (b). (b)(d): evolution of the equivalent control under 0.3 m/s wind velocity, from reference sphere (a) and from the higher conductivity sphere in (d). In red, from ISFDM. In magenta and cyan from 8-poles DR. Equivalent control is computed by averaging 50 consecutive sampling periods in the case of ISFDM and 5000 in the case of diusive representation using 1kHz clock frequency 51 100 120 140 160 180 200 Time [s] 0.02 0.022 0.024 0.026 0.028 0.03 Equivalent control [mW] (b) Reference SFDM closed loop response detail SFDM, wind 0.3 m/s 8 poles DR, wind 0.3 m/s 100 110 120 130 140 Time [s] 7.5 7.6 7.7 7.8 7.9 8 8.1 ∆Temp [K] (c) Higher conductivity - SFDM closed loop response detail SFDM, wind 0.3 m/s 8 poles DR, wind 0.3 m/s 100 120 140 160 180 200 Time [s] 0.02 0.022 0.024 0.026 0.028 0.03 Equivalent control [mW] (d) Higher conductivity - SFDM closed loop response detail SFDM, wind 0.3 m/s 8 poles DR, wind 0.3 m/s 0 0.5 1 0 0.2 0.4 0.6 0.8 1 ∆Temp [K] 85 90 95 100 105 Time [s] 7.6 7.7 7.8 7.9 8 (a) Reference SFDM closed loop response detail SFDM, wind 0.3 m/s 8 poles DR, wind 0.3 m/s Figure 30: Transition details of closed loop response. (a)(b): ∆ temperature evolution under 0.3 m/s wind velocity, from reference sphere (a) and from the higher conductivity sphere in (b). (b)(d): evolution of the equivalent control under 0.3 m/s wind velocity, from reference sphere (a) and from the higher conductivity sphere in (d). In red, from ISFDM. In magenta and cyan from 8poles DR. Equivalent control is computed by averaging 50 consecutive sampling periods in the case of ISFDM and 5000 in the case of diusive representation using 1kHz clock frequency Fig. 31, also shows closed loop time response under 0.8 m/s wind velocity. In (b), in dashed green line, equivalent control is drawn from 8-poles DR of ISFDM under 0.8 m/s wind velocity and reference sphere, with initial conditions to 0. In red, equivalent control is drawn from ISFDM, from 0 to 300s under 0.3 m/s wind velocity, and from 300s to 500s under 0.8 m/s wind velocity. Wind velocity changes at 300s, and initial conditions at this point, for ISFDM are the temperatures in each elementary cell (g. 14). The same representation as in red is done in magenta for 8-poles DR, except that in this case initial conditions are unknown. The drawback here is that diusive representation provides a diagonalised version of thermal operators under each wind velocity, and therefore the values of the state variables are lost. Then, initial conditions should be obtained from ttings of experimental measures. However, surprisingly in this case, using state variables of the previous wind velocity as initial conditions, the curve is matched with the ISFDM curve. (b) also shows in green, the closed loop response from 8-poles DR but always under 0.8 m/s wind velocity, to check that close loop response tends to the same stationary value when the approximated initial conditions are used. In (d), the same is done, but in the case that ideal sphere has twice conductivity. 52 100 150 200 250 300 350 400 Time [s] 5.5 6 6.5 7 7.5 8 8.5 ∆Temp [K] (a) Reference SFDM closed loop response SFDM, wind 0.3 m/s, in 300s change to wind 0.8 m/s 8 poles DR, wind 0.3 m/s, in 300s change to wind 0.8 m/s 8 poles DR, wind 0.8 m/s 0 100 200 300 400 500 600 Time [s] 0.01 0.015 0.02 0.025 0.03 0.035 0.04 Equivalent control [mW] (b) Reference SFDM closed loop response SFDM, wind 0.3 m/s, in 300s change to wind 0.8 m/s 8 poles DR, wind 0.3 m/s, in 300s change to wind 0.8 m/s 8 poles DR, wind 0.8 m/s 100 150 200 250 300 350 Time [s] 6 6.5 7 7.5 8 8.5 ∆Temp [K] (c) Higher conductivity - SFDM closed loop response SFDM, wind 0.3 m/s, in 300s change to wind 0.8 m/s 8 poles DR, wind 0.3 m/s, in 300s change to wind 0.8 m/s 8 poles DR, wind 0.8 m/s 0 100 200 300 400 500 600 Time [s] 0.01 0.015 0.02 0.025 0.03 0.035 0.04 Equivalent control [mW] (d) Higher conductivity - SFDM closed loop response SFDM, wind 0.3 m/s, in 300s change to wind 0.8 m/s 8 poles DR, wind 0.3 m/s, in 300s change to wind 0.8 m/s 8 poles DR, wind 0.8 m/s Figure 31: Closed loop response from 0 to 300s, with wind velocity 0.3 m/s (except for green and yellow colours that starts with wind velocity 0.8 m/s) and initial conditions to 0. And closed loop response from 300s to 600s, with wind velocity 0,8 m/s and initial conditions from previous wind velocity state. (a)(b): ∆ temperature evolution under 0.3 m/s wind velocity from 0 to 300s and then wind velocity is changed to 0.8 m/s, from reference sphere (a) and from the higher conductivity sphere in (b). (b)(d): evolution of the equivalent control under 0.3 m/s wind velocity and then changed to 0.8 m/s wind velocity at 300s, from reference sphere (a) and from the higher conductivity sphere in (d). In red, from ISFDM. In magenta and cyan from 8-poles DR. In green and yellow from 8-poles DR, but always under 0.8 m/s wind velocity, to check that close loop response tends to the same stationary value when the approximated initial conditions are used. Equivalent control is computed by averaging 50 consecutive sampling periods in the case of ISFDM and 5000 in the case of diusive representation using 1kHz clock frequency 4 8-Poles Diusive Representation from Laboratory Data 4.1 8-Poles Diusive Symbol from Laboratory Data The laboratory thermal anemometer consists of a hollow silver sphere split in two hemispheres separated by a PCB board. Each hemisphere has a sensor/heater attached to its internal surface and the PCB board places a sensor/heater in the core of the sphere. A PRBS of 104 s and sample time 5·10−2 s is applied to one heater of one hemisphere of the laboratory sensor while it is under two dierent wind velocities, 0.3 m/s and 0.8 m/s. Therefore, the bandwidth where this system can be correctly identied is: 2π tK ,2π 2max{∆tk}=6.28 ·10−4,6.28 ·10+1 (58) 53 This PRBS is applied to one heater of one hemisphere when the other sensor/heaters are disconnected (other heaters OFF) and when the other heaters are kept to a higher constant temperature of ∆Temp =8.3K (other heaters ON). The conguration with all the heaters ON is that actually used when the anemometer is running. The diusive symbols obtained in these two cases are presented in g. 32. It can be seen in g. 32 that the dierences between DSs are at low frequencies, where convection phenomena are present, while, at high frequencies, where conduction phenomena are present, DSs are practically immutable. 10-2 100102 Frequency (rad· s -1) 0 100 200 300 400 500 600 700 η(ξ)k [K/J] (a) Diffusive symbols from laboratory data TLADS 8 Poles DR, wind 0.3 m/s TLADS 8 Poles DR, wind 0.8 m/s 10-3 10-2 10-1 100101 Frequency (rad· s -1) 0 20 40 60 80 100 η(ξk) [K/J· rad· s -1] (b) Interpolated diffusive symbols from laboratory data TLADS 8 Poles DR, wind 0.3 m/s TLADS 8 Poles DR, wind 0.8 m/s Figure 32: 8 diusive symbol displaced 20 times obtained from laboratory sensor data under two winds of 0.3 m/s (magenta and cyan) and 0.8 m/s (green and yellow), when the other heaters are OFF (magenta and green) and when the other heaters are ON (cyan and yellow). TLADs are plotted in dierent color 54 4.2 Frequency Response Using TLADSs obtained in previous section, 8-poles DR is performed and its frequency response drawn in g. 33. In these Bode plots, as in previous section from interpolated DSs (g. 32-(b)), variations at low frequencies are patent, whereas, at high frequencies, there is little variation. It is also patent that, these Bode plots resemble those seen before from ISTQM and ISFDM, in sec. 3.3. 10-3 10-2 10-1 100101 Frequency (rad· s -1) 25 30 35 40 45 50 Magnitude [K/W] (dB) (a) Laboratory thermal impedance 0.3 m/s, other heaters OFF 0.8 m/s, other heaters OFF 0.3 m/s, other heaters ON 0.8 m/s, other heaters ON 10-3 10-2 10-1 100101 Frequency (rad· s -1) -45 -40 -35 -30 -25 -20 -15 -10 -5 0 Phase (degree) (b) Laboratory thermal impedance 0.3 m/s, other heaters OFF 0.8 m/s, other heaters OFF 0.3 m/s, other heaters ON 0.8 m/s, other heaters ON Figure 33: Frequency response transfer function from 8 diusive representation symbols obtained from laboratory sensor data under two winds of 0.3 m/s (magenta and cyan) and 0.8 m/s (green and yellow), when the other heaters are OFF (magenta and green) and when the other heaters are ON (cyan and yellow) 55 4.3 Open Loop Time Response In order to verify that TLADS 8-poles DR identied in the previous sections describe correctly the laboratory system, outputs of the previous systems from DR are computed. This is done by using these TLADS 8-poles DR in the DR block diagram shown in g. 1 and coded in Simulink R  . These outputs are drawn along with the output from laboratory data in g. 34. It is seen that DR output and laboratory outputs are perfectly overlapped, which indicates that the laboratory system is well identied. 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 Time [s] 0 2 4 6 8 10 12 ∆Temp [K] (a) Prototype heaters OFF - PRBS time response Laboratory data 0.3 m/s wind DR 8 poles, 0.3 m/s wind Laboratory data 0.8 m/s wind DR 8 poles, 0.8 m/s wind 505 510 515 520 525 530 535 540 545 550 Time [s] 6.5 7 7.5 8 8.5 9 9.5 10 ∆Temp [K] (b) Prototype heaters OFF - PRBS time response detail 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 Time [s] 0 2 4 6 8 10 ∆Temp [K] (c) Prototype heaters ON - PRBS time response Laboratory data 0.3 m/s wind DR 8 poles, 0.3 m/s wind Laboratory data 0.8 m/s wind DR 8 poles, 0.8 m/s wind 500 510 520 530 540 550 Time [s] 5 6 7 8 9 ∆Temp [K] (d) Prototype heaters ON - PRBS time response detail Figure 34: Laboratory output time response. And output time response TLADS 8-poles DR obtained from laboratory sensor data under two winds of 0.3 m/s (magenta and cyan) and 0.8 m/s (green and yellow), when the other heaters are OFF (magenta and green) and when the other heaters are ON (cyan and yellow) 4.4 Closed Loop Time Response The systems shown in the previous sections were also simulated in the laboratory in the closed loop sigma-delta ( Σ−∆ ) conguration seen in sec. 3.3.3.3. These simulations were done by using Σ−∆ clock frequency 20 kHz, and the equivalent control power is computed by averaging power samples obtained from 1000 consecutive sampled periods. Fig. 35 shows the case when one sensor/heater works in a closed loop sigmadelta ( Σ−∆ ) conguration, under two dierent wind velocities, setting ∆ Temp=12K, the other heaters are disconnected, and zero initial conditions. In blue and red, the temperature of the heater/sensor and the equivalent power injected to this heater/sensor from laboratory data, and in green and magenta their respective temperature and the equivalent power from the TLADS 8-poles DR that was identied in the previous sections. The equivalent power in the case of the identied TLADS 8-poles DR was computed by using Σ−∆ clock 56 frequency 10 kHz, and by averaging 500 consecutive sampled periods. It can be seen that curves match quite well, which shows that the identied TLADS 8-poles DR can be used to model an anemometer prototype in a close loop conguration. 100 200 300 400 500 600 700 800 900 1000 Time [s] -2 0 2 4 6 8 10 12 14 16 18 ∆Temp [K] Closed loop reponse - other heaters OFF - ∆Temp=12ºK Laboratory data, wind 0.3 m/s Laboratory data, wind 0.8 m/s 8 poles DR, wind 0.3 m/s 8 poles DR, wind 0.8 m/s 0 100 200 300 400 500 600 700 800 900 Time [s] 0.035 0.036 0.037 0.038 0.039 0.04 0.041 0.042 0.043 0.044 0.045 Equivalent control [mW] Closed loop response - other heaters OFF - ∆Temp=12ºK Laboratory data, wind 0.3 m/s Laboratory data, wind 0.8 m/s 8 poles DR, wind 0.3 m/s 8 poles DR, wind 0.8 m/s Figure 35: Closed loop Σ−∆ when the other heaters are OFF. Comparison between laboratory data and TLADS 8-poles DR Fig. 36 shows the case when one sensor/heater works in a closed loop Σ−∆ conguration, under two dierent wind velocities, setting ∆ Temp=8.3K, while the other sensor/heaters are connected in the same Σ−∆ conguration from the beginning, and with zero initial conditions in the analysed sensor/heater. This conguration, with all the heaters ON, is the conguration that is actually used when the anemometer is running. In blue and red, the temperature of the heater/sensor and the equivalent power injected to this heater/sensor from laboratory data, and in yellow and blue their respective temperature and the equivalent power from the TLADS 8-poles DR that was identied in the previous sections. The equivalent power in the case of the identied TLADS 8-poles DR was computed by using Σ−∆ clock frequency 10 kHz, and by averaging 500 consecutive sampled periods. It can be seen again that curves match quite well, which shows that the identied TLADS 8-poles DR can be used to model an anemometer prototype in a close loop conguration. Figure 36: Closed loop Σ−∆ when the other heaters are ON. Comparison between laboratory data and TLADS 8-poles DR 57 10-3 10-2 10-1 100101 Frequency (ω) 0 50 100 150 200 250 300 η(ξk) wind 0.3 m/s wind 0.8 m/s X: 0.007616 Y: 160.6 X: 0.004971 Y: 243.4 Figure 39: Interpolated diusive symbols computed one pole being the rst intersection point frequency from table 9 surrounded by other geometrically scaled random poles until completing 11 poles Notice that the value of interpolated DS at delta pole for 0.8 m/s wind velocity is smaller than the value of interpolated DS at delta pole for 0.3 m/s wind velocity due to interpolation, since at high frequencies, interpolation reduction is bigger. In g. 40, where the diusive symbols are not interpolated, it is shown that the value at DS delta pole for 0.8 m/s wind velocity is actually a little higher than the value at DS delta pole for 0.3 m/s wind velocity. 10-3 10-2 10-1 100101 Frequency (ω) 0 5 10 15 20 η(ξ)k wind 0.3 m/s wind 0.8 m/s X: 0.007616 Y: 1.966 X: 0.004971 Y: 1.945 Figure 40: Diusive symbols computed using as pole the rst intersection point frequency from table 9 surrounded by other geometrically scaled random poles until completing 11 poles Fig. 41 shows that when the nearest poles are brought closer to the delta singularity, the value of the interpolated diusive symbol grows in the delta singularity whereas the other two nearest poles remain close to zero. That shows that this singularity is very localised compared to other values of the diusive symbol. Again, this growth is due to interpolation since the closer the poles are the higher the diusive symbol is. 64 10-3 10-2 10-1 100101 Frequency (ω) 0 500 1000 1500 2000 η(ξk) wind 0.3 m/s wind 0.8 m/s Figure 41: Interpolated diusive symbols computed using as pole the rst intersection point frequency from table 9 surrounded by two close poles and other geometrically scaled random poles until completing 11 poles Diusive symbols in g. 41 are shown in g. 42 without interpolation, becoming evident that the growth is due to interpolation. Moreover, it is also shown that the values at singularity points decrease (compared with the values at singularity points in g. 40) due to correlation with the nearest points which in this case is greater. 10-3 10-2 10-1 100101 Frequency (ω) 0 5 10 15 20 η(ξ)k wind 0.3 m/s wind 0.8 m/s X: 0.007616 Y: 1.705 X: 0.004971 Y: 1.819 Figure 42: Diusive symbols computed using as pole the rst intersection point frequency from table 9 surrounded by two close poles and other geometrically scaled random poles until completing 11 poles Previous results suggest a method to nd singularities from ideal data since the higher correlation from one isolated pole seems to aect directly the closest poles. This means singularity could be found moving three nearby geometrically spaced poles along the bandwidth where it is expected to be, until recovering a triangular shape with the value of DS at the base close to 0. This method is shown in g.43, where several diusive symbols without interpolation are plotted near singularity. Therefore, singularity is found near the central pole when it reaches a relative high value and the other two poles have similar values and these values are near zero (See green line in g.43). Fig. 43 shows the values of three nearby poles in dierent positions, matching singularity in green, and the other three colours displaced from singularity to the left. The central point positions are labelled. It can also be seen that, when singularities are not 65 matched, the singularity value is shared between the nearby poles according to the correlation of DR identication matrix. Moreover, the closest point to singularity is higher than the other. 2 3 4 5 6 7 8 9 10 11 12 Frequency (ω)×10-3 -0.5 0 0.5 1 1.5 2 2.5 η(ξ)k X: 0.007616 Y: 1.705 X: 0.0083 Y: 1.353 X: 0.0087 Y: 0.7194 X: 0.0093 Y: -0.7774 Figure 43: Method to search theoretical delta singularity with three nearby poles The previous method suggest a simplied method to nd singularities with only two nearby poles which would also reduce the correlation between them. This method is shown in gs. 44 and 45 for 0.8 m/s wind velocity, which has the isolated pole at 7.62·10−3 rad/s. The two poles are moved along bandwidth where it is expected to nd the isolated pole whereas the other poles remain in the same place, and the diusive symbol is computed. The isolated pole would be found where the two values of DS at these two poles change their sign. This is found at 7.64 ·10−3 rad/s, the dierence with the real position is for sure due to the correlation between those and the other poles. 10-3 10-2 10-1 100 Frequency (ω) -100 -50 0 50 100 150 η(ξ)k X: 4.288 Y: 21.11 X: 4.288 Y: -123.5 X: 0.006275 Y: -29.27 X: 0.008239 Y: -14.11 Figure 44: Method to search delta singularity with two nearby poles 66 10-2 100 Frequency (ω) 0 5 10 15 20 25 30 35 η(ξ)k (a) 7.6 7.62 7.64 7.66 7.68 Frequency (ω)×10-3 -5 0 5 10 η(ξ)k (b) X: 4.288 Y: 21.34 X: 4.288 Y: 21.36 X: 0.007651 Y: 0.5265 X: 0.007636 Y: -0.5222 X: 0.007644 Y: 1.455 Figure 45: Method to search delta singularity with two nearby poles. Pole found near 7.64 ·10−3 rad/s. Unfortunately, this method to nd isolated singularities is useless in the laboratory since there are not isolated deltas as in ideal case. This is because the convective heat transfer coecient, h , has a wide range of values depending on the wind incidence over the sphere surface. Moreover, the laboratory sphere, which is split in several sectors, is more complex than the ideal sphere. 67 Glossary Σ−∆ : Sigma-delta. CO 2 : Carbon dioxide. DR: Diusive representation. DS: Diusive Symbol. PCB: Printed circuit board. PRBS: Prseudorandom binary sequences. ISFDM: Ideal Sphere nite dierences model. ISTQM: Ideal Sphere thermal quadrupole model. TLADS: The lowest projection angle diusive symbol. 68