Microsoft Word - 476hernandez.docx CHEMICAL ENGINEERING TRANSACTIONS VOL. 43, 2015 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Chief Editors: Sauro Pierucci, Jiří J. Klemeš Copyright © 2015, AIDIC Servizi S.r.l., ISBN 978-88-95608-34-1; ISSN 2283-9216 A General Method for the Solution of Inverse Problems in Transport Phenomena Marco Vocciantea, Andrea P. Reverberib, Vincenzo G. Dovì*b aDipartimento di Ingegneria Civile, Chimica e Ambientale, Università di Genova, Via Montallegro 1, I-16145 Genova bDipartimento di Chimica e Chimica Industriale, Università di Genova, Via Dodecaneso 31, I-16146 Genova dovi@istic.unige.it The typical inverse problems in transport phenomena are given by partial differential equations with unknown boundary conditions, which are to be estimated from measurements corresponding to solutions of the PDEs or of their gradients. The resulting problem is an ill-posed problem, which can’t be solved unless it is adequately regularized, because arbitrarily small errors on the data can give rise to very large deviations in the reconstruction of the unknown boundary conditions. In other words the estimated solution does not depend continuously on the data. The method proposed in this paper is the generalization of methods already applied to a number of problems (diffusion, heat transfer, percolation). The unknown boundary condition is replaced by a piecewise constant (or linear) functions with unknown coefficients. This approximation makes it possible to solve the resulting equation analytically and to estimate the coefficients by comparing theoretical and experimental values by means of linear least squares methods. The ill-posedness of this class of problems (resulting in singular matrices of the least squares problems) can be tackled using the well-known method proposed by Tikhonov. The value of the parameter contained in Tikhonov’s method is estimated from the statistical noise on data by means of the Morozov’s discrepancy principle. If additional information on the structure of the solution is available, specialised algorithms can be developed on a case-by-case basis. The presence of additional unknown non-linear parameters in the partial differential equations (such as diffusion coefficients, conductivities or adsorption coefficients of plant roots) leads to a particular model known as separable least squares. This procedure makes use of the conditional optimality principle for carrying out the overall identification problem in terms of the non-linear parameters only. 1. Introduction Conservation laws in continuous media take the general form ( )Q Q q t ∂ + ∇ ⋅ + = ∂ Φ v (1) where Q is the quantity being conserved, Φ the flux of Q in the absence of fluid transport, v the bulk velocity and q a source or sink of Q . Conservation laws are generally coupled with constitutive equations, which provide the necessary analytical expressions for Φ and q (Bagatin et al., 2014). As far as mass and energy conservation is concerned, the well-known Fick’s (Fourier’s) law is generally used for Φ , i.e. k Q= − ∇Φ . Furthermore, the velocity vector is often determined from separately solved fluid dynamics problems. The conservation of momentum uses Navier-Stokes equations (with laminar and turbulent coefficients) for Φ with Q ρ= v . For the sake of illustration simplicity only mass and energy conservation are considered in this article. This makes it possible to rewrite Eq. (1) as: DOI: 10.3303/CET1543270 Please cite this article as: Vocciante M., Reverberi A., Dovi V., 2015, A general method for the solution of inverse problems in transport phenomena, Chemical Engineering Transactions, 43, 1615-1620 DOI: 10.3303/CET1543270 1615 ( )| ,Q Q t t ∂ = ∂ L x | θ (2) where θ is the set of unknown parameters contained either in Φ or q and L is a differential operator of the second order. Equation (2) requires the definition of one initial condition and two boundary conditions: ( ) ( ) ( ) ( ) ( ) 0 (1) 1 ( 2) 2 0, | , 0 | , 0 Q Q f Q t Q t = = = = I B B x x x | x | θ θ (3) In inverse problem the exact functional form of one of the two boundary conditions has to be estimated from experimentally measured data. It will be supposed (without any loss of generality) that the unknown boundary condition is the one related to the boundary ( 2 )x , i.e. ( ) ( ) ( ) ( ) 0 (1) 1 0, | , 0 - Q f Q t dQ k h Q f t dx = =B + =(2) (2) x= x x= x x x x | θ (4) Inverse problems generally arise from the necessity of estimating either the very unknown boundary conditions (Solisio et al., 2012) or single physical-chemical parameters in experiments in which boundary conditions are not under the control of the experimenter (Reverberi et al. 2013). Occasionally, both boundary conditions and non-linear parameters are estimated from the data (Palazzi et al., 2014). In all these cases the geometry of the problem is known and is generally selected so as to take advantage of symmetry factors (such as thin disks, cylindrical rods, spherical shells). Thus, the spatial coordinates x reduce to the scalar variable x in the case of small rectangular slabs, to the variable r in spherical shells and long circular conduits. Furthermore the experimental conditions over small intervals of the dependent variables (concentrations or temperature) make it possible to assume a constant coefficient for the Fick’s (Fourier’s) coefficient, so that L can be assumed to be a linear differential operator of the second order. The linearity of L makes it possible to apply Duhamel’s theorem for the solution of Eq. (1) subject to conditions (4). In its simplest formulation Duhamel’s principle can be stated as follows (Carslaw and Jaeger, 1959): suppose ( )* ,Q x t λ| |θ is the solution of problem (1)-(4) when ( )f t is replaced by ( )f λ , i.e. a time independent parameter. The complete solution ( ),Q t x λ| |θ is then given by: ( ) ( ) ( )*0 0, , t Q x t f x Q x t dλ λ λ= + −| | |θ θ (5) The general procedure for the estimation of the unknown boundary condition, which is shown in the next section, is based on a suitable parameterisation of ( )f t and on the use of Eq. (5). 2. The general procedure for the estimation of the unknown boundary condition The assumptions made (linearity, one dimensionality, and time-independence of the unknown boundary condition) make it possible in many cases to evaluate an analytical solution for ( )* ,Q x t λ| |θ . A suitable procedure for the evaluation of the integral contained in Eq. (5) is to parametrise the function ( )f λ as a piecewise polynomial function whose coefficients satisfy the usual smoothness conditions at the boundary of each interval into which the overall range is divided. If M equal intervals of width h=T/M are assumed, ( )f λ can be rewritten as: 1 ( ) ( ) M i i i f c fλ λ = =  (6) 1616 where fi(λ) are the polynomials and ci are adjustable coefficients to be determined by regression. If the fi(λ) are j-th order polynomials, (j-1) continuity/smoothness conditions provide the additional conditions 1 1( ) ( ) q q i i i iq q d d f f d d λ λ λ λ+ + = at each interval boundary. (7) The simplest (and generally adequate) assumptions are to consider either a piecewise constant function (in which case no continuity condition is required) or a piecewise linear function. In the latter case introducing the continuity condition provides the general relation: 1 1 1 1 ( ) (0) ( ) , m i m m m m i f f h c cλ λ λ λ λ λ − − − = = + + − ≤ ≤ (8) where mc is the angular coefficient of the linear approximation in the range 1m mλ λ λ− ≤ ≤ . Since the functional form of ( )f t is reconstructed from comparison with measured data, a physical dependence of f on Q is also allowed, i.e. ( ) ( )( ),f t Q t tψ= . As shown in previous articles (Dovì et al., 1994) the error introduced by this approximation can be made arbitrarily small by increasing the number of intervals. Thus, not only do traditional numerical schemes lack a general strategy to deal with unknown boundary conditions and to regularize the solution, but they also aren’t substantially more accurate. The overall procedure is exemplified for three important special cases in the next section. 3. Important special cases a) One dimensional time-dependent heat/mass diffusion between two planes at a distance L 2 2 T T K t x ∂ ∂ = ∂ ∂ 0( , 0)T x T= /2 0 x L T x = ∂  = ∂  ( )( , / 2) ,T t L f T t= The analytical solution corresponding to a piecewise constant approximation of ( ),f T t is given by (Dovì et al. 1994) ( ) ( ) ( ) ( ) ( ) ( ) ( ){ } ( ) ( ) ( ) ( ) 2 2 1 2 1 1 42 1 2 1 1 1 1 1 4 2 2 42 1 1 4 2 1 / 2, / 2 1 32 / 2 1 , 96 4 i i r r n t t n t t r r i i i i n n t t r r n T L t t t c t t e e n L K c t t e n = K L β β β α π π α α β π − − − ∞ − + − − + − − − = = ∞ − + − − =  ≤ ≤ = − + − + +      − + + − =     b) One dimensional time-dependent heat/mass diffusion in cylinders of infinite length or (equivalently) steady state heat/mass diffusion in cylinders of finite length. The velocity distribution ( )u r can follow either a laminar or a turbulent distribution law. 1 ( )p T T C u r K r x r r r ρ ∂ ∂ ∂ =  ∂ ∂ ∂  (0, ) 1T r = 0 0 r T r = ∂  = ∂  [ ]( , ) r R T K H T x R r = ∂ − = ∂  1617 The analytical solution corresponding to a piecewise linear approximation of ( )f t is given by ( , ) ( , ) ( , )i i i T r x r x c r xω ξ= +  , where the analytical expression of the functions ( , )r xω and ( , )i r xξ for both the laminar and turbulent case can be found in (Reverberi et al, 2013). c) One dimensional percolation through an unsaturated homogeneous soil under a time-dependent surface flux. For practical reasons (Yuan and Lu, 2005), the problem is modelled in terms of capillary pressure ( ) = − ( ) + 1 ( ,0) = ( ) + 1 = − ( ) (0, ) = ( , ) It can be shown that the analytical solution corresponding to a piecewise constant approximation of ( , ) is given by (z, t < < t ) = (z,t), ( , ) = ( ) + ( )( ) + ( , ) + ( − )ℎ ( ) where is the hydraulic conductivity at saturation, the soil pore-size distribution parameter and ( ), ( ), ( ) and ℎ ( ) are known function of z an t. 4. Ill-posedness and regularisation techniques As shown in the previous section experimental measurements are generally predicted by the sum of M terms dependent on the unknown nonlinear parameters, each multiplied by the coefficient mc that defines the unknown boundary condition. If the nonlinear parameters are known, the determination of the coefficients reduces to ordinary linear least squares, i.e. they can be determined by the minimization of ( ) 2 2 k i i k k i x,t y c G x ,t  ′ = −      k y - G( )c =min (9) where ( )i i k k i c G x ,t is the value of the measured value ky predicted by the solution (5) when Eq. (6) is used for ( )f λ . The least squares solution of (9) is given by ( ) 1− ′T Tc = G G G y (10) provided the symmetric matrix ( )TG G is full rank. However, due to the ill-posed character of inverse problems of PDEs, (Kabanikhin, 2008), full rankness cannot be assumed for ( )TG G . Thus the determination of c requires the use of adequate regularization techniques for preventing arbitrarily small errors on the data from giving rise to very large deviations in the reconstruction of the vector c . This can be accomplished by replacing ( ) 1− ′T TG G G y by + ′G y , where +G is the pseudo-inverse of G . A precise formulation of +G can be established using the well-known method proposed by Tikhonov and Arsenine (1976). The value of the parameter contained in Tikhonov’s method is estimated from the statistical noise on data by means of the Morozov’s discrepancy principle (Morozov, 1984), as shown in section 4.1. If additional information on the structure of the solution is available, specialised algorithms can be developed on a case-by-case basis. An example is considered in section 4.2. 1618 4.1 Tikhonov-Morozov method The matrix ( )G H+ is approximated by ( ) 1T T TG G G−+Ψ Ψ where the matrix Ψ is selected so that ( )T TG G +Ψ Ψ is invertible and the vector c depends continuously on the data. This is equivalent to replacing the original minimisation problem 2 x,t′y - G( )c =min with the problem 2 2 x,t′ +y - G( )c cΨ =min (11) A frequent choice for Ψ is Ψ =αI, where the parameter α is selected so as to satisfy some statistical properties. In Morozov’s discrepancy method α is iteratively increased until 2x,t′y - G( )c is equal to the variance of the data multiplied by the number of data points. Thus, the solution provided, in addition to limiting the norm of c through the introduction of the matrix Ψ , is consistent with the statistical properties of the data. Additionally, some smoothing can be introduced into the solution by introducing tridiagonal and pentadiagonal matrices (Twomey, 1963), although this heuristical method does not rest on a sound statistical basis. 1 1 0 0 1 2 1 0 0 1 2 1 0 0 2 5 4 1 0 0 0 1 2 1 0 0 1 4 6 4 1 0 0 1 4 6 4 1 0 0 0 1 2 1 0 0 0 1 2 1 0 1 4 5 2 0 0 1 1 0 1 2 1 − −       − − −       − − −    = + + − −       −     − − −       − −            IΨ α ε ε The tridiagonal matrix tends to reduce the values of the slope, whereas the pentadiagonal matrix tends to hamper oscillations. The value of ε is adjusted until the solution is sufficiently smooth. 4.2 Unimodal regression Sometimes additional constraints dependent on the particular conditions of the experimentation can be added to the estimation task. A frequent condition is the unimodal distribution of the vector c . In other words, there exists an index i, such that k-1 k k k+1c c for k < i c c for k i ≤ ≥ ≥ (or analogously k -1 k k k+1c c for k < i c c for k i≥ ≤ ≥ ). Clearly if i=1 or i=M this condition implies a monotonous distribution of c . The minimization of Eq. (9) is subject now to a set of linear constraints, which change the optimization task into a QP (quadratic programming) problem. If, as is often the case, the value of i is not known, the QP problem has to be repeated for all i ( 1 i M≤ ≤ ) until the value of Eq. (9) is consistent with Morozov’s principle. 5. Estimation of non-linear parameters If the non-linear parameters θ are unknown the minimisation problem (9) changes to ( ) 2 2 k i i k k i x,t y c G x ,t  ′ = −      k y - G( )cθ θ =min (12) While in principle the problem can be solved by using general non-linear minimization algorithms with respect to both θ and c , the high dimensionality of the resulting task, the difficulty of determining suitable starting values for c and the possibility of multiple local minima reduce the reliability of the solutions determined. However, these difficulties can be overcome if the principle of relative optimality is applied, which in our case 1619 provides optimal values of the vector c for every value of the non-linear parameter θ , i.e. using separable least squares (Golub and Pereyra, 1973). For every value of θ the vector c is computed as ( )G+ ′=c yθ and the computation of the non-linear parameter can be carried out by minimising the objective function { } 22 2 ( ) ( ) ( ) ( ) ( )+ + ⊥′ ′ ′ ′− = − = Gy G G y I G G y P yθ θ θ θ θ (13) where ( ) ( )⊥ += −GP I G Gθ θ is the projector on the orthogonal complement of the column space of G with respect to the non-linear parameters only. The typical Gauss-Newton step is given by ( )n+1 n n n nτ D ( ) ( ) + ⊥ ⊥ ′ ′= −  G GP y P yθ θ θ θ (14) where the operator D is the Fréchet derivative of matrices and nτ is an interpolation-extrapolation factor used to locate the minimum of the objective function in the direction ( )n nD ( ) ( ) + ⊥ ⊥ ′ ′ G GP y P yθ θ . Furthermore ( )2 nD ( ) + ⊥ ′ GP yσ θ (where 2σ is the variance of the data) provides an estimate of the variance of the parameter θ , once convergence at θ * has been attained. 6. Conclusions The general method presented is capable of solving a large number of estimation problems related to transport phenomena when the boundary conditions are only partially known. This situation is frequently encountered in real life problems and is generally dealt with using approximate models or ad-hoc assumptions. Using more realistic models leads to the estimation of more reliable boundary profiles in critical situations (De Rademaeker et al., 2014) and provides more accurate estimations of physical- chemical parameters which can greatly add up to the quality of design of plants and processes. Acknowledgements Financial support by the EU FP7 project EFENIS is gratefully acknowledged. References Bagatin R., Klemeš J.J., Reverberi A.P., Huisingh D., 2014, Conservation and improvements in water resource management: a global challenge, J Cleaner Prod, 77, 1-9. Bro, R., Sidiropoulos, N.D., 1998, Least Squares Algorithms under Unimodality and Non-negativity Constraints, J. Chemometrics, 12, 223-247. Carslaw H.S., Jaeger, J.C., 1959, Conduction of heat in solids, 2nd Editions, Clarendon Press, Oxford, UK. De Rademaeker E., Suter G., Pasman H.J., Fabiano B., 2014, A review of the past, present and future of the European loss prevention and safety promotion in the process industries, Proc. Saf Environ, 92, 280-291. Dovì, V.G., Reverberi, A.P., Meshalkin, V.P., 1994, A general procedure for the estimation of diffusivities in solids, Canadian Journal of Chemical Engineering, 72, 1042-1046. Golub G.H. and Pereyra, V., 1973, The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate, SIAM J. Numerical Analysis, 10, 413–432. Kabanikhin, S. I., 2008, Definitions and examples of inverse and ill-posed problems, J. Inv. Ill-Posed Problems, 16, 317–357. Morozov, V.A., 1984, Methods for Solving Incorrectly Posed Problems, Springer, New York, US. Palazzi E., Currò F., Reverberi A., Fabiano B., 2014, Development of a theoretical framework for the evaluation of risk connected to accidental oxygen releases, Process Saf Environ, 92, 357-367. Reverberi, A.P., Fabiano, B., Dovì, V.G. 2013, Use of inverse modelling techniques for the estimation of heat transfer coefficients to fluids in cylindrical conduits, Int. Comm. Heat and Mass Transfer, 42, 1-94. Solisio, C., Reverberi, A.P., Del Borghi, A., Dovì, V.G., 2012, Inverse Estimation of Temperature Profiles in Landfills Using Heat Recovery Fluids Measurements, Journal of Applied Mathematics, ID 747410. Tikhonov, A., Arsenine, V., 1977, Solution of Ill-posed Problems, Winston & Sons, Washington, USA. Twomey, S., 1963, On the numerical solution of Fredholm integral equations of the first kind by inversion of the linear system produced by quadrature, J. Assoc. Comput. Mach., 10, 97-101. Yuan F. and Lu Z., 2005, Analytical Solutions for Vertical Flow in Unsaturated, Rooted Soils with Variable Surface Fluxes, Vadose Zone J., 4, 1210-1218. 1620