Inverse Estimation Problem of Determining the Unknown Timewise-varying Strength of a Primer Rapid Heat Source

Available online at www.sciencedirect.com

ScienceDirect Procedia Engineering 79 (2014) 295 – 304

37th National Conference on Theoretical and Applied Mechanics (37th NCTAM 2013) & The 1st International Conference on Mechanics (1st ICM)

Inverse Estimation Problem of Determining the Unknown Timewise-varying Strength of a Primer Rapid Heat Source Yang-Hsiung Ko*, Jiu-Zhang Lu, Kai-Tai Lu, Tsao-Fa Yeh CCIT, National Defense University, No.75, Shiyuan Rd., Daxi Township, Taoyuan County 33551, Taiwan (R.O.C.)

Abstract The research presents an input method to estimate the unknown time-varying strength on the sources of a primer applied at this on-line inverse heat conduction problem. We use the input method to recursively predict both the impulse strength acting on the front of thin copper foil and the transient temperature-time curves from the thin copper foil with the rapid combustion of different charge designs for primers. The input method includes the extended Kalman filter and weighted recursive least squares estimator. The purpose of the extended Kalman filter is to generate the recursive innovation sequence and the weighted recursive least squares estimator is used to identify the unknown time-varying strength. The results of simulation and experiment show that the input method can estimate the unknown time-varying heat flux on-line with more accuracy in the inverse heat conduction problem mode. The test experiments, namely different conductive thickness designs for primers with time-varying heat flux input, are solved to illustrate the effectiveness and good accuracy of the presented method. © 2014 2013 Elsevier The Authors. Published by Elsevier Ltd. under the CC BY-NC-ND license © Ltd. This is an open access article Selection and peer-review under responsibility of the National Tsing Hua University, Department of Power Mechanical (http://creativecommons.org/licenses/by-nc-nd/3.0/). Engineering. Selection and peer-review under responsibility of the National Tsing Hua University, Department of Power Mechanical Engineering Keywords: input method; primer; inverse heat conduction problem; Kalman filter

Nomenclature B H I

sensitivity matrix measurement matrix identity matrix

* Corresponding author. Tel.: +886-928891399; fax: +886-3-3906385. E-mail address: [email protected]

1877-7058 © 2014 Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/3.0/).

Selection and peer-review under responsibility of the National Tsing Hua University, Department of Power Mechanical Engineering doi:10.1016/j.proeng.2014.06.346

296

Yang-Hsiung Ko et al. / Procedia Engineering 79 (2014) 295 – 304

K Kalman gain Kb steady-state correction gain M sensitivity matrix N total number of spatial nodes P filter’s error covariance matrix Pb error covariance matrix q heat flux R measurement noise covariance rs measurement location s innovation covariance t time T temperature ǻt sampling interval v measurement noise vector w process noise vector Z measurement vector Greek symbols Į thermal diffusivity Ȗ forgetting factor ī input matrix į Dirac delta function ț thermal conductivity Ȟ residual of innovation sequence ȥ a vector of delay input and output ĭ state transition matrix ȍ coefficient matrix Subscripts ̝ estimated by filter ^ estimated T transpose of matrix 1. Instruction The mortars are used where high angles of fire are desired for plunging fire behind hills or into trenches, emplacements, or foxholes. The mortar primer (Fig.1) is designed to ignite the charge of ignition cartridge. It consists essentially of an arms primer to which a charge of black powder has been attached. The primer cap and the black powder charge are usually assembled in a metal tube. When used with the mortar, this tube is force fitted into the base of the cartridge case at the time of manufacture. Heat flow to the next charge occurs by radiation and conduction. The net radiation surrounds one body and the other in the case. Luminous flames which contain large numbers of solid particles in suspension radiate much more intensely than non-luminous flames. This accounts in large measure for the superiority of black powder for use in primers because the explosion products contain large amounts of solids such as potassium carbonate and sulfate, which radiate intense heat, in contrast to the nonluminous flames from smokeless powders.

Fig.1 A primer

Yang-Hsiung Ko et al. / Procedia Engineering 79 (2014) 295 – 304

Substantial theories of the inverse heat conduction problems have been developed and applied in recent years in the heat transfer field. Meng [1] presented a modified input estimation method for tracking a manoeuvring target. Scarpa and Milano [2] applied the Kalman smoothing technique in inverse heat transfer estimation and accurate results were obtained. Polat [3] used the wall function and the corrected presentation of the shearing stress to explain the influences on material conduction on impact surfaces. Linton and Agonafer [4] utilized the software applied to calculate the fluid dynamic problems to explore the heat-conducting phenomenon on a flat type plate. They compared the simulated results with results from actual experiments. Sathe [5] made use of the numerical simulation method to investigate the flow field and heat transfer by conduction on vertical square-column fins under the impact-emitting flow cooling process. Jonsson and Palm [6] found that the heat resistance value was estimated using the empirical formula from the actual experiment. Maveety and Jung [7] used several flow fields with different Reynolds number values to conduct numerical simulations and compare the results with experimental results. In the practical use, the parameters should be determined in real time. The input estimation method is composed of the Kalman filter and the recursive least square estimator. The amount of heat coming along with rapid primer combustion will expand rapidly. The temperature strength can be estimated when the ultrafast heat flux produced by the primer can be sufficiently measured in real time. The study proposes input estimation method to estimate the unknown ultrafast heat flux inversely using the measured temperature in the rear of a copper plate. The distribution of the estimated temperature curve can be obtained accurately. The input estimation method has the properties of rapid tracking and small affection owing to the noise. Most of input estimation researches were using numerical simulation with assumed known conditions. This study will verify the availability of the input estimation method by using the copper plate with 4 different thicknesses. These tests are ignited by using the standard black powder heat source from the primer. This study presents an online input method to estimate unknown ultrafast heat flux from the transient temperature data measured by a fine gage thermocouple on the outer surface of a copper plate. The influence owing to the processing noise variance Q is investigated to show the robustness and efficiency of the input estimation method. 2. Methods The standard black powder heat source with ultrafast heat flux of q(t ) is applied at x 0 . The copper plate is thin enough to assume that it is a homogeneous thermal conductor with thickness L. The boundary condition of x=L is insulated. The temperature is measured using a fine gage K type thermocouple. By equipping the thermocouple sensor at the position x L and measuring the surface temperature of the copper plate, the heat-conducting model is formed as shown in Fig.2.

Fig.2 Physical model of one dimensional inverse heat conduction problem(geometry and coordinates).

297

298

Yang-Hsiung Ko et al. / Procedia Engineering 79 (2014) 295 – 304

The 1-D heat-conducting governing equations are as follows: w 2T wT (1) k 2 Ucp 0  x  L ,0  t d tf wt wx wT k q (t ) x 0, t ! 0 (2) wx wT k 0 x L, t ! 0 (3) wx (4) T ( x, 0) T0 0 d x d L (5) z (t ) T ( x, t ) Q (t ) x L where T0 is the initial temperature. z (t ) is the temperature measurement. T(x,t) represents that the temperature is a function of time t and the position x . v(t) is the measurement noise which is assumed to be the Gaussian white noise with zero mean. The proportionality constant k is a transport property known as the thermal conductivity (W/mÂK) of a copper plate. Density ȡ and specific heat cp are two such properties used extensively in thermodynamic analysis. The product ȡcp (J/m3ÂK), commonly termed the volumetric heat capacity, measures the ability of a material to store thermal energy. The Equations (1)~(5) are the heat-conducting governing equations of the copper plate. The Equation (5) is the measured equation. The copper plate is separated into N-1 equal portions with the length of 'x . “i=1” is marked at x=0. The temperature is T1 . “ i N ” is marked at x L . The temperature is TN . By using the D`Souza [8] central differential method to substitute Equation (1) in the space derivative and the boundary condition equations, the heat-conducting governing equations can be transferred as following differential equation. The deduced process can be shown as follows: wTi wx when i 1 x

T1 (t )

Ti 1  Ti 1 w 2T , 2 2 'x wx

wT1 wt

k w 2T1 U c p wx 2

Ti 1  2Ti  Ti 1 'x 2

k T2  2T1  T0 'x 2

Ucp

From the boundary condition in Equation (2), it is clear that T T wT k 1 k 2 0 q (t ) wx 2'x 2q (t )'x T0 T2  k By substituting Equation (9) to (7) can be rearranged as follows: x 2k (T2  T1 ) 2q (t ) T1 (t )  U c p 'x 2 U c p 'x when i

2,3,! , i,! , N can be obtained as the following equation. k Ti 1  2Ti  Ti 1 Ti (t ) ( ) 'x 2 Ucp x

when i

2,3," N , the Equation (11) can be rearranged as follows: k (T  2Ti  Ti 1 ) Ti (t ) U c p 'x 2 i 1

when i=N, the Equation (11) can be rearranged as follows: k TN (t ) (TN 1  2TN  TN 1 ) U c p 'x 2 From the boundary condition in Equation (5), it is clear that

(6)

(7)

(8) (9)

(10)

(11)

(12)

(13)

299

Yang-Hsiung Ko et al. / Procedia Engineering 79 (2014) 295 – 304

wTN T T k N 1 N 1 0 2'x wx TN 1 TN 1 By substituting Equation (15) to (13) can be rearranged as follows: 2k 2k TN (t ) T  T U c p 'x 2 N 1 U c p 'x 2 N k

(14) (15) (16)

From the Equations (10-12, 15) and the fictitious process noise input, the one-dimensional continuous-time state equation can be obtained as the following: T (t ) AT (t )  Bq(t )  Cw(t ) (17) T (t )

Let D

A

B

T

^T1 (t ), T2 (t ),! , TN (t )` k

, then, the state matrix, A is shown as follows:

U c p 'x 2 ª 2D « D « « 0 « « # « 0 « « # « « « # « « 0 «¬ 0

(18)

2D 2D

D

0

D 2D %

! " " " 0 0

" 0

!

2D

0

0

0

D

2D

0

0

D

2D

D 0 0

D

"

%

" "

"

0 " 0

" 0

"

0 " 0 " 2D 0

0

D

%

%

%

%

0

"

0

2

"0 º " 0 »» "0 » » » "0 » » "0 » "0 » » "0 » » "0 » 2D »¼

(19)

ª 2 º « U c 'x » « p » « 0 » « » « 0 » « # » « » «¬ 0 »¼

(20)

where A is the state matrix, B and C are the input matrices. w  t is assumed to be the Gaussian white noise with zero mean and it represents the modeling error. The continuous-time state Equation (17) can be discretized with the The discrete-time state equation and its relative equations are shown as follows: sampling time 't .

!

(21)

T (k  1) ĭ[(k  1)'t , k 't ]T (k )  *(k  1)q(k )  w(k  1)

where ĭ[( k  1) 't , k 't ]

*(k  1) w(k  1)

³ ³

( k 1) 't

k 't ( k 1) 't

k 't

e A' t

ĭ[(k  1)'t ,W ]B(W )dW ĭ[(k  1)'t ,W ]G (W ) w(W )dW

w  k is the processing error input vector assumed to be the Gaussian white noise with zero mean and with the

variance E{w(k ) wT ( j )} QG kj . G kj is the Dirac delta function. The discrete-time measurement equation is shown below. Z (k ) HT (k )  v(k ) where Z( k ) is the observation vector at the kth sampling time. H

>0 0.......1@

(22) is the measurement matrix. v(k ) is the

300

Yang-Hsiung Ko et al. / Procedia Engineering 79 (2014) 295 – 304

measurement error vector assumed to be the Gaussian white noise with zero mean and with the variance E{v(k )vT ( j )} RG kj . The input estimation method has two parts: one is the Kalman filter without the input term and the other is the fuzzy weighted recursive least square estimator. The system input is the unknown ultrafast heat flux. The Kalman filter is operating under the setting of the processing error variance Q and the measurement error variance R. It is to use the difference between the measurements and the estimated values of the system temperature as the functional index. Furthermore, by using the recursive least square estimator, the ultrafast heat flux can be precisely estimated. The equations of the Kalman filter are shown as follows: (23) X  k / k  1 )X  k  1/ k  1 P (k / k  1)

s (k ) K (k )

)P(k  1/ k  1))T  Q

(24)

HP(k / k  1) H T  R

(25)

P (k / k  1) H T s 1 (k )

(26)

P (k / k ) [ I  K  k H ]P(k / k  1)

(27)

Z (k )

(28)

Z (k )  HX (k / k  1)

X (k / k )

X (k / k  1)  K (k ) Z (k )

(29)

The recursive least square estimator: B (k ) H [)M (k  1)  I ]* M (k ) [ I  K (k ) H ][)M (k  1)  I ] 1

T

1

T

(30) (31) 1

Kb (k) J Pb (k 1)B (k) ª¬B(k)J Pb (k 1)B (k)  s(k)º¼

Pb (k )

> I  Kb (k ) B(k )@ J 1 Pb (k  1)

(32) (33) (34)

qˆ (k ) = qˆ (k  1) + K b (k ) [ Z (k )  B(k ) qˆ (k  1) ] Xˆ  k / k  1 = )Xˆ  k  1/ k  1  *qˆ (k )

(35)

Xˆ  k / k = Xˆ  k / k  1 + K  k [ Z  k - HXˆ  k / k  1 ]

(36)

where P is the filtering error covariance matrix. s(k) is the residual covariance. Z(k) is the bias innovation produced by the measurement noise and the input disturbance. B(k) and M(k) are the sensitivity matrices. Kb(k) is the correction gain. Pb(k) is the error covariance of the estimated input vector. J is the weighting constant or weighting factor in the range 0  J d 1 . qˆ(k ) is the estimated input vector. After the state equation is obtained, the inverse estimation process is carried out using the input estimation method system based on the Kalman filter and the recursive least squares estimation algorithm. 3. Experimental methods

The purpose of experiment is to inversely estimate the temperature and ultrafast heat flux of the primer by using the temperature measurements of the copper plate surface. The entire experiment modular includes the impulse heat source, the sensor, the high-speed data acquisition device, and the computer module. The copper plate is heated by adopting the standard impulse heat source and the fine gage K type thermocouple is set in the rear of the thin copper plate. The DAQ is linked with the thermocouple to measure the temperature data. In the ultrafast heat flux for explosives, the measurement method must cope with rapidly changing burning chiefly with the transient temperature measurements. The experiment devices and the elements used are illustrated as follows (Fig.3).

Yang-Hsiung Ko et al. / Procedia Engineering 79 (2014) 295 – 304

301

Fig.3 Schematic diagram of the experimental facility

We used fine gage K type thermocouples whenever fast, accurate temperature measurements were required. The fine wire diameters enable accurate temperature measurements without disturbing the base temperature of the body by keeping heat transfer via the leads to a minimum. They are available in wire sized 0.125 mm in diameter. Intel processor 1.6G computer, the signal software, and the Matlab programming language can be used to process the signal data. The input estimation method can be programmed by using the Matlab programming language. The temperature measurements are then regarded as the inputs into the input estimation method which is to estimate the ultrafast heat flux in the front of the copper plate. 4. Results and discussions

To verify the estimated performance of the input estimation method, an inverse heat conduction problem of the copper plate is considered. The ultrafast heat flux in the front and temperature are estimated inversely by measuring the temperature in the rear. The test is heated by the standard impulse heat source with the fixed power (three black powder pellets: KNO3, 74% ; Sulfur, 10.4%; Charcoal, 15.6% ). The rear of the copper plate is insulated with the insulating rubber. The transient temperature data are measured using a fine gage K type thermocouple in the rear of the copper plate whose 4 different thicknesses are 0.15 mm, 0,20 mm, 0.25 mm, and 0.30 mm. Using a total time period of t f =1 sec and a sampling interval of 't =0.0005 sec (20 kHz), find the measurement temperature curves of 4 different tests shown in Fig.4. It shows that the measurement error of the thermocouple is approximately ±0.01% (with the measurement noise covariance, R=10-4). Using the space interval of ǻx=L/10, the processing noise covariance matrix of Q= 103 and the measured temperatures in the rear of different samples with 4 different thicknesses (L=0.15, 0.20, 0.25 and 0.30 mm), ultrafast heat flux vs. time (Fig.5(a)~(d)) may be obtained from Fig.4. Fig.5(a)~(d) look very much alike and are in good agreement. Therefore, it appears that the estimated ultrafast heat flux is the heat flux generated by the primer on the copper plate front. 200 180 Copper Thickness & Max. Temp.

o

Measured Temperature, C

160

o

0.15 mm & 159.5 C o

0.20 mm & 149.6 C

140

o 0.25 mm & 142.5 C

120

o 0.30 mm & 137.1 C

100 80 60 40 20 0 0.0

0.2

0.4

0.6

Time, sec

0.8

1.0

302

Yang-Hsiung Ko et al. / Procedia Engineering 79 (2014) 295 – 304 Fig.4 Experimentally measured temperature- time curves at x=L for properly ignited three pellets of black powder.

20000

20000 20000

2

16000

15000

14000 12000 10000

10000

8000 6000 4000 2000

5000

0 0.00

0.05

0.10

0.15

0.20

0 0.2

0.4 0.6 Time, sec

0.8

14000 12000 10000

10000

8000 6000 4000 2000

5000

0 0.00

1.0

0.0

(a) The estimated ultrafast heat flux (L=0.15 mm).

0.05

0.10

0.15

0.20

0.2

0.4 0.6 Time, sec

0.8

1.0

(b) The estimated ultrafast heat flux (L=0.20 mm).

20000

20000 20000

20000

t=0.012 sec Max. SHF =19194

18000

2

16000

15000

14000 12000 10000

10000

8000 6000 4000 2000

5000

0 0.00

0.05

0.10

0.15

t=0.012 sec Max. SHF =19101

18000

0.20

0

Estimated SHF, kW/m

2

16000

15000

0 0.0

Estimated SHF, kW/m

t=0.0125 sec Max. SHF =19439

18000

Estimated SHF, kW/m

2

Estimated SHF, kW/m

20000

t=0.0125 sec Max. SHF =19324

18000

16000

15000

14000 12000 10000

10000

8000 6000 4000 2000

5000

0 0.00

0.05

0.10

0.15

0.20

0 0.0

0.2

0.4 0.6 Time, sec

0.8

(c) The estimated ultrafast heat flux (L=0.25 mm).

1.0

0.0

0.2

0.4 0.6 Time, sec

0.8

1.0

(d) The estimated ultrafast heat flux (L=0.30 mm).

Fig.5 Estimated ultrafast heat flux–time curve at x=0 by input estimation method with measured temperature in Fig. 4 as input.

Fig.6 shows the increasing modeling error as Q=10-1, 103 , and 106 . The Kalman filter is operating under the setting of the processing error variance Q and the measurement error variance R. It regards the renewal values produced by using the difference between the temperature estimates and the temperature measurements as the functional index, and utilizes the real-time least square algorithm to precisely estimate the ultrafast heat flux. The measurement error variance of the thermocouple is assumed as a given value R 104 . The processing error variance Q is adjusted between 101 to 106 . The temperature measurements are then regarded as the inputs into the input estimation method which estimates the ultrafast heat flux in the rear of the copper plate, L=0.15 mm. When the processing noise variance Q increases, it will accelerate the estimation speed and produce the better estimation result. When the processing noise variance Q decreases, the error covariance matrix will decrease and the Kalman gain

303

Yang-Hsiung Ko et al. / Procedia Engineering 79 (2014) 295 – 304

K (k ) will decrease. The main reason is that the correction need is decreasing and the smaller Kalman gain K (k ) is needed to offer that correction. When the processing noise variance Q increases, the error covariance matrix will increase. The main reason is that the correction need is increasing and the larger Kalman gain K (k ) is needed to offer that correction. As a result, the Kalman filter will have faster correction performance but have larger oscillation in the estimation process. When the processing noise variance Q increases ( Q 101 ~ 103 ), the input estimation method makes the new measurement significant for the correction of estimation. The relative root mean square error is reduced. If the modeling error is too large ( Q 106 ), the Kalman filter will produce larger oscillations which cause the increase of the relative root mean square error value. As a result, the relatively smaller value of relative root mean square error is making the parameter Q 103 the most excellent option. Q 101 , 103 , and 106 in Fig.6 are used to show that when the error is larger, the result with the larger oscillation and faster response will be produced. The processing noise variance of Q 101 in Fig.7 shows that the estimated temperature curve can’t be estimated accurately. 20000 Estimated SHF -1

Q=10

3

Q=10

15000

6

Estimated SHF, kW/m

2

Q=10

10000

5000

0

0.0

0.2

0.4 Time, sec

0.6

0.8

1.0

Fig.6 Estimated ultrafast heat flux curves at x=0 by means of input estimation method with the processing noise variance of 103 , and 106 as input.

Q 101 ,

200 Process noise variance Q

180

Measured Temperature -1 Estimated Temp. (Q=10 )

o Estimated Temperature, C

160

-3 Estimated Temp. (Q=10 ) 6 Estimated Temp. (Q=10 )

140 120 100 80 60 40 20 0 0.0

0.2

0.4

0.6

Time, sec

0.8

1.0

Fig.7 Estimated temperature-time curve at x=L by direct heat conduction problem with estimated ultrafast heat flux in Fig.6 as input.

304

Yang-Hsiung Ko et al. / Procedia Engineering 79 (2014) 295 – 304

5. Conclusions

The input estimation method is utilizing the transient measured temperature data to estimate the ultrafast heat flux in the front of the copper plate at the end of the primer. The influence on the estimation caused by the processing noise variance Q is investigated by utilizing the experimental verification. The results reveal that if the processing noise variance Q increases, the faster estimation convergence and larger oscillations will be produced. The experimental verification shows that the input estimation method has the properties of better targat tracking capability and more effective noise reduction, and that it is an efficient inverse estimation method for the estimation of the unknown ultrafast heat flux of the primer. References [1] H.H. Meng, Aircraft Maneuver Detection Using an Adaptive Kalman Filter, M. Thesis, Naval Postgraduate School, Mon. CA, USA. 1989. [2] F. Scarpa, G. Milano, Kalman smoothing technique applied to the inverse heat conduction problem, Num. Heat Trans. B, 28 (1995) 79-96. [3] S. Polat, A.S. Mujumdar, A.R.P, Heiningen, W.J.M. Douglas, Numerical model for turbulent jet impinging on a surface with through-flow, Int. J. Thermophysics. 5, No.2 (1991) 172-180. [4] R.L. Linton, D. Agonafer, Coarse and detailed CFD modelling of a finned heat sink, IEEE Trans. on Comp. Pack. and Manu. Tech. A, 18, No.3 (1995) 517-520. [5] S. Sathe, K.M.K. Kelkar, C. Karki, C. Tai, C. Lamb, S.V. Patankar, Numerical prediction of flow and heat transfer in an impingement heat sink, J. of Elect. Pack. 119, No.1 (1997) 58-63. [6] H. Jonsson, B. Palm, Thermal and hydraulic behaviour of plate fin and strip fin heat sinks under varying bypass conditions, IEEE Trans. on Comp. Pack. and Manu. Tech. 23, No.1 (2000) 47-54. [7] J.G. Maveety, H.H. Jung, Design of an optimal pin-fin heat sink with air impingement cooling, Int. Comm. on Heat and Mass Tran. 27, No.2 (2000) 229-240. [8] N. D’Souza, Numerical Solution of One-dimensional Inverse Transient Heat Conduction by Finite Difference Method, ASME 75-WA/HT81, 1975.