Approximation of slow moving interface phase change problems using a generalized Fourier series and the CVBEM

Approximation of slow moving interface phase change problems using a generalized Fourier series and the CVBEM T. V. HROMADKA II Director of Water Resources Engineering, Williamson and Schmid, 17782 Sky Park Boulevard, Irvine CA92714 USA Research Associate, Princeton University, New Jersey, USA

C. C. YEN Hydrologist, Williamson and Schmid, lrvine, California, USA Many important engineering problems fall into the category of being linear operators, with supporting boundary conditions. In this paper, a new inner-product and norm is developed which enables the numerical modeler to approximate such engineering problems by developing a generalized Fourier Series. The resulting approximation is the "best" approximation in that a least-squares (L 2) error is minimized simultaneously for fitting both the problem's boundary conditions and satisfying the linear operator relationship (the governing equations) over the problem's domain (both space and time). For slow moving interface phase change problems where the heat flux balance can be adequately described by the Laplace equation, the generalized Fourier series technique results in a highly accurate solution. PURPOSE O F PAPER In this paper, the mathematical development of the approximation procedure and the application of the technique to a heat transfer problem are presented. A detailed derivation of the technique and the application to several simple problems are given in Hromadka et alt,2 Because this technique shows considerable promise in many engineering applications, the greater computational effort involved over that needed with a finite element or finite difference method solution may be offset by the mathematical attractiveness of a convergence (in L 2 sense) which can be linearly programmed. INNER PRODUCTS FOR THE SOLUTION OF LINEAR OPERATOR EQUATIONS The general setting for solving a linear operator equation with boundary values by means of an inner product is as follows: Let f~ be a region in R" with boundary F and denote the closure of f~ by cl(t)). Consider the Hilbert space L 2(cl(t)), dlO, which has inner product (f, O) = ~fgdlt. (This is a real Hilbert space. For the complex version, use the complex conjugate of the function 0 in the integral.) The way to construct the necessary inner product for the development of a generalized Fourier Series is to choose the measure/z correctly; that is let/z be one measure/~1 on f~ and another measure/z 2 on F. One natural choice for a plane region would be for/zx to be the usual two dimensional Lebesgue measure dV on f~ and for/z 2 to be the usual arc length measure ds on F. Then an inner product is given by 3

(f,,g)=fnfgdV+~rfOds. Accepted May 1987. Discussion closes August 1988

© 1988

Computational Mechanics

Publications

(1)

Consider a boundary value problem consisting of an operator L defined on domain D(L) contained in L2(f~) and mapping into L2(Q), and a boundary condition operator B defined on a domain D(B) in LZ(t~) and mapping it into LZ(F). The domains of L and B have to be chosen so at least for f in D(L), Lf is in L2(fl), and for f in D(B), B f is in L2(F). For example we could have L f = V2f, and Bf(s) equal the almost everywhere (ae) radial limit of f at the point s on F, with appropriate domains. The next step is to construct an operator T mapping its domain D(7) = D(L)nD(B) into L2(cl(Q)) by (for example, Davis and Rabinowitz 4) Tf(x) = Lf(x) for x in f~

Tf(s) = BJ(s) for s on F. (2) From (2), there exists a single operator T on the Hilbert space L2(el(D)) which incorporates both the operator L and the boundary conditions B, and which is linear if both L and B are linear. An application of this procedure using the Complex Variable Boundary Element Method (CVBEM) is given in Hromadka et al s. In that study, L f = V2f and B f is the radial limit of f on F. Other applications are contained in Hromadka et al 1'2. Consider the inhomogeneous equation L f = at with the inhomogeneous boundary conditions B f = 02- Then define a function g on cl(Q) by g = gt o n g = g2 on F Then if the solution exists for the operator equation Tf=# the solution fsatisfies V2f = 01 on ~ and f = 02 on F in the usual sense of meaning that the radial limit of f i s 02 on F. One way to attempt to solve the equation T f = g is to look at a subspace D. of dimension n, which is

Engineering Analysis, June 1988, VoL 5, No. 2 95

Approximation of slow moving interface phase change: T. V. Hromadka H and C. C. Yen contained in D(T), and to try to minimize UTh -011 over all the h in D. such as developed in Hromadka, et al. 6 .

Definition of Inner-Product and Norm Given a linear operator relationship L(~b) = h on D., ~ --- 4~b on F

(3)

defined on the problem domain £~ with auxilliary conditions of ~ = ~bb on the boundary F(see Fig. 1). Here f~ may represent both time and space, and ~bbmay be both initial and boundary conditions. It is assumed that the working space is sufficiently restricted (see following) such that 4~ is a unique almost everywhere (ae) solution to (3). Choose a set of m linearly independent functions (fj).,, and let S" be the m-dimensional space spanned by the elements of ( f i ) m. Here, the elements of (f~)" will be assumed to be functions of the dependent variables appearing in (3). An inner-product is defined for elements of S" by (u, v) where for u, v e S"

(u,v)=rUOar+faLuLvm

(4)

It is seen that (u, v) is indeed an inner-product, because for elements u, v, w in STM (i) (u, v) = (v, u) (ii) (ku, v) = k(u, v), for L a linear operator (iii) (u + v, w) = (u, w) + (v, w) for L a linear operator

(iv)

u) = J'v (u)2dr + l'. (Lu) 2 m t> 0

(v) (u, u) = 0 =~ u = 0 ae on F, and Lu = 0 ae over fl The above restrictions on the operator L imply that L is linear (see (ii) and (iii) in the above definition); if Lu = 0 ae over f~ and u = 0 ae on F, this must imply that the solution u = [0], where [0] is the zero element over t)UF; and for the inner-product to exist, the integrals must exist. For the inner-product of (4) to exist, the integrands must be finite. Additionally, each element u 6 S~ must satisfy

~r u2dF < oo. For the above restrictions of L an the space S', the inner-product is defined and a n o r m " II II" immediately follows, II u II --- (u,

(5)

The generalized Fourier series approach can now be used to obtain the "best" approximation ~b,,eS" of the function ~b using the newly defined inner-product and corresponding norm presented in (4) and (5). The next step in developing the generalized Fourier series is to construct a new set of functions (g j ) " which are the orthonormal representation of the (f~)m.

Orthonormalization Process The functions (gj)., can be obtained by the wellknown Gramm-Schmidt procedure (Kantorovich and Krylov, pg. 45, 7) using the newly defined norm of (4). That is, Ox = f,/II A II

o., =

I f . , - (f.,, ol)o, . . . . .

IIf . , - ( f ' , 01)01 . . . . .

( f . , g.,- O g . - 1]/

(f,., 0.,-1)0,,-1 II

96 Engineering Analysis, June 1988, l/ol. 5, No. 2

(6)

1"

Fig. 1 Definition of Problem Domain, £L and Boundary, F. Hence, the elements of ( g y " properties that (gj, gk) =

satisfy the convenient

0, i f j ~ k 1, ifj = k

(7)

In a subsequent section, a simple one-dimensional problem illustrates the orthonormalization procedure of (6). The elements (g j)" also form a basis for S" but, because of (7), can be directly used in the development of a generalized Fourier series where the computed coefficients do not change as the dimension m of (ai>., increases. That is, as the number of orthonormalized elements increases in the approximation effort, the previously computed coefficients do not change. Each element ~b, e S" can now be written as 4'., = ~. ?jgj, 4', ~ S., j=l

(8)

where ~j are unique real constants.

Generalized Fourier Series The ultimate objective is to find the element ~b,,e S" such that II 4~..- 4~ II is a minimum. That is, we want II 4'., - 4' II 2 to be a minimum, where II 4 , ' - 4, II 2 = ~r

?jg j - Ldp) 2df~

(9)

Remembering that L is a linear operator, and L~b = f by the problem definition of (3), we have that (9) can be rewritten as j=l

Thus, minimizing II 4 ~ ' - 4' II 2 is equivalent to minimizing the error or approximating the boundary conditions and the error of approximating the governing operator relationship in a least-square (or L z) sense. Because the (g j ) " are orthonormalized and the innerproduct (,) is well-defined, the coefficients 7i of (8)

Approximation of slow moving interface phase change: T. V. Hromadka H and C. C. Yen are immediately deterimined by the generalized Fourier constants, ~', where r* = (gj, ~b),j = 1, 2. . . . , m

(11)

Thus ? * = ~ ?*gl= ~ (gj, ~b)gj j=I j=1

(12)

is the "best" approximation of q~, in the space S'. Because the generalized Fourier series approach is used, several advantages over a matrix solution (for the generalized Fourier series coefficients) are obtained: 1. Elimination of the need for solving large, fully populated, matrices such as occurs when solving the normal equations. 2. Elimination of the instability which typically arises in a matrix solution for Fourier coefficients (i.e., higher powers of the expansion basis functions assumed). 3. The generalized Fourier series coefficients do not change as additional functions are added (i.e., as the dimension m of the space S" is increased). 4. Generalized Fourier series theory applies; hence, error analysis can be conducted using Bessel's inequality as discussed in the next section.

Approximation Error Evaluation Due to the generalized Fourier series approach and the definition of the inner-product, Bessel's inequality applies. That is, for any dimension m in

(¢~,¢~)~ ~,, (gj, ~)2= ~ r~ 2 j=l

afforded by the CVBEM solution is the error analysis by use of the "approximate boundary" technique (see Hromadka9). In this paper, the CVBEM trial functions are used in order to eliminate the second integral in the innerproduct of (4). Thus,

(14)

j=l

where

(u, ,,) = fr

uvdF

(16)

becomes the inner-product for the generalized Fourier series development.

Modeling Approach The modeling approach (the governing equations and modeling assumptions are given in Hromadka, 19866) initiates by developing a CVBEM approximator 1° &I(Z) and &t(z) for the frozen and thawed domains, respectively. The numerical technique determines the analytic function t3(z) which satisfies the boundary conditions of either normal flux or temperature specified at nodal points located on the probem boundary, F. Because 6~(z) is analytic throughout the interior domain ~ which is enclosed by F, then the real and imaginary parts of ~(z) + i if(z) both exactly satisfy the Laplace equation over f~. For the steady state condition, the governing heat flow equations reduce to the Laplace equations. Consequently, an ~ z ) determined for both the frozen and thawed region satisfy the Laplace equations exactly, leaving only errors in satisfying the boundary conditions. To develop a CVBEM steady state solution, an &(z) is developed for each of the separate regions. Initially, both tSf(z) and &,(z) are defined by Co:(z) = C.}, z E

(~,qb)=fr(Ck)2dF+fa(L~)2dfl=;r ~b2dl-"+ f2dfl (15) Equation (15) is readily evaluated and forms an upper bound to the sum of(0j, ~b)2 as the dimension m increases. Consequently, one may interact with the approximation effort by carefully adding functions to the ( f j ) " in order to best reduce the difference computed by Bessel's inequality. The technique of reducing Bessel's inequality can be programmed by choosing additional basis functions which provide the greatest reduction in (14) from the set of basis functions available. APPLICATION TO SLOW MOVING INTER-FACE PHASE CHANGE PROBLEMS Many freezing/thawing phase change situations fall into the category of heat transfer problems where the heat flux along the phase change boundary is adequately estimated by assuming the Laplace equation. For example, in soil-water phase change in freezing soils, Hromadka and Guymon s successfully use the Laplace equation to compute heat flux quantities along a freezing soil column freezing front in order to propogate the front due to soil-water phase change. In another application Hromadka 6 uses the Complex Variable Boundary Element Method or CVBEM to extend the soil-water phase change solution to two-dimensions. A distinct advantage

69,(z) = 6~, z ~ t)

(17)

where t ) = ~ Uf~: is the global domain, and the first order CVBEM approximators are based on the entire domain. This procedure results in simply estimating 'the 0°C isotherm location for the homogeneous problem of f~ being entirely frozen or thawed. Let C 1 be the contour corresponding to this 0°C isotherm. The second iteration step begins by defining f ~ and D~z based on the mutual boundary of C 1. CVBEM approximators &~ and &~ are then defined for fl~ and f~z, respectively. Examining the stream functions ~} and ff~, estimates of the discrepancy in matching the flux rates along the interface between t2t and f~I can be evaluated. The ~ function is now used to determine the next location of the 0°C isotherm. This is accomplished by determining a new ~ with the stream function values of &~ (and modified by conductivity) superimposed at the nodal values of C 1. Next, a new 0°C isotherm C* is located for ¢5". The next estimated location for the 0°C isotherm, ~ , is located by averaging the y-coordinates of the nodal points between C 1 and C*. Figure 2 illustrates this procedure. The third iteration step proceeds by defining D~ and f ~ based on the mutual boundary of C z and the above procedure is repeated. The iteration process continues until the final estimates of f~: and f~t are determined with corresponding 59: and

Engineering Analysis, June 1988, Vol. 5, No. 2 97

Approximation of slow moving interface phase change: T. V. Hromadka H and C. C. Yen

oo -

6-

PfJ ~

"~

t

J

Of

u~oCONSTANT

4 --

~0)

~ • 1.0c,m.

"

Fig. 3 The Approximate Boundary f ' / and the Closenessof-fit to the Problem Boundary F s.

Fig. 2 Iterative Estimation of Freezing Front Location. &, approximators such that I k~d~s/ds - k,d~Jdsl < ~, s e C

(18)

where ks and k, are thermal conductivities for frozen and thawed regions, respectively. Using the Approximate Boundary As discussed previously, the subject problem reduces to finding a solution to the Laplace equation in f~s and O.t where f~/ and D.t coincide along the steady state freezing front location, C. The CVBEM develops approximators 6)s and &, which satisfy the Laplace equation over f ~ / a n d ~t, respectively. Consequently, the only numerical error occurs in matching the boundary conditions continuously on Fs, F,, and C. The generalized Fourier series develops the best CVBEM approximation which minimizes the norm in (15), where due to use of analytic functions as basis functions, fr (4), 4)) = 4)2 dF (19) To evaluate the precision in predicting the freezing front location, an approximate boundary is determined for each subproblem domain ofDs, ~ . The approximate boundary results from plotting the level curves of each CVBEM approximator (i.e., &s and &t) which correspond to the boundary conditions of the problem, For example, in flj, the thermal boundary conditions for a roadway embankment (Fig. 3) are defined on the problem boundary Fy by 4) = -10°C, z ~ top surface 4) = 0°C, z ~ freezing front ~b = 0, z e left side (symmetry)

Thus, the CVBEM modeling error is directly evaluated by the closeness-of-fit between 1~'/and F s. However, in this application, the approximate boundary concept is used not only to examine the closeness-of-fit to the boundary conditions, but possibly more crucial, the closeness-of-fit of matching the estimated freezing front location between fly and D~ along the contour, C. Should t~s and f~, match C continuously, then &j, and 69, equate thermal flux continuously along C. Applications Figure 3 depicts an application of the geothermal model for a roadway embankment problem and the use of the approximate boundary. Figure 4 illustrates the two-dimensional steady state freezing front location for a geothermal problem involving a buried subfreezing 3-meter diameter pipeline. An examination of the approximate boundaries indicate that a good CVBEM approximator was determined by use of a 26-node CVBEM model. The maximum departure 6 between the approximate boundaries and the problem boundary F occurred along the top of the pipeline and had a value of approximately 3.5 era. The average departure ~ is estimated at less than 1 cm. The freezing front maximum departure is approximately 4 cm and occurred at the right-hand side. Average departure on C is less than 2 cm. The example problems presented illustrate the usefulhess of the CVBEM in predicting the steady state freezing front location for two-dimensional problems. Possibly the most important result is the accurate determination of the approximation error involved in using the CVBEM. The usual procedure in estimating the freezing front is to use a finite element or finite difference numerical analog. A hybrid of these domain methods is to include a variable mesh in order to better accommodate the interface. However, none of these methods provide the error of approximation. In comparison, the CVBEM model provides the approximation error not only in matching the boundary conditions, but in predicting the interface location between t'ls and ~ . And this error is

~b = constant, z ~ right side (zero flux) After developing an tbj, and [l~r from the CVBEM, the approximate boundary 1~r is determined by plotting the prescribed level curves. The figure also includes l"y superimposed with r s. Because ¢bs is analytic within the interior of the approximate boundary and satisfies the prescribed boundary conditions on the boundary 1~/,then &I is the exact solution of the boundary value problem redefined on l='sand its interior, ~r. Should l~r completely cover F~,, then &s is the exact solution to the subject problem.

98 Engineering Analysis, June 1988, Vol. 5, No. 2

: _

~ ~

6-

? .I

,-

"[

e

~-

f

"

l~.-go '

f

~"AX~MU. O E P * ~

:

I

e~.5~

o, /-c~*-01 . a f, " • - :

.P.N'

~ / ~ Z " : -\' ~ ~

o - _-~-[ (MErERs) -- V_~ - f)t ~t'~. ~] M ( ~ - ~ ~ - ~ ~

I q,.co.~,~r

.J /k"*'-a° ~,~sr~cr I.V.. .=

,.u Fig. 4 Application of the CVBEM Geothermal Model to Predict Steady-State Conditions. 'f~"

Approximation of slow moving interface phase change: T. V. Hromadka H and C. C. Yen tional effort employed by the C V B E M analysis is adequate for this case study. The figure shows a variation in the approximate boundary location as the solution progresses in time; however, the variation is of less than 1.0 cm in magnitude. CONCLUSIONS In this paper, the C V B E M basis functions are utilized in a generalized Fourier series which uses the inner-product of (15) to determine the basis function coefficients. Because analytic functions are used, the inner-product reduces to a least-square fit of the boundary conditions. The C V B E M is used to approximate a slowly-moving interface between two quasi-potential problem solutions. The case study considered is soil-water phase change in freezing soils. The approximate boundary technique is used to demonstrate the C V B E M modeling error in achieving the prescribed boundary conditions as the time-stepped advancement in time is approximated.

11"51

ii•ol / / ~ M E

(HOURS| l

x

LEGEND •

NODAL POINT TRUE BOUNDARY APPROXIMATE BOUNDARY

Fig. 5 Approximate Boundary Evolution for TimeStepped Problem Solution (see _Fig. 4 for Domain Definition). simple to interpret as an a p p r o x i m a t e boundary displacement from the true problem boundary, and the displacement between ~ s and D~ along the freezing front contour, C.

Time-stepped Approximate Boundary By plotting the several C V B E M generated approximate boundaries, the time evolution of approximation error is readily seen. Figure 5 demonstrates the C V B E M modeling error in the time sequence of approximations developed for the pipe solution isolated from the Fig. 4 problem. F r o m Fig. 5, it is concluded that the computa-

REFERENCES 1 Hromadka II, T. V., Pindcr, G. F. and Jots, B. Approximating a linear operator equation using a generalized Fourier series: development, EngineeringAnalysis, 1987, 4, 2, 82 2 Hromadka II, T. V., Yen, C. C. and Pinder, G. F. Approximating a linear operator equation using a generalized Fourier series: applications, EngineeringAnalysis, 1987, 4, 4, 214 3 Birkhof,G., and Lynch,R. Numerical solution of elliptic problems, SIAM Studies in Appli~l Math, 1984 4 Davis, P. J. and Rabinowitz, P. Advances in orthonormalizing computation, Academic Press, 1961 5 Hromadka II, T. V., Lai, Chintu and Yen, C. C. A complex boundary element model of flow-fieldproblems without matrices, EngineeringAnalysis, 1987, 4, I, 25 6 Hromadka II, T. V. Predicting Two-dimousional steady-state freezing fronts using the CVBEM and an approximate boundary, Journalof Heat Tragsfer,ASME, 1986 7 Kantorovich, L. V. and Krylov, V. I. Approximate methods of higher analysis, Interscience Publishers, New York, 1964 8 Hromadka II, T. V. and Guymon, G. L. Simple model of ice segregation using an analytic function to model heat and soil-water flow, ASME Third Int. Symposium on Offshore Mechanics and Artic Engineering, New Orleans, LA, 1984a 9 Hromadka II, T. V. The complex variable boundary element method: development of approximative boundaries, Engineering Analysis, 19841>,1, 4, 218-222 10 Hromadka II, T. V. The complex variable boundary element method, Springer-Vedag, 1984c

Engineering Analysis, June 1988, VoL 5, No. 2 99