Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Generalised Stefan problems : linear analysis and computation Donaldson, Roger David 2003

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_2004-0103.pdf [ 4.36MB ]
JSON: 831-1.0080048.json
JSON-LD: 831-1.0080048-ld.json
RDF/XML (Pretty): 831-1.0080048-rdf.xml
RDF/JSON: 831-1.0080048-rdf.json
Turtle: 831-1.0080048-turtle.txt
N-Triples: 831-1.0080048-rdf-ntriples.txt
Original Record: 831-1.0080048-source.json
Full Text

Full Text

Generalised Stefan Problems Linear analysis and computation by Roger David Donaldson B . A . S c , The University of British Columbia, 2001 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in The Faculty of Graduate Studies Department of Mathematics Institute of Applied Mathematics We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A December 2003 © Roger David Donaldson, 2003 In presenting this thesis in partial fulfilment of the requirements for an ad-vanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permis-sion for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics Institute of Applied Mathematics The University Of British Columbia Vancouver, Canada Abstract i i Abstract We consider elliptic problems in which the domain is separated into two regions by a free boundary, on which mixed Dirichlet-Neumann conditions are speci-fied. Led by the classical Stefan condition applied to change of phase models, we consider idealised numerical methods which evolve the interface by a trial method that uses the error in one of the boundary conditions as the normal velocity of the boundary. Using linear perturbation analysis of simple cases, we show exactly which interfacial conditions lead to well-posed problems and which choices of velocities lead to convergent methods. Moreover, some veloci-ties lead to methods having superior numerical properties. Analysis of numerics representing the free boundary by a cubic spline fit is presented, followed by an example computation. Related ongoing work is introduced. Contents iii Contents Abstract ii Contents iii List of Tables v List of Figures vi Preface vii 1 Numerical Methods for Free Boundary Problems 1 1.1 Physical examples 3 1.2 Ice-water model 9 1.3 Solutions to free boundary problems 11 1.4 Convergence of explicit schemes 14 2 Linearised Problem 18 2.1 Base solutions 18 2.2 Criteria for well-posedness 20 2.3 Generalised Stefan velocities 24 2.4 Good velocity choices 26 3 Analysis of the Implementation 32 3.1 Fourier transform of the cubic spline 32 3.2 The cubic spline in generalised Stefan schemes 35 3.3 Smoothing Stefan velocity data 36 4 Computer Examples 39 5 Continuing Work 46 5.1 Generalised Stefan velocities for a level set method 46 5.2 Extensions to general operators and base boundaries 50 5.3 Boundary contact points 54 5.4 Summary 56 Bibliography 58 Contents iv A C a l c u l a t i o n De ta i l s 61 A . l Ellipticity of the porous media flow problem 61 A.2 Calculation of A 63 A.3 The Fourier integral of the cubic spline 64 A.4 The Fourier transform of periodic functions 65 A.5 Series solution to model ice-water problem 66 List of Tables v List of Tables 2.1 Well-posedness for feedback conditions 24 4.1 Maximum stable step sizes for each discretisation 43 List of Figures v i List of Figures 1.1 The dam seepage problem 4 1.2 The Alt-Caffarelli problem 5 1.3 Riabouchinsky flow 5 1.4 Flow of a Bingham fluid 7 1.5 Falling droplet 7 1.6 Two-phase flow in porous media 9 1.7 A model of the ice-water system 10 1.8 Forward Euler method stability 15 2.1 Normal to a curve with r]'(x) <C 1 21 2.2 Effect of velocity choice on the ice-water problem 30 2.3 Effect of velocity choice on the feedback problem 30 3.1 Sketch of sN(n) = s(n/N) 35 4.1 Evolution of the ice-water boundary by the Neumann method . . 41 4.2 Evolution of the ice-water boundary by the Dirichlet method . . 42 4.3 Mean square error of numerical solutions 43 4.4 Oscillation in the Neumann method 44 4.5 Accuracy of the Neumann and Dirichlet methods 45 5.1 Boundary contact point geometry 55 Preface vii Preface The original motivation for this work is a problem in fuel cell water manage-ment. The fuel cell system facilitates a controlled reaction between hydrogen and oxygen to create electric energy and water, and understanding the forma-tion and flow of this water in the fuel cell system plays a role in optimising fuel cell performance. This motivating problem seeks to determine the boundaries between liquid water, water vapour, and mixed liquid-vapour regimes within the porous media that permits water transport within the fuel cell. This problem is a free boundary problem, where not only state variables, but also phase bound-aries, are unknown a priori. While the work of this thesis does not address the fuel cell problem directly, it analyses a novel strategy for handling the general class of free boundary problems of which the fuel cell problem is a member. I wish to acknowledge my supervisor Brian Wetton for not only bringing this problem to me as a masters thesis project, but also for his ideas and support in its development. Work presented here is joint with Dr Wetton. I also wish to acknowledge N S E R C , Ballard Power Systems Inc, and The University of British Columbia Department of Mathematics for providing financial support for this work. I look forward to seeing the application of ideas contained herein to fuel cell development in the near future. Roger Donaldson Vancouver, Canada August 2003 Chapter 1. Numerical Methods for Free Boundary Problems 1 Chapter 1 Numerical Methods for Free Boundary Problems Partial differential equations (PDEs) often describe physical laws, such as con-servation of mass, momentum, or energy. The solution of a P D E system gives information about the physical situation by describing how state variables, such as velocity or temperature, vary throughout a domain. In many cases, the boundaries of the domain are known, and boundary conditions, together with the P D E system, provide sufficient information to solve for the state variables. However, some physical situations do not specify the location of the domain boundaries, instead providing additional boundary information such that the state variables and the boundary location may be determined. Such problems are referred to as free boundary problems. A typical free boundary problem statement includes a fixed physical bound-ary, but within the fixed domain, boundaries with unknown locations separate subdomains where different equations may apply to different sets of dependent variables. Even if the P D E s that apply in the subdomains are linear, free boundary problems are, on the whole, nonlinear, and are often difficult to solve numerically, let alone analytically [19, ch 7]. Free boundary problems can be time-dependent, and study of the time-dependent case has led to interest in so-called front tracking problems. However, efficiently determining steady-state solutions of free boundary problems is alone challenging, and the purpose of this thesis is to discuss a novel numerical method for solving steady-state free boundary problems, motivated by the classical Stefan velocity for phase change models. While a variety of techniques exist for numerically solving free boundary problems, we look in particular at trial methods, also called fixed point methods. At a free boundary, in addition to the number of boundary conditions required to specify the problem were the boundary fixed, one further boundary condition must be provided to specify the free boundary position. Trial methods solve for the free boundary of a steady-state system as follows: 1. Choose an initial approximation to the free boundary. 2. Solve the resulting fixed boundary problem using all but one of the free boundary conditions. 3. The remaining condition will not be satisfied unless the position of the free boundary is accurate: use the residual in the remaining condition to Chapter 1. Numerical Methods for Free Boundary Problems 2 adjust the location of the boundary. 4. If the adjustment required is sufficiently small, stop; otherwise, adopt the new boundary and return to step 2. Our focus is on the manner in which step 3 is carried out. We treat the error in the remaining condition as a normal boundary velocity, and use this veloc-ity to track the boundary toward its steady state solution. In addition, we recognise that the choices of which boundary conditions are applied to solve the fixed boundary problem, and which is used as the velocity, are arbitrary. Some choices give the trial method superior numerical properties, while others give unstable methods. While we have already stated that free boundary prob-lems are inherently nonlinear, addressing questions of uniqueness is beyond the scope of this thesis, and we instead confine our discussion to issues surrounding stability of the trial method using various velocity choices. The work of this thesis uses tools of applied mathematics and computation accessible to senior undergraduate students of mathematics, physics or electrical engineering. Some of the signal processing tools may be unfamiliar to mathe-matics students, while some of the linear algebra and asymptotic analysis may be unfamiliar to students of physics or engineering. While this thesis is intended to be self-contained, the tools used, together with some suggested references, are as follows: The overall structure of the stability analysis is adapted from the theory of numerical solutions to ordinary differential equations (ODEs): for a discussion of numerical methods for, ODEs, cf. Dormand's text [10]. See Okendon et. al. [19, ch 7] for an introduction to free and moving boundary problems, including the form of our basic stability analysis. Motivation for the perturbation method used to evaluate stability can be found in Gustafson et. al. [15], and various smoothing and Fourier transform techniques can be found in signal analysis texts, cf. Hamming's book [16]. The physics of heat transfer may be studied from any number of engineering texts on the subject, but the development of the heat equation in any introductory P D E s book, such as that by Boyce and DiPr ima [5] should suffice as a reference. Linear algebra is used extensively in our arguments: the text by Strang [28] is an excellent resource. For a comprehensive text on the finite element method, balancing theory with practical discussion, see Oden and Reddy [20]. Finally, Crank's text [8], con-tains an analytic treatment of classic free boundary problems, with its closing chapters discussing numerical methods. We now proceed to several physical examples of free boundary problems. The first examples motivate the solution of free boundary problems in gen-eral, while the last of these examples provides motivation for our numerical method. This introductory chapter closes with a discussion of methods for solv-ing free boundary problems, with focus on the convergence of our trial method. Chapter 2 evaluates the stability of the trial method applied to a linearised version of the free boundary problem, showing which free boundary problems are well-posed, and which choices of velocity conditions give superior numerical methods. Chapter 3 shows that the discretisation and smoothing operations in practical schemes preserve the numerical properties predicted by the idealised Chapter 1. Numerical Methods for Free Boundary Problems 3 linear analysis, and Chapter 4 presents the performance of the practical scheme in example computations. Chapter 5 discusses continuing work, introducing a hybrid scheme combining the benefits of a level set method with good numerical qualities of our trial method, and presents questions regarding points of contact between free boundaries and the physical boundary of the problem. Finally, an appendix provides details of some of the longer calculations. 1.1 Physical examples Our analysis deals only with linear, constant coefficient boundary conditions. We also restrict our attention to the case of problems in two spatial dimensions where the elliptic operators on either side of the interface are the vector Lapla-cian, with some mild extensions possible as outlined in Chapter 5. However, we believe our analysis and the insights into robust and fast numerical methods that can be gained from it are extensible to much more general situations. We outline below some of the examples of the type of problems in our general class that motivated our work and to which our methods could be applied. The examples are roughly in order of increasing complexity. Having a single dependent variable, the first four examples are the lowest order of the class of free boundary problems we consider. They consist of a scalar elliptic problem in a domain with a free boundary at which two mixed Dirichlet-Neumann conditions are applied. The next two examples have two state variables on either side of the free boundary, a total of four dependent variables, necessitating five free boundary conditions. Despite the range of complexity, all of these problems belong to the single class of free boundary problems we may treat using our trial method. Dam seepage A simple free boundary value problem is outlined in some detail in [8, ch 2]. Several techniques are also outlined for the solution of this problem in later chapters. It is the problem of finding the seepage region Cl (the region in which liquid water will be present) in a porous dam separating two reservoirs of water at different levels. A diagram of the situation is given in Figure 1.1. The reservoir levels j / i and 2/2 and the dam width L are given. The free boundary V (y — f(x)) must be determined. The equations that govern the system are Acj) <P dy 2/i on x — 0 2/2 on x = L 0 i n l l (1.1a) (1.16) (1.1c) 0 on y = 0 ( l . l d ) y on 7 y on r ( L i e ) (1-1/) Chapter 1. Numerical Methods for Free Boundary Problems 4 Figure 1.1: The dam seepage problem. The reservoir levels y\ and j/2 and the dam width L are given. The free boundary r (y = f(x)) must be determined. ^ = 0 o n T (l.lg) where <j> is the modified velocity potential. Equation (1.1a) is the elliptic problem for the scalar quantity (j>. One boundary condition (1.16-l.le) is given on each physical boundary. Two boundary conditions (1.If,l.lg) (one Dirichlet and one Neumann) are given on the free boundary T. Thus, this is an example of the kind of general problem we consider in this work. We note that additional constraints on the problem need to be specified on 7 and where the free boundary T meets the physical boundary at x — 0 and x = L to obtain a physically reasonable solution. We refer the reader to [8] for details. Alt-Caffarelli problem Another scalar example is the Alt-Caffarelli problem [3], also known as the exterior Bernoulli problem [11]. Several applications in which this problem arises (fluid flow, optimal insulation, electrolytic galvanisation) are discussed in [11]. The schematic of this problem is given in Figure 1.2. Here, the equations that apply are 0 in 1) (1.2a) 1 o n E (1.26) 0 o n T (1.2c) Q on T (1.2d) where Q < 0 is given. Here, (1.2a) is the elliptic problem for the scalar u, (1.26) is the single boundary condition on the physical boundary E and (1.2c,1.2d) are the two conditions (one Dirichlet and one Neumann) on the free bound-ary T. The conditions on the free boundary are very similar to those for the Au = u = u = du dn Chapter 1. Numerical Methods for Free Boundary Problems 5 r Figure 1.2: The Alt-Caffarelli problem. The physical boundary E is given. The free boundary T must be determined. Figure 1.3: The Riabouchinsky flow past circular disks. Here, the disk size R and spacing L are given as well as the strength of the external flow past the disks. The free boundary T is to be determined. dam problem above (1.1/, 1.1s)- In f a c t the only difference at T is that the Alt-Caffarelli problem is forced through the Neumann condition and the dam problem through the Dirichlet condition. The physical boundary conditions are also different. Riabouchinsky flow past circular disks We discuss here a problem presented in [14] on the inviscid, irrotational flow past two disks, separated by a distance L, positioned perpendicular to the x-axis with their centres on this axis. A free surface Y is to be found that connects the perimeters of the disks. This surface separates the external fluid flow from a region of cavitation. The situation is shown pictorially in Figure 1.3. The flow is axisymmetric and can be described by a stream function ip(x,r). r / x r Chapter 1. Numerical Methods for Free Boundary Problems 6 The equations that apply here are 0 in fi (1.3a) 0 on the disks (1.36) 0 on T (1.3c) r on T (1.3d) where here fi is the exterior flow domain. As before, we have one scalar el-liptic equation in the domain (1.3a), one boundary condition at each physical boundary (1.36), and two conditions at the free boundary (1.3c, 1.3d). Far field conditions must also be given to solve the problem. They can be considered physical boundary conditions at "co". r > = dn Hagen-Poiseuille flows of a Bingham fluid Bingham fluids are particular non-Newtonian fluids that have a yield stress, a minimum stress below which they behave like solids, cf. [13]. We consider here the steady, pressure-driven flow of a Bingham fluid through a pipe of constant cross-sectional shape E , taken to be in the xy plane. We consider flows in which fluid moves only in the direction along the pipe with speed w(x, y). It is possible that an unyielded plug forms away from the pipe walls as shown in Figure 1.4. This plug moves as a solid with common speed U, to be determined along with its shape I \ It is also possible that unyielded regions form stuck to the walls of the pipe depending on the shape of E . These unyielded regions will have zero velocity. This is just an indication that for most free-boundary value problems from applications, the topology of the free boundary is not know a priori as in the preceding examples. For just a single unyielded region away from the pipe boundary, the following describe the problem: - A in fi (1.4a) 0 on E (1.46) U onV (1.4c) 0 on T (1.4(f) where B is the yield stress of the fluid (a given physical parameter), —A is the given, imposed pressure gradient and U is the speed of the plug. Although nonlinear, the problem has the same structure as the previous ones: (1.4a) is a scalar elliptic equation, (1.46) is a single condition on the physical boundary, and (1.4c, lAd) are conditions on the free surface. The velocity U is a global unknown to be determined as part of the problem. V B \Vw\ w dw dn Chapter 1. Numerical Methods for Free Boundary Problems 7 Figure 1.4: Pressure-driven Hagen-Poiseuille flow through a pipe of constant cross-sectional shape E . A n unyielded region inside T can develop. Here, the material moves as a solid with constant speed U. Both the shape of T and U must be determined as part of the problem. U air Figure 1.5: The falling droplet in a moving coordinate system. Falling droplet We consider the problem of a steady state droplet falling with speed U. Wi th the speed given, the mass of the droplet falling with that steady state speed will be determined. We present the situation for a two-dimensional droplet to avoid the complexity of axisymmetric coordinates but the system has similar character in that more realistic case. We consider the droplet in a moving frame of speed (0, —U) so it will appear stationary as shown in Figure 1.5. The equations that apply here can be derived from equations for unsteady flow, cf [7]. pu • V u + Vp — / J A U = 0 in both air and water (1.5<z) V • u = 0 in both air and water (1.56) Chapter 1. Numerical Methods for Free Boundary Problems 8 pn + p du dn [u] u • n ->• (0,U) far from droplet (1.5c) = Kofi on T (1.5d) = 0 o n T (1.5e) = 0 o n T (1.5/) where p is the density and p the viscosity, both different for air and water; u is the vector velocity; p is the pressure; [•] denotes the difference in quantities from the water to the air side of the free surface; n is the unit normal direction of the free surface; K is the curvature of the interface; and a is the surface tension. We can show that this system falls into our class of free boundary problems. First, it is appropriate to consider (1.5a, 1.56) as a rank 2 system of elliptic equations in both air and water; it is possible to rephrase the problem in the right function spaces eliminating the pressure, see [30]. Equation (1.5 c) represents two physical (far-field) boundary conditions for the air velocity. If we consider the surface tension as a forcing term, the stress-continuity equations (1.5d) are two mixed Dirichlet-Neumann conditions on T. Continuity of velocity (1.5e) gives two more Dirichlet conditions and the steady-state condition (1.5/) gives another. Thus, we have 5 total mixed Dirichlet-Neumann conditions for this system of four dependent variables, making this another example of our general class of problems. T w o - p h a s e f l o w i n p o r o u s m e d i a We consider below a problem of two-phase flow in porous media, a two-dimen-sional version of the problem studied mathematically in [6], and by experiment in [34]. The problem considers a closed system containing a sand pack, water vapour and liquid water. If the system is heated at the top, it is possible to form two distinct zones separated by a free boundary. In the lower, two-phase zone, both vapour and liquid water can be found. In the upper zone, there is only water vapour. Wi th low temperatures or enough total water in the system it is possible to have another free boundary to a third zone of pores filled with water, but we do not consider that case here. We consider the steady, two-dimensional case shown in Figure 1.6, taking the solution periodic in x for simplicity. The variables pressure Pv(x,t) and temperature T(x,t) describe the flow in the upper region, and are governed by A T = 0 (1.6a) V-(Cu) = 0, (1.66) where C and u are functions of T and Pv. Flow in the lower region is described by temperature and saturation s(x,y). s is the fraction of pore space occupied by liquid water. The equations are A T = -KTV{piui) V(piui+pvuv) = 0. (1.6c) (1.6d) Chapter 1. Numerical Methods for Free Boundary Problems 9 T2(x) given vapour only zone 2-phase zone \; y Tl(x) given Figure 1.6: Two-phase flow in porous media. The upper temperature Ti(x), the lower temperature T2{x) and the total water content W in the domain are given. The two-phase zone boundary Y is to be deter-mined. where uj, uv, and pv are functions of T and s. KT is a constant. The rela-tionships between the intermediate quantities and state variables give elliptic problems for T,P and s. Details are given in Appendix A . l . On the free boundary, we have the following conditions: s = 0 (1.7a) [T] = 0 (1.76) [Pv] = 0 (1.7c) (pvuv)+-n = (pvuv + piui)~ • n (1.7d) 1^ = Kripm)-.. (1.7e) These are five Dirichlet-Neumann conditions required for our problem of four variables. Wi th additional conditions on the fixed boundaries, as well as a global condition, the problem is well-posed. This is too complex an example to develop our numerical scheme, but this problem demonstrates the range of problems to which our scheme applies. This problem is not addressed further in this thesis, but will be solved by the methods contained herein in future work. 1.2 Ice-water model This final example is a two-dimensional model of temperature in an ice-water system with a phase interface, Y. This example is a two-sided problem, but with Chapter 1. Numerical Methods for Free Boundary Problems 10u~(x, - 1 ) < 0 Figure 1.7: A model of the ice-water system. Two Dirichlet and one Neumann condition must be satisfied at T at steady state, although these conditions can be generalised to problems with other free boundary conditions. a scalar elliptic problem on each side of the boundary for a total of two dependent variables, hence three conditions at the free boundary. From this last example, we develop our framework for making good choices for the velocity condition in our trial method. Motivated by the physics, we take our first step in developing our numerical scheme for solving free boundary problems. Figure 1.7 illustrates the problem of a box filled with water, with temperature, u, held at u~(x, —1) < 0 on the bottom boundary and u+(x,l) > 0 on the top boundary. A phase boundary, where u = 0, appears as y = n(x), dividing the domain into an upper, water-filled subdomain, D+, and a lower, ice-filled subdomain, D~. A t steady state, the temperature obeys A u * = 0, where u+ and u~ represent the temperatures above and below T, respectively. Thus, at T, we have the two conditions u+ = u~~ = 0. Since T is a free boundary, an additional condition is again required to completely specify the problem, and for this third condition, we look to the time-dependent problem. In the time-dependent problem, the phase boundary moves locally as ice melts to water, or as water freezes. The normal velocity of the boundary, vn, is given by the classical Stefan condition on heat flux normal to the boundary, a is the ratio of thermal conductivities of ice to water and K is a physical con-stant of proportionality 1. At steady state, vn = 0, and we use this as the third 1In supercooled liquid interface problems, a small surface tension term is added to the RHS of (1.8). This term, while negligible in our situation, is needed in the supercooled case to obtain physically reasonable (well-posed) interface motion [24]. Chapter 1. Numerical Methods for Free Boundary Problems 11 condition on T. This representation suggests a particular trial method for it-eratively solving the free boundary problem by explicit time-stepping: Assume an approximation for T, solve the problem for u+ and u~ in the resulting sub-domains using the conditions u+ = u~ — 0 at V, and then move T normally proportional to the Stefan velocity given by equation (1.8). When we discre-tise T, as is necessary in a numerical solve, this method treats the behaviour of each of the points of T as the solution to a time-dependent O D E . The iter-ation discretises "time" as the boundary evolves toward a steady state, yet as we are solving a steady-state free boundary problem, we give time no physical interpretation, but refer to it only in the sense of the progression of the solve. This method converges, but, as we shall show, the stability of this method is dependent on the discretisation we choose for T; finer discretisations for T re-quire that shorter time steps be taken to maintain stability. The fixed-boundary problem must be solved each time T is adjusted, and, as this solve can be ex-pensive, it is desirable to take time steps as large as possible, independent of the discretisation of the boundary. Fortunately, there is another option. We can choose to set the Stefan velocity equal to zero for each fixed-boundary solve in the iteration. We enforce one other boundary condition, say, u+ = u~, and use the error in the remaining boundary condition, the value of u+ = u~, which we expect will ultimately tend toward zero, to determine the normal velocity of T. This is the essence of our approach: If a good choice is made for the velocity condition, the stability of the method can be made independent of the boundary discretisation. Note that out of deference to the classical Stefan velocity, we have adopted the term "velocity" to refer to the error in the boundary condition used to adjust the boundary at each iteration. We return to the ice-water system to compare the numerical performance of the trial method using the physical Stefan velocity with the that using the alternative "Dirichlet" velocity in Chapter 4. We will use this system as moti-vation for considering the stability of problems with more general conditions at F. Our analysis will be restricted to the case where A u = 0 is the problem in each subdomain. However, we expect that our stability predictions will hold for other elliptic operators, an important consideration for the other applications suggested. Before presenting a study of our numerical scheme, let us look at methods available for solving free boundary problems, with particular attention to their convergence. 1.3 Solutions to free boundary problems We are not aware of any analytic theory covering our general class of free bound-ary problems, although there'has been progress on individual problems. Refer-ences to early work can be found in [8], including a lengthy discussion on the dam problem. Existence and regularity of the Alt-Caffarelli problem is shown Chapter 1. Numerical Methods for Free Boundary Problems 12 in [3], extensively using a reformulation of the problem as a functional minimi-sation. A n exact solution for the ice-water example of the previous section does exist: the interface is the zero level set of a harmonic function after a reformula-tion shown in Appendix A.5 . In [6] semi-analytic solutions and their regularity for the two-phase zone problem are presented, but only for a one-dimensional problem. Given the difficulty in formulating analytic solutions for free boundary prob-lems, seeking good numerical schemes for solving these problems is desirable. Several techniques are discussed in the literature, of which we introduce cap-turing techniques and front tracking. Front tracking encompasses both shape optimisation and our trial method, the latter of which we expand upon with respect to its stability. Capturing techniques A useful method for solving the type of problems we have considered is to reformulate the problem so that modified equations are solved in the whole domain (both sides of the interface). The interface is then recovered from the resulting solution. The Extended-Lagrangian formulation of Bingham fluids [13] and the entropy formulation for the Stefan problem [8] are examples of this approach. Not all problems can be solved using this idea, although a variety of related approaches are available. The falling droplet problem can be approached using the methods developed in [7]. In their method, they discretise both the air and water as a uniform fluid with a transition zone in density and viscosity at the interface. This transition zone automatically approximates the interface conditions. The level set method of Osher and Sethian [23] is used to capture the interface location and provide the curvature for the surface tension term. The level set method is able to capture the topological changes of merging droplets and bubbles. In total, this is a delicate and sophisticated method. We return to the level set method in Chapter 5, where we combine it with our trial method to capture the good qualities of each approach. It may be possible to reformulate the two phase flow problem and cap-ture the interface by relaxing the saturation assumption in the two-phase zone and allowing condensation to happen at a finite rate proportional to the over-saturation [12]. However, this is not a straight-forward process either, and further study would be required to extend this idea to more general problems. Front tracking and shape optimisation Another numerical approach is to discretise the free boundary directly, forming a problem where the mismatch of an approximate boundary to the desired one is minimised. Such techniques are called front tracking, for their analogues in time-dependent free boundary problems. This is exactly the form of our trial method, whereby each point on the boundary discretisation is moved normally according to a chosen velocity. Chapter 1. Numerical Methods for Free Boundary Problems 13 To our knowledge, work exploring the effects the choice of velocity condi-tion has on a trial method reduces to work by Garabedian [14], in which he addresses the one-sided Riabouchinsky flow problem (1.3a-1.3d). Garabedian's idea is instead of using the first or second boundary condition alone in the fixed boundary P D E solve, to choose a combination of the two conditions that does not change along the normal to T, a stationary combination. Garabedian uses the error in the Dirichlet condition u = 0 to adjust the boundary. For this particular problem, just as in Newton's method, Garabedian's method gives a quadratic rate of convergence. If u = 0 is used to solve the fixed boundary problem, with error in the Neumann condition used to adjust the boundary, the convergence rate is linear. In their review paper on free boundaries of liquid in capillary tubes, Cuvelier and Schulkes [9] reproduce Garabedian's idea from a similar perspective. They determine how to adjust the free boundary by linearising solutions of the free boundary problem near its supposed location. The linearised version of the problem local to the boundary gives Garabedian's method, which, owing to its derivation, they refer to as total linearisation. They apply total linearisation to a one-sided free boundary problem similar to that considered by Garabedian. Several other authors have since reproduced Garabedian's idea using free boundary adjustment by shape optimisation [11, 18, 31-33]. Shape optimisation solves free boundary problems by a variational approach, minimising an integral functional that integrates a cost function over the variable domain having a free boundary. Shape optimisation uses shape calculus, which takes derivatives of an integral functional with respect to changes in its domain of integration. Ultimately, then, shape optimisation uses Newton's method to determine the boundary location where the shape derivative is zero. For an introduction to the shape calculus, cf [26, 27]. Garabedian's idea manifests itself in the context of shape optimisation when the freedom in the variational formulation is used to keep one combination of free boundary conditions stationary as the boundary is adjusted. This is ex-actly Garabedian's approach, and results in the error in u being, from the vari-ational perspective, the correct condition to use to adjust the free boundary. Again, these authors consider the one-sided free boundary problems addressed by Garabedian, and report the superior convergence properties of the Garabe-dian trial method. To avoid the calculation of shape derivatives, which is really only feasible for fairly simple systems, other techniques using only the values of the interface residuals have been developed, cf. [8]. In certain situations, such as in [14], they can be made to work almost as well as shape optimisation (Newton) methods. In some situations they are able to move a poor initial estimate to the desired interface. In this work, we consider a general framework to analyse value-type methods and identify those with superior numerical properties. Our method looks generally at choices for the velocity condition, and while rate of convergence, the primary motivation of the above authors, is important, we look instead at the stability of methods using different velocity conditions. We also generalise the problem to systems having many dependent variables Chapter 1. Numerical Methods for Free Boundary Problems 14 on either side of the free boundary, a feature particularly important in our application to phase change models. We treat the free boundary iteration as a semi-discretisation, tracking the evolution of the boundary in time. This allows us to evaluate the convergence of free boundary solves in the framework of numerical methods for ordinary differential equations (ODEs). 1.4 Convergence of expl ic i t schemes The trial method is an explicit scheme: the future location of T in the iteration depends only on its current (and, in a more general setting, previous) location. Such schemes are convenient to compute, and often converge rapidly, but are prone to instability if not applied carefully. That explicit iterative solutions to free boundary problems are prone to instability is known, cf. [25], and we must ask how stability depends on features of the problem. This section formalises the compromise between stability and convergence rate in terms of the problem stiffness. We begin with a simple example. Consider the O D E <"> where y(0) is given. A is the growth rate of the problem, and constraints on how we solve this problem depend on the size of A. The lowest order explicit scheme we can use for this problem'is the forward Euler method, for which we have the approximation yn+1 =yn + Xkyn, (1.10) where k is the time step size, and if y° = y(0), yn approximates y(nk). If Re(A) < 0, the exact solution decays, so, too, we should expect the ap-proximation to decay. Hence, we apply the criterion that the growth rate, \yn+i/yn\ < 1, so that | l + AJfc|<l. (1.11) If this condition is violated, the numerical solution will grow, an unacceptable representation of the decaying exact solution. We use condition (1.11) as the stability condition for the forward Euler method applied to this problem, and we note that our choice for k depends on the intrinsic growth (decay) rate of the problem at hand. Figure 1.8 is a sketch of the stability regime for the forward Euler method as a function of Xk. If A is purely real, we require — 2 < Xk < 0, and if we would like yn to decay to zero as quickly as possible, we choose k = -j (1.12) In which case y1 = 0. This is different than considering how to choose k such that our numerical solution best approximates the exact solution 2. In our trial 2 H i g h e r order R u n g e - K u t t a schemes for approximating (1.9), those that obtain better accuracy for a given k, do exist, but, in determining steady-state solutions, because we are not concerned with the accuracy of the decay dynamics , higher order schemes provide us no advantage over the forward Euler method. Chapter 1. Numerical Methods for Free Boundary Problems 15 Re(Xk) Figure 1.8: Stability regime of the forward Euler method. The method is stable as long as Xk lies inside the circle. Note the allowance for both real and complex A. method, y represents the error in the free boundary location, and we wish to see this error reduced to zero in as few iterations as possible. This simple development applies to solving for the scalar y = y(t). Let us extend this idea to our trial method. Consider our numerical method. Let y be a vector representing the y-locations of N points discretising the free boundary. At each iteration, we adjust each j/-value of the boundary to a new location by y»+! = y" 4- kvn. (1.13) Again, vn is the normal velocity of the boundary, as determined by the velocity condition. (Strictly speaking, we should be moving each boundary point nor-mally to r, but this is a good approximation if T is fairly flat. We will make this assumption rigorous in the next chapter.) vn = An(y), A™ : RN ->• RN, is a function that computes the normal velocities of the boundary points based on the solution to the associated fixed boundary problem at iteration n . A n(y), therefore, encompasses the P D E solve at iteration n and depends on the choice of velocity condition. Let us assume that A" is linear. This implies that the solution to the free boundary problem is y = 0. This assumption is carried through the analysis of succeeding chapters and, despite the fact that free boundary problems are Chapter 1. Numerical Methods for Free Boundary Problems 16 inherently nonlinear, choosing linear A n does serve as a good model for the behaviour of our iterative scheme. We omit the n superscript from here on, with the implicit understanding that A may change from iteration to iteration. In a practical problem, the solution will not likely be y = 0, but if we instead regard y as the error, the difference between our approximation and the exact solution, the assumption that A is linear for small error is reasonable and using this model makes sense. Suppose that A is negative definite, so that the velocity decays as y ap-proaches the solution, y = 0. Furthermore, suppose that A has eigenvalues Ajvf = Ai < A2 < ... < AAT = A m < 0. Thus, equation (1.13) becomes y»+! = y« + kAy? = yf(l + k\i) (1.14) for each eigenvalue A; and each eigenvector yf. We would like to choose k that keeps the solve stable while advancing it as quickly as possible. That is, we wish to know the value of k that minimises the maximum growth rate r of the problem, dependent on {Aj}. Hence, r = m i n { max {|l + fcA|}l. (1.15) fc>o 1 A G { A , } 1 J To compute r and the k at which r is attained, we observe that all values | l + fcA| lie between |1 + k\m\ and |1 + fcA^I- Since A M < Xm < 0, the minimising k occurs when - ( l + fcAm) = l + kXM, k = , ~ 2 , , (1.16) and so r = (1.17) K + l where K = A M / A t o is the condition number of A . K is also referred to as the problem stiffness, large values of K representing stiff problems. Thus, equa-tion (1.17) shows that when K 3> l , r « 1 and y™ + 1 RS y " . Since linearity of A implies that y = 0 is the solution, such stiff problems converge slowly. However, when K RJ l , r C 1 and for any y " , y™ + 1 w 0, indicating a rapidly converging problem. Having the freedom to choose the velocity condition for the trial method is synonymous with having the freedom to choose stiffness qualities of A . Equa-tion (1.16) tells us how to best choose k for a given A , but more importantly, expression (1.17) for r tells us what the eigenvalues of optimal A should be. The closer K is to l,.the faster the convergence to y = 0 will be for optimal k. The central idea of this work is to capitalise on the freedom in our choice of velocity condition to minimise the stiffness of the problem. In the next chapter, we calculate the eigenvalues of the problem with general free boundary conditions where the calculation is done for a linearised problem Chapter 1. Numerical Methods for Free Boundary Problems 17 over an infinite domain. In this setting, the problem has an infinite number of eigenvalues, but we can immediately anticipate how the distribution of eigen-values affects the stiffness of related discrete problems. Chapter 3 applies the results from the linear analysis to a practical scheme with a discretised free boundary, hence a discrete set of eigenvalues, which more closely follows the motivation from this introductory discussion. Chapter 2. Linearised Problem 18 Chapter 2 Linearised Problem Our trial method for tracking boundary points is an ODE-type semi-discre-tisation. We use a linearised version of the problem to determine conditions for stable numerical schemes that use explicit semi-discretised front tracking methods to evolve the free boundary. The stability analysis takes a known base solution, perturbs the free boundary in that solution by a single wave-mode, and checks for growth or decay of that mode. The wave-modes carried by T depend on its discretisation, but for this initial analysis, we compute how the stiffness, hence maximum step size, depends on wavenumber, and treat the discretised boundary in Chapter 3. This procedure is equivalent to a spectral analysis of the linearised problem. We shall see several major results in this linear analysis. Before proceeding to the stability analysis, we examine the existence of base solutions for general linear boundary conditions. This investigation will indicate that the behaviour of the free boundary can depend not only on the boundary conditions them-selves, but also on global conditions of the problem. Using these base solutions, we use the form of the perturbation analysis to show which boundary and global conditions give well-posed problems: conditions at the free boundary cannot be specified arbitrarily. Finally, after developing the relevant formalism, we shall come to the point of our analysis having direct impact on the numerics: we pre-dict the effect that velocity choice has on the stability of the explicit iterative scheme for determining the free boundary location. 2.1 Base solutions We consider the problem in E 2 . While the ice-water problem has only two dependent variables, u+ and u~, we are by no means bound by this restriction, and instead of 2 elliptic problems, we consider m elliptic problems for each of m dependent variables. Each of the variables may apply variously above and below the free boundary. We assume that A u = 0 in both D+ and D~ (the upper and lower subdomains, respectively), and that the free boundary, T, can be represented as a function, y = i](x). We consider the domain to be infinite in y and in x. We assume that the m+1 conditions at V are linear in u and u „ , that is, Gii = b (2.1) where G is an m+l-by-2m matrix. Hence, G accounts for the additional con-dition needed to allow for the freedom in the boundary location. We assume Chapter 2. Linearised Problem 19 no redundancy in the boundary conditions, so G is full rank, b is a constant vector representing any driving conditions at T, and u = [ u n , u ] , the normal derivatives and values of u at the boundary. Boundary normals are taken in the positive y-direction. In the case of the ice-water system, for example, G = 1 - a 0 0 0 0 1 0 0 0 0 1 (2.2) b = 0, and u = [u+,u~,u+,u~]T. Equations (2.1) may be row-reduced without losing boundary data, so this form for linear boundary conditions makes it clear why conditions for one-sided free boundary problems, those for which m = 1, ultimately reduce to those studied by Garabedian: with only two boundary conditions to specify and only one dependent variable, the only choices are u = const and un = const. We assume asymptotic expansions for the boundary y = n(x) and the asso-ciated solution u : r\ ~ r]o+en1 +... (2.3) u ~ u 0 + e u i + . . . . (2.4) Our goal is to determine base solutions, u n , for the case where the base boundary is y = r]o{x) = 0. Our differential operator is the Laplacian in each subdomain, so base so-lutions are of the form u = yr1 + r ° , where the constant vector r = [ i - 1 , ! - 0 ] 7 solves Gf = b. (2.5) f = iu + fp , where the particular solution, fp , solves G T G r P = G T b , (2.6) and GrH = 0, (2.7) so in belongs to the nullspace of G. G is m+l-by-2m, so its nullspace has dimension m-i. Including m 4 far-field conditions on u specifies the base solution completely. As a result, stability will depend on global far-field conditions as well as local conditions at T. Note that if b = 0, fp = 0, and only m—2 far-field conditions are required for a stability analysis, as these will determine the base solution to within multiplication by a constant. At worst, this constant may be chosen to have the incorrect sign, and if an instability owing to the base solution is determined, changing the sign of any resulting boundary velocity will be the obvious corrective action. This is the case for the ice-water system, where m = 2 and the global condition is the total energy flux across the boundary. In problems with more than two dependent variables, however, the base solution depends on a linear combination of null vectors of G and local stability (ie at T) is affected by the global constraints (global fluxes) that determine the Chapter 2. Linearised Problem 2 0 correct combination of these null vectors for the problem. The local problem cannot determine these fluxes: they are determined by global conditions, and we consider them to be given. 2.2 Criteria for well-posedness Since we have generalised our boundary conditions, it is not clear which condi-tions G and global fluxes give well-posed problems. A problem is well-posed if a solution exists, that solution is unique, and if the solution depends continuously on the data. We well-posedness spectrally, by driving the boundary conditions by a single wave-mode, in which case the boundary conditions ( 2 . 5 ) become G u = b + efe i a a : . (2.8) We seek a perturbation solution of form ( 2 . 4 ) . Given the form of the driving function, we suppose that rji = helax, and, since A u = 0 , u i = c e ^ ' V 0 " , ( 2 . 9 ) where we choose the appropriate sign for |a | to ensure decay as y —> ±oo . This reduces solving the problem with boundary condition (2.8) to solving for c and h. Thus, we suggest the following definition: When the solution [c, h] exists and is unique for all driving conditions defined by f and all a, we state that the problem with boundary conditions G and given global fluxes is linearly well-posed. Following the development in [ 1 9 ] , we use the asymptotic expansions ( 2 . 3 ) and ( 2 . 4 ) with the assumption that e <C 1 to approximate the conditions that must hold at 770 = 0 given that we wish to satisfy the boundary conditions on n(x). First, we expand u(x,w(x)) as a Taylor series, recalling that rio(x) = 0 . We shall suppress the t-dependence of n for clarity. u(x,n(x)) = u(x,r)0(x) + ern(x) + ...) = ufoO) +eni(x)uy(x,0) + 0(e2). ( 2 . 1 0 ) Now we substitute the asymptotic expansion for u , and our assumed 771 (x), so that at the boundary u(x,r)(x)) = uo(a;,0) +e[u1(x,0) + heiaxuOy(x,0)] + 0{e2). ( 2 . 1 1 ) We also need u „ in order to approximate i i . Now, u n = V u • n, so we need to compute n, the unit normal vector to T. Referring to the sketch in Figure 2 . 1 , we note that the slope rj'(x) = en'1(x) C 1 , so the unit tangent vector is i oc (1,77'(a;)) and so h oc (—n'(x), 1 ) . The ^-component of u „ is 0(e) while the ^/-component is 0 ( 1 ) . Hence, to lowest order, u n = Uj,, so un(a;,J7(a:)) = u0y(x, 0 ) + e[uly(x, 0 ) + heiaxu0yy(x, 0 ) ] + 0(e2). ( 2 . 1 2 ) Chapter 2. Linearised Problem 21 yr+rjW x Figure 2.1: Calculation of the normal to a function y = n(x) showing that when n'(x) <SC 1, i ~ (1,0) and n ~ (0,1). The r](x) sketched here is excessively curved in order to clearly distinguish t and h. In our analysis, however, n(x) is nearly flat. Setting u = [ u n , u ] r , equation (2.8) becomes, to 0(e), G u 0 + cG[ i i i + heiaxp0] = b + efeiax, (2.13) where po = [uoyy, uoy]T is dependent on uo, the base solution, and iio balances the 0(1) terms by equation (2.5). i i i depends linearly on c, so where ui(a; ,0)=e , O I Kc, ±\a\ (2.14) K = 1 (2.15) is a 2m-by-m matrix dependent on the specific elliptic differential operators at hand, the Laplacian in this case. The sign of each ±\a\ term is chosen such that the 0(e) solution decays to zero far from V. Hence, looking again at equation (2.9), this means that — \a\ applies to u above T, and +\a\ applies to u below. Thus, the 0(e) equation becomes G K c + hGpo = f. (2.16) This system is full-rank, hence solvable for c and h, as long as Gpo does not lie in the columnspace of G K . G K is m+l-by-m and is full rank, hence its left Chapter 2. Linearised Problem 22 nullspace is spanned by a single vector w = ke r [ (GK) T ] thus for well-posedness, we require that wTGp0 ^ 0 So for solvability, and (2.17) for each a, where w depends on a. Recall that po depends on the boundary conditions denned by G and b, as well as the given global fluxes. Let us consider two examples. The first will use the boundary conditions for the ice-water model, equations (2.2), a problem we know to be well-posed, and the second will use a set of conditions directly relating the values of the dependent variables at the boundary to their normal derivatives in order to make some connection between our mathematics and physically reasonable problems. T h e ice-water p r o b l e m . The ice-water problem has two dependent vari-ables, one above the boundary and one below, hence we use the operator matrix K = - | a | 0 0 |a | 1 0 0 1 (2.18) Recalling the boundary conditions G for the ice-water model (2.2), G K -\a\ —a\a\ 1 0 0 1 (2.19) so that w = [1, |a | , a |a | ] T . The base solution for the ice-water problem is, as previously stated, represented by the nullspace of G which is [1,1/a, 0 ,0] T . Thus, with the global flux a constant, we may take po = [0,0,1, l / a ] T , and our well-posedness condition is w TGp0 =2\a\ ^0 . (2.20) This is certainly true for \a\ > 0. (The fact that the solvability condition does not hold for a = 0 is not concerning or surprising, as in such case, the perturbation is ef = const with is parallel to the base driving condition b, and may be added to the base driving condition prior to determining the base solutions.) That the solvability condition holds here is an indication that the problem is well-posed, as we should expect from its physical interpretation. Feedback condi t ions . Again, we consider a system having two dependent variables, one above T and the other below. However, let us investigate another set of boundary conditions, 6+ 0 1 0 0 a r 0 1 0 0 1 - 1 (2.21) where again a > 0 is the ratio of thermal conductivities of the upper and lower regions. G again applies to i i — [u+, , u+, u 1 T . W i t h 5+ = ± 1 and Chapter 2. Linearised Problem 23 5~ — ± 1 , the first two conditions are feedback conditions on the variables, where the normal derivative of the solution at the boundary (often interpreted as proportional to a rate of energy flow) depends on the value of the solution. Two feedback types, positive feedback, where an increase in u increases un, and negative feedback, where an increase in u decreases un, can be represented, depending on the signs of 5 ± . Given our convention that the n always points upward, these first two conditions represent negative, or stabilising feedback if 5+ = —1 and 5~ = 1, and represent positive, or destabilising feedback if S+ = 1 and 8~ = — 1. A discussion following the calculation clarifies the physical interpretation and addresses the other two cases. The third condition is relatively arbitrary: in this case we enforce continuity of the solution across V. Leaving the physical argument briefly, let us see what the mathematics says about what the signs of 5^ must be to satisfy the well-posedness condition. We apply the Laplacian in each subdomain, so we use the same operator matrix, K , as for the ice-water example, equation (2.18), so that G K = H<5+ + 1 0 0 a\a\5-+\ 1 - 1 (2.22) and so a\a\S~ + 1 \a\6+ - 1 (a\a\S- + l)(\a\S+ - 1) (2.23) The nullspace of G is [-S+, -5~/a, 1,1] T so p0 = [0,0, -<5+, - ( T / a ] T . Working through the algebra, including identifying (S+)2 = (5~)2 = 1, we have the well-posedness criterion wTGp0 = \a\[(6+ - a5~)\a\ - 2] ^ Ofora ^ 0. (2.24) Hence, requiring that zeros of equation (2.24) lie on the negative real line is sufficient for well-posedness. That is, p{6+,6-) :=5+ -aS~ <0. (2.25) There are four possible combinations of 5+,5~, and we summarise these com-binations in Table 2.1. The results in this table, while a direct result of the mathematics, are con-sistent with the physical interpretation of feedback conditions. Suppose do indeed represent temperature. If temperature decreases, we should expect heat to flow into the boundary. Rate of heat flow out of a region is proportional to — V i i • h*, where n* is the outward unit normal. Hence well-posed boundary conditions should satisfy something like un + u = const. The third case in our table sets this very condition on both u+ and u~, remembering that for u+, the direction of the outward normal to the boundary appears reversed with respect Chapter 2. Linearised Problem 24 6+ S~ comment 1 1 1 - a well-posed for a > 1 1 -1 1 + a never well-posed -1 1 - 1 - a well-posed for all a -1 -1 - 1 + a well-posed for a < 1 Table 2.1: Well-posedness for various feedback conditions. Recall that a > 0. Physical interpretations of these results are in the text. to the region D+, since by our convention, n always points upward. This sit-uation is stabilising, and is called negative feedback. If the opposite condition is imposed, we have positive feedback, the second case in the table, and the problem is ill-posed. The two conditional cases in Table 2.1 are a result of the coupling between the two regions by the continuity condition at the boundary. These cases are well-posed for certain values of a, because if too much energy leaves or enters the boundary from one region for a given u, energy can leave or enter the boundary from the other region to maintain negative feedback overall. Since the degree of this compensation depends on a, so do these conditional cases. In the case of feedback of a single dependent variable at a fixed boundary, a similar analysis shows that our well-posedness condition (2.17) is the same as the requirement that a/9 > 0 for the Robin condition aun + f5u = g to provide a unique solution to the usual elliptic boundary value problem, cf. [19]. 2.3 Generalised Stefan velocities With the linear well-posedness of problems for constant coefficient boundary conditions established, we can discuss the stability of numerical methods for free boundary problems. In particular, motivated by our discussion of stability in the introductory chapter, we wish to compute the decay rate or stiffness of the problem. The stiffness, A, is a function of the wave-modes carried by the boundary, hence, we determine stability by computing the decay rate of perturbations to the boundaries of base solutions of the problem. Recalling that the base boundary is y = n0(x, t) = 0, we set n(x, t) = eqi(x,t) =eh(t)eiax. Our goal is to determine conditions for the decay of the amplitude h(t), dependent on our choice of velocity in our explicit solve. The perturbed problem is linear, so we assume exponential decay, h'(t) = Xh(t), and seek conditions for Re(A) < 0. To satisfy the conditions (2.5) at the boundary, we apply the same calcula-tion as that for determining well-posedness where specifically we choose f = 0, the perturbation now being in the boundary location rather ^than in a driving function at the boundary. Wi th this change, the 0(e) equation (2.16) becomes G K c + hGfio = 0. (2.26) Solving m of these equations determines c. We set the error in the remaining Chapter 2. Linearised Problem 25 equation to equal h'(t), which, to lowest order, is vn, the normal velocity of the boundary. Our freedom now is in which equation, or linear combination of equations, we choose from equations (2.26) as the velocity for adjusting T at each iteration. To lowest order, vn equals the vertical boundary velocity, so dt dt v 1 We specify our velocity choice by the m+l-unit vector v, so h obeys v T G K c + hvTGfj,0 = ~ = \h. (2.28) at For example, if the classical Stefan velocity (1.8) is used for the velocity in the ice-water problem, the velocity vector, v = [1,0,0] T . If instead we choose the error in u+ = u~ as the velocity, the velocity vector is v = [0,1,1]. We examine these particular choices further in the computer examples in Chapter 4. Finally, we isolate A. If the conditions from G used to solve the fixed-boundary problem are orthogonal to v, using orthonormal matrices to manipu-late v T and w T from equation (2.17) to eliminate G K gives A = ^ ^ . (2.29) w J v If the conditions from G used to solve the fixed-boundary problem are not orthogonal to v, we can write the collection of velocity conditions and fixed-boundary solve conditions as M T (2.30) where the m columns of A represent the combination of boundary conditions in G used to compute the fixed-boundary problem. Note that M is full rank. Defining e, = [1 ,0 , . . . , 0 ] T , A = ^ ^ . (2.31) w J M 1 e i Details of this calculation, including the modification required when v is not orthogonal to all of the fixed-boundary conditions, can be found in the appendix. A is a ratio of polynomials in |a | . Since we require that all modes decay, A must not change sign; restricting the numerator and denominator never to equal zero is sufficient to guarantee this. The numerator of (2.29) and (2.31) is exactly the expression (2.17), and, as we recall, depends only on features of the problem, namely the free boundary conditions, defined by G and b, derivatives of the base solution, / J O , which depend again on G as well as on global flux conditions, and on the operator acting on u, the Laplacian in this case, contained in K and hence in w. Requiring that the numerator never equal zero is the same as condition (2.17), a satisfying consistency in our analyses. Chapter 2. Linearised Problem 26 The denominator of A depends on the choice made for the velocity. Since the sign of A(|a|) indicates the stability of the method with the chosen veloc-ity, equation (2.29) can be used to evaluate the effects of our velocity choice independently from the underlying problem. For a well-behaved denominator, one that does not change sign, we must choose a good scheme, which amounts to choosing a good velocity. Moreover, v affects the order of the denominator and, hence, the order of growth of A. We now discuss the effects of our velocity choice. 2.4 Good velocity choices We are now in a position to present general rules for choosing good velocities, that is, those that are single-signed giving A = 0(1) in \a\. Since the numerator of A is fixed by the underlying problem, we reduce our discussion to asking when the denominator of equation (2.29) is of highest order. At the very least, this minimises the growth rate of A(|a|). We can accomplish this by determining the elements of w and choosing v to select its highest order elements. Consider the calculation of w. Since the boundary conditions are linearly independent, we can row-reduce G to be upper triangular without affecting the problem. The bottom row of such row-reduced G is a pure Dirichlet condition. We write this row-reduced G as a sum of its Neumann and Dirichlet components, G = [N,D] = G N + G D , G w = [ N , 0 ] , G c = [0 ,D] , (2.32) and do likewise for K : KN = K = K J V + Kj), ±\a\I 0 (2.33) where I is the m-by-m identity matrix. (We have abused notation somewhat here, in that elements in the upper half of K may be variously +\a\ or — \a\.) With these definitions, we observe that GN~KD = G ^ K j v = 0, so G K = G N K N + G D K D = ± | a | N + D (2.34) (again with an abuse of notation for ± | a | ) , where not only is N upper triangular, but, because G is row-reduced, the bottom row of N is zero. W i t h this in mind, we see that G K has at most m elements of 0( |a | ) , none of these elements appearing in its bottom row. Thus, its left nullspace, w, which by definition is orthogonal to its columnspace, has elements of at most 0( |a | ) in all but its bottom row, where it is 0 ( | a | 9 ) , 1 < q < m; q is the number of different C( |a | ) terms in G K . We see that the numerator of A is in general 0 ( | a | ' ) , and its denominator is also 0 ( | a | 9 ) as long as v = [0,0,1] T . We revisit our two Chapter 2. Linearised Problem 27 examples, the ice-water problem, and the feedback problem, to illustrate this calculation. Recall that for the ice-water problem, G = 1 - a 0 0 0 0 1 0 0 0 0 1 (2.2) and w = [1, |a | , a |a | ] T . We see the highest-order elements of w in the positions corresponding to the rows of G representing pure Dirichlet conditions. A choice for v that maximises the denominator of A will be non-zero in at least one of its last two elements. Note that ice-water problem is somewhat unique among problems of two dependent variables in that it has not one, but two pure Dirich-let conditions at the free boundary. This means that there are two highest-order elements of G , and hence two basis directions for v that maximise the order of the denominator. Although one consequence of this is that the denomina-tor can be at most 0( |a | ) rather than 0 ( | a | 2 ) , as is the case for the general m = 2 problem, so, too, is the numerator limited to 0( |a | ) , the order of the highest-order elements of w. For example, if we calculate the stiffness of the trial method using the Neu-mann and Dirichlet velocities introduced on page 11, we find that A = 2|a| for the Neumann velocity, where v = [1,0,0] T , and 2y/2 A = 1 + a (2.35) (2.36) for the Dirichlet velocity, where (normalised) v = [0, l / \ / 2 , l / \ / 2 ] T - Indeed, the Neumann velocity is 0( |a | ) and the Dirichlet velocity is 0(1). Again, the effects of the wave-number dependence of the stiffness for each of these v will be seen in the computer experiments in Chapter 4. For a second example, consider the unconditionally well-posed (negative) feedback problem, where S+ = — 1 and S~ = 1, so that G = - 1 0 1 0 0 a 0 1 0 0 1 - 1 (2.37) and a\a\ + 1 w = | - ( H + l) . (2.38) - ( | a | + l ) ( a | a | + l) . The highest order element of w again corresponds to the problem's Dirichlet condition in the third row of G , and choosing a velocity such that v is non-zero in its third element is enough to ensure that the denominator is the same order as the numerator in \a\ and, hence, that A(|a|) = 0(1) for large \a\. Chapter 2. Linearised Problem 28 Thus, regardless of the possibility that there may be more than one pure Dirichlet condition represented in G , the preceding argument and example in-dicate that using Dirichlet conditions from the row-reduced G as the velocity gives a A that is 0(1) in |a | . This is a general rule: In order to guarantee that A be 0(1), a Dirichlet condition must be included in the velocity. A similar analysis shows that choosing pure Neumann conditions for the ve-locity gives an O d a l 9 - 1 ) denominator, 1 < q < m as defined following (2.34), whereupon A = 0( |a | ) . For a fixed step size, the method becomes unstable for some sufficiently large wave-number, or, equivalently, for some sufficiently fine discretisation. Indeed, the effect velocity choice has on stability manifests itself in previous work by Karkkainen and Tiihonen, where one of their (shape optimisation) algorithms fails to converge for certain fine discretisations [18]. Moreover, this behaviour is consistent with the well-known gain of regular-ity incurred in deducing Dirichlet boundary data from Neumann data. The good behaviour of Dirichlet velocities agrees with the converse result, that the Neumann-to-Dirichlet boundary map tends to improve data regularity. Regular-ity relationships between boundary data are well-known in the inverse problems community, cf. [29], a direct consequence of the trace theorem for functions in Sobolev spaces, cf. [2, ch 7]. Garabedian's approach to the Riabouchinsky flow problem (1.3o-1.3rf) is to use dnu + KU = 1 to solve the fixed-boundary problem. If we replace the differential equation (1.3a) with Aip(x,r) = 0, and linearise around r = R = 1 (this changes the LHS of (1.3d) to 1), we have a paradigm problem similar to that treated by Garabedian that fits our analysis. In Garabedian's approach, K is the curvature (not the stiffness), and can be considered constant to first order in our linearised treatment. (Garabedian makes this choice because this is a boundary condition that does not change under small normal perturbations of T.) We set the error in u (actually, — u is the correct direction) as the velocity. The problem is one-sided (m = 1) so we solve for u in the half-plane r > 1 and take K = [ - | a | , l ] T . G is the identity matrix, b = [ -1 ,0] T (the negative sign is a result of our choice of h as the upward, not outward, normal). Note that m = 1 is a special case for which G has no nullspace. Instead, we must have nonzero b so that there is a nonzero particular solution, ip. In this example, f = rp = b, giving po = [0,1] T. Garabedian's approach amounts to M = 0 - 1 1 K (2.39) Using equation (2.31), we have A = + H (2.40) and in Garabedian's method, A is bounded for large wave-numbers. Instead, if we consider the original approach, where u = 0 is applied in the fixed boundary Chapter 2. Linearised Problem 29 solve and the error in un is used to adjust the boundary, M = 1 0 0 1 (2.41) and A = - | a | . (2.42) A is unbounded with \ct\ for this standard velocity choice. If we use the Garabe-dian velocity on our linearised problem, the curvature of the base boundary is zero, so A = — 1 for all a. Referring back to our convergence discussion on page 14, we see that this gives a growth rate G = 0 with an optimal k = 1. This is exactly the value of k that Garabedian successfully uses in his method. One point remains: while the numerator of A will not change sign with \a\ as long as the well-posedness condition in equation (2.17) is satisfied, it is possible that choices of v give A that change sign for \a\ > 0. We must be wary of this possibility in choosing v. A random choice for v , even one that maximises the denominator of A, can destabilise the solve of an otherwise well-posed problem. It is possible to choose v giving Re(A) > 0 for all \a\. However, these v are straightforward to correct by choosing —v for the problem. Our concern, then, is limited to those v that cause a change sign in w T v , regardless whether the sign change is from negative to positive or positive to negative. We define v to be workable if w T v does not change sign with \a\. W i t h this caveat, we can construct plots of the behaviour of A for each choice of v. For example, if our problem has two dependent variables (m = 2), v is a 3-vector, and, normalised, sits on the unit sphere in E 3 . Figures 2.2 and 2.3 plot the effect of velocity choice on the workability and stability of solutions to the ice-water system and unconditionally well-posed feedback problem, respectively. The axes for these plots represent the radial and azimuthal angles of v and the area shading of the plots indicates the effect of v on each problem. The velocity choice plots show that in the case of the ice-water problem, some combination of the pure Dirichlet velocities, v = a[0,0, ± 1 ] T + 6[0, ± 1 , 0 ] T , gives a method that is 0(1) and workable. In fact, the only workable 0(|c*|) velocity is the classical Stefan velocity! In the case of the negative feedback problem, the method using the pure Dirichlet velocity, v = [0,0, ± 1 ] T , again gives an 0(1) method that is workable, as do combinations of the Dirichlet and Neumann conditions. While there are relatively few unworkable velocities for the negative feedback problem, identifying which velocities to avoid would be difficult without the linear per-turbation theory. In Garabedian's example, represented by equation (2.40), it happens that the curvature K is always positive, so his method is workable. In all three cases presented here, the ice-water problem, the negative feed-back problem, and Garabedian's example, the pure Dirichlet velocity gives a workable method having 0(1) stiffness. This suggests the conjecture that for well-posed problems, choosing the pure Dirichlet condition for the velocity al-ways gives a workable 0(1) trial method. Such a result would be especially useful for cases where m > 2, and the construction of two-dimensional stability maps such as figures 2.2 and 2.3 is not possible. Unfortunately, this conjecture Chapter 2. Linearised Problem 30 •&nf2 71/2 Tt G 371/2 2TI Figure 2.2: Effect of v on the ice-water problem, o = 10. White areas represent unworkable v, grey areas represent workable v where A = 0(1), and the velocities v = [ ± 1 , 0 , 0 ] T , not distinguishable at this plot resolution, are workable v where A = 0( |a | ) . ^ T t / 2 0 TC/2 7c 3TC/2 2n e Figure 2.3: Effect of v on the unconditionally well-posed negative feedback problem, a = 10. White areas represent unworkable v, grey areas represent workable v where A = 0(1), and black areas represent workable v where A = 0(|cv|). Chapter 2. Linearised Problem 31 is at this time unproven, and with the results presented here, we can only test given velocity choices for their workability and effect on stiffness. The treatment of the free boundary problem thus far has assumed infinite domains where V may support an infinite collection of wave-modes. In a prac-tical numerical solve, however, T will be represented by a set of points on a finite domain. While our analysis of stiffness and stability is indicative of the behaviour of a practical scheme, before proceeding to the computer examples, we must analyse stability in the context of the discretised problem. Fortunately, with some adaptation, our work on the stability of the infinite linearised problem applies. Chapter 3. Analysis of the Implementation 32 Chapter 3 Analysis of the Implementation Although velocity choice is our primary point of discussion, an implementation of the explicit scheme requires further treatment of the boundary in the solve. T is represented by discrete points, with a spline fit to represent the continu-ous boundary. This lends some independence between the numerics of the free boundary adjustment and the numerics of the P D E solves in D+ and D~. In this way, we analyse a semi-discrete scheme in which we consider the P D E to be solved exactly. In addition, we smooth the data used in computing the local velocity of T at each iteration by averaging the P D E solutions (and their deriva-tives normal to F) by appropriate windowing functions (also called mollifying functions) along the arc length of the spline. Both of these operations, spline fitting and smoothing by windowing, can be treated as signal filtering processes. Windowing, especially, is an integral part of signal analysis, and is treated in the engineering literature as low-pass filtering (low-frequency signal selection). In light of existing tools that treat signals in the frequency domain, we make extensive use of Fourier transform space in the succeeding analysis, consistent with our frequency-dependent treatment of stiffness. This section both i) verifies that the spline fit and smoothing operations do not destabilise schemes deemed otherwise stable by the preceding linear stability analysis; and ii) shows that the stiffness predictions from the idealised linear analysis extend to a practical implementation having a discretised boundary. Our domain is now 27r-periodic in x. 3.1 Fourier transform of the cubic spline Weighing the costs of implementation complexity and the effective benefits, we choose' a piecewise cubic polynomial spline fit to approximate F. Since our linear stability criterion (2.29) assumes an infinite spectrum of modes carried by F, we first seek a relationship between the necessarily finite modes carried by the discretisation of the free boundary, and the modes carried by the spline fit of the discretisation. We continue to assume that T may be represented as a discretisation of a function y = n(x) However, n(x) is now a 27r-periodic function in x. We discretise y — n(x) by N points (xj,yj) evenly spaced by Aa; = Xj+i —Xj Chapter 3. Analysis of the Implementation 33 (3.1) so that its 27r-periodic cubic spline interpolation is y(x) = Y,fjo' y3{x)x[xi,xi+l)(x), where X[Xj,xj+1)(x) is the characteristic function for the interval [XJ,XJ+±), and *i = ^ Bi = l- Aj Cj=A3j-Aj, Dj=B]-Bi. This representation of the cubic spline fit is taken from [22]. For a fixed discreti-sation, [XJ], y(x) depends linearly on the parameters [yj] and [y'J]. However, [y'J] is the solution to the linear system E[j/"] = F[j/j], where, for periodic T, _ Xj + i—X ^3 — Ax ' (3.2) E = F = 4 1 1 4 1 - 2 1 1 (3.3) (3.4) are N-by-N matrices, so y(x) depends linearly on the data [yj]. We wish to compute the Fourier transform of y(x). Since this transform is also linear, we represent the data as a sum of its discrete Fourier modes, N-l yj = /Ty(n)wN, 3=0 (3.5) where W^f := exp(i^-nj) is the kernel of the discrete Fourier transform. We compute the Fourier transform of y(x) with the specific data Sj = W^J and sum the results over all modes, n. Now, since Sj are also eigenvectors of E , bfF, we can assume y'J bW^fJ, where b is determined from the linear solve, E[bW#] = F[W#], N (Wy+4 + WR)[bW^] = (W^ cos(2im/N) 2 + W^)[W^j N J ' b = (3.6) cos(2im/N) + 2 We wish to compute the Fourier transform of the 27r-periodic spline fit y{x). From the variety of definitions for the Fourier transform, we choose y{n) = ~ j 2 J y{x)e-inxdx. (3.7) Chapter 3. Analysis of the Implementation 34 The interpolation (3.1) is C2 in its 27r-periodic domain, so we can integrate by parts, observing that y(0) = y(2w), y'(0) = y'(27r), and y"(0) = J / " (2TT ) . Now, ^ = ^-Tv I*'y'"We-inXdx- (3-8) 27r (in)* J0 Furthermore, y(x) is C°° a.e., that is, except at the finitely many data points. So with convergence of the Fourier series for y(x), we compute its continuous Fourier transform piecewise: ^ ^ h j k ^ C 1 ^ 6 ^ ' (3-9) V ' j=0 JX3 where from the definition of the cubic spline, » W = ( £ j s ( W i + i - » i ) - ( 3 - 1 0 ) Wi th our pure Fourier data, substitute y" = b W^j1, so yi{x) = (iwm"1} w" = const' (3'n) Evaluation of the integrals is now trivial, and, after some manipulation, the Fourier transform of the cubic spline fit of a pure Fourier mode is , s i i~\T\ 3 sin4(7rn/AT) sN(n) = s(n/N) = — . . ' . . (3.12) v ' v 1 ' cos(2im/N) + 2 (nn/N)'1 v ' Details of the calculation of the Fourier integral can be found in Appendix A.3 . Given the Fourier representation of the data as in (3.5), we can write the Fourier transform of the cubic spline fit of general data as -t \ 3 sin4(7rn/A^) y(n> = 7^ 1An , o~7 TTJ^y(n> ~ SN{n)y(n), (3.13) cos(27rn/A') + 2 (im/Np and the behaviour of sjy(n) indicates how high-order modes not present in the discrete data appear in the continuous cubic spline fit. This relationship assumes that n takes all integer values, yet y(n) is defined only for 0 < n < N — 1. Thus, we make sense of (3.13) by noting that WjJ1 is A^-periodic, so we define y(n — rN) = y(n) for integer r. Equation (3.13) shows that stability depends on the response of our choice of generalised Stefan velocity to all modes in y(x), not just the low-order modes in the discretisation [yj]. Three qualities of sjv(n) are immediately apparent: i) for n —> oo, sjv(n) ~ 0 ( n ~ 4 ) ; ii) s^r(n) is non-negative; and iii) the argument of SN(TI) is the ratio of the mode number to the number of points in the discretisa-tion, and does not depend on n explicitly. The first quality is a consequence of the degree of smoothness of the piecewise cubic fit, and serves as a check of the calculation; the remaining two features are important for the stability analysis that follows. Figure 3.1 gives a sketch of sjv(n). Chapter 3. Analysis of the Implementation 35 s 0.5 0 i i i i i 0 0.5 1 1.5 2 n/N Figure 3.1: A sketch of Sjv(n) = s(n/N), showing its non-negativity and rapid decay. Siv(n) is even, so only n/N > 0 is shown. 3.2 The cubic spline in generalised Stefan schemes Using the Fourier transform of the cubic spline fit, we identify the continuous spectrum in \a\ from Chapter 2 with our discrete spectrum in n, and seek a decay rate analogous to A(|a|) that indicates the stability of the scheme. At each iteration in the computation, we determine the error in the residual boundary condition at each Xj, and adjust yj proportional to this error. As determined by the linear stability analysis, the adjustment to each point of T is proportional to A(|a|) for each pure wave-mode, so on our 27r-periodic domain, the total adjustment is the sum over integer modes, oo y(x) = Y,. Kn)y(n)einx. .(3.14) n = —oo The notation y(x) is suggestive of the velocity of the boundary. y(n) is again the continuous Fourier transform of the spline fit. As we are interested only in the adjustments of the points (xj,yj), we recognise Xj = jLj and rewrite (3.14) so oo N-1 y(xj)= Y ~ rN)y{n - rN)WJN(n~rN). (3.15) r = — oo n=0 Next, using (3.13), we substitute for the Fourier transform of the spline fit. Recalling that at worst, A(|a|) = C( |a | ) , and s^ = 0(n~4), we are guaranteed Chapter 3. Analysis of the Implementation 36 uniform convergence, and can exchange the sums. Wf}n r is iV-periodic, so J V - 1 / OO \ y(xj) = I A ( n -rN)sN(n - rN)y(n - rN) Wtf. (3.16) n = 0 \ r = —oo / We see that the expression in brackets is the D F T of [y(xj)]. By our definition, y(n) is iV-periodic, hence y(n) = 7jv(n)j7(«), (3.17) where oo 7jv(n) = Kn-rN)s{n/N -r), (3.18) r= — oo reverting to notation for s that reflects its single argument, n/N. Just as A(|c*|) does for continuous spectrum T, 7JV(H ) indicates the stability of the solve with spline fit data having ./V points in its discretisation. We can use equation (3.18) to assess the stability of our scheme. If A ~ C( | a | « ) , then X(n-rN) ~ AN^n/N -r\q, (3.19) for some constant, A, and oo •yN(n)~ANq ^ \n/N - r\9s{n/N - r). (3.20) r== —oo Here we observe the benefits of the spline fit in s(n/N). s > 0 for all n, so if our Stefan velocity gives A < 0, 7 < 0 also. Since s(n/N) decays as 0(n~4), the sum only converges as long as q < 3. However, for our ice-water example, q < 1, so 7jv(n) ~ 0(Nq), just as A(|a|) ~ 0 ( | a | « ) . The behaviour of the spline fit implementation coincides with that of the problem in the idealised setting. 3.3 Smoothing Stefan velocity data Although we now have rules for generalised Stefan schemes that will work in practise, we have neglected an important element of the method: the P D E solve of the domains D± that occurs at each time step introduces error to the velocity computation. The velocity computed directly from the P D E solution (and its derivatives normal to T) has an error according to the P D E solver used. One way of attenuating this error is to smooth solution data along the boundary by convolution of that data with a windowing function. (This technique is also called mollification.) Wi th the framework of the spline fit implementation, we are in a position to incorporate smoothing into our analysis. Convolution in real space becomes multiplication in Fourier space. Thus, with the appropriate modification, the velocity in (3.14) is 0 0 y(x)=2rr £ X(n)f(n)y(n)einx. (3.21) 71= — OO Chapter 3. Analysis of the Implementation 37 / is the Fourier transform of the windowing function, and the factor of 2w comes from our definition of the Fourier transform (3.7) applied to convolution. Wi th this modification, oo 7 Ar(n)=27r Y Hn-rN)f(n-rN)s(n/N -r). (3.22) r= — oo We begin discussion of appropriate windowing functions with an example. Although not computationally efficient, an example windowing function is the 27r-periodic Gaussian, r = —oo v v ' normalised such that the area under f(x) on [0, 2TT] is 1. o controls the degree of smoothing by f(x). To see the effect of varying a, we compute the Fourier transform / according to our convention (3.7): /(n) = ^ e x p ( - ^ ) , (3.24) whereupon we see the usual spreading of data in the frequency domain as f(x) narrows in real space. Details of this calculation, including a general result for computing the Fourier transform of periodic functions, can be found in Appendix A.4. We may consider the velocity data provided by the P D E solver to be the sum of an exact solution and noise. The spectrum of the noise depends on the P D E solver used, and the amplitude of the noise depends on the order error of the solver. While we cannot eliminate this error altogether, we see that the effect of the windowing function is to attenuate the high-frequency error. Of course, the smoothing function also attenuates high-frequency components of the solution; an ideal smoothing function attenuates only modes having error and solution values of the same magnitude. We control the degree of smoothing by varying a. In real space, a represents the width in x over which the window averages data; in Fourier space, a controls the largest frequencies retained by the scheme. Thus, we expect that low-frequency modes will remain unaffected by the smoothing. If we fix <r, in turn smoothing over a fixed x-interval, as we increase N, an ever smaller fraction of modes, n with respect to N, are left unaffected by smoothing, and we see dissipation of the fine detail we may wish to capture by a fine discretisation. Instead, we can set a = a Ax = a2^-, affecting smoothing over approximately a data points. Now, / ( n -rN) = 1j- exp (-2ir2a2(n/N - r)2) , (3.25) and a constant fraction of low-frequency modes are relieved of dissipation, such that we can recover the fine boundary detail we desire. Here we explicitly see the trade-off between consistency and stability endemic to smoothing techniques. Chapter 3. Analysis of the Implementation 38 Wi th this windowing function, we have oo lN(n) = A ( n - r A 0 e x p ( - 2 7 r V ( n / A T -r)2)s{n/N -r). (3.26) r= — oo The incorporation of this rapidly-decaying smoothing function does more than damp error from the P D E solve; it also guarantees convergence of the sum. As previously mentioned, for the ice-water example, A does not grow quickly enough relative to the cubic spline decay for convergence of the sum to be an issue. However, for higher-order problems (problems with many dependent variables, where m > 2) or lower-order splines (ie piecewise linear) applying smoothing can force convergence of 7jv(n) and aid stability. In addition to being easy to work with analytically, this example window-ing function has several good practical features. Most importantly, its Fourier transform is everywhere positive. Hence, smoothing by windowing here cannot destabilise the solve by introducing sign changes in T J V ( T I ) . f(x) has unit mass, so limy \ \v-f*v\\ L r = 0, v£Lp, (3.27) where * denotes convolution. f(n) is everywhere real, indicative of the even symmetry of f(x): when averaging about a point, data to the left of that point should weigh equally with data to the right. f(x) is periodic, consistent with the problem layout. Finally, the degree of smoothing, hence the error introduced, is easy to control. However, our choice for f(x) falls short in that it does not have compact support as is often desirable in computation. Furthermore, it may be desirable to choose smoothing functions tailored to the spectrum of the error produced by the specific P D E solver employed. In addition, we have chosen to smooth over x rather than boundary arc length. Because the perturbation is small in our analysis, the two metrics are nearly identical, but this will not be true in general and may be important for highly irregular boundaries. A l l the same, we ignore these shortcomings and use the Gaussian window to smooth over x in the example computations. For other windowing functions, the reader is invited to consult [16] and like texts on signal analysis. Theoretical treatments of windowing (mollifying) functions appear in many numerical texts, cf. [20]. Wi th our practical scheme established, we proceed to the computer exam-ples. Chapter 4. Computer Examples 39 Chapter 4 Computer Examples We now demonstrate the effect of velocity choice on the ice-water model system of Figure 1.7. For this problem, we use the free boundary conditions (2.2) for G , the differential operator matrix K representing the Laplacian, and adjust the boundary at each iteration by two different choices of velocities, v. In particular, we wish to compare the velocities presented on page 11, a Neumann velocity (the classical Stefan velocity, in this case) and a velocity containing Dirichlet components. That is, we compare v = [1,0,0] T , which we refer to as the Neumann method, and v = [0, l / \ / 2 , l / \ / 2 ] T , which we refer to as the Dirichlet method. As we have already investigated, there exist choices for v that give unworkable methods. However, looking again at the map in Figure 2.2, we note that our velocity choices do not make for unworkable schemes. Our domain is the rectangle —it < x < w, —1 < y < 1, with periodic boundary conditions at x = ±n. Our boundary conditions are u(x, —1) = —1 on the bottom boundary, and f - —- < r < -u(x 1) = < ? 2 - ^ - 2 (41) v ' ' (I elsewhere We set a = 10 and our initial T as y = 0, discretised and cubic spline fit by N points [(xj,yj)]. To solve the fixed boundary problem at each iteration, we use the two re-maining free boundary conditions orthogonal to v. For the Neumann method, we impose u+ = u~ = 0, and for the Dirichlet method, we impose i t + = u~ and u% — au~ = 0. Wi th these conditions imposed, at each iteration, we solve A n = 0 in each subdomain using the commercial finite element method (FEM) package F E M L A B , version 2.2 produced by Comsol Inc [17]. The F E M solve uses first-order (piecewise linear) Lagrange elements on mesh elements having diameter D « 0.3 near y = ± 1 and D « 2ir/3N near V. This ensures that to within the resolution of the boundary, our trial method uses a good approx-imation to the exact solution of Au^ = 0 near T. We again emphasise the simplicity of this implementation: the use of this particular F E M package to solve the fixed-boundary problem is independent of our method of adjusting the boundary. Very little original code had to be written in the implementation of this scheme. Following each F E M solve, we compute the residual v(x) in the appropriate third boundary condition along the boundary, the residual in u+ — au~ for the Neumann method, and the residual in (1/2)(u+ + u~) for the Dirichlet method. We smooth the residual by convolution of v(x) with respect to x by Chapter 4. Computer Examples 40 the normalised 27r-periodic Gaussian function of our example, equation (3.23). We choose o = n/N, so that this convolution computes the smoothed residual Vj(x) at each boundary point (xj,y.j) by smoothing over approximately b — 2 data points to either side of each point in question. The algorithm moves each boundary point normally to T, so where k is the time step size, and returns a cubic spline fit of this revised T to the succeeding F E M solve. Note that our implementation deviates from our analysis on this point, as the boundary is moved only in the y-direction in the analysis. The implementation is consistent with the more general case where T is not nearly flat. Wi th the boundary discretised and spline fit by N = 50 points, the largest stable time step for the Neumann method is k = 0.01, and for the Dirichlet method is k = 0.2. The smaller time steps required by the Neumann method corresponds to the stiffness predictions for the ice-water problem on page 27. Figures 4.1 and 4.2 for the Neumann and Dirichlet methods, respectively, illus-trate the progress of these solves showing u in each subdomain, as well as the development of the free boundary, T. Note that each method solves a different problem, but converges to the same steady-state boundary. Furthermore, the initial boundary for each case is y = 0, a poor approximation to the solution, yet both methods converge. The effects of the stiffness of each method can be seen as N is increased, which we observe by comparing our numerical solution with the exact solution. (Indeed, with this problem, we are in the position of having the exact solution to the free boundary problem at hand. Essentially, by setting v(x, - 1 ) = au(x, - 1 ) and leaving v(x, 1) = u(x, 1), we solve At; — 0 in the entire domain. This change of variables in the lower subdomain ensures that the flux is continuous everywhere in the domain, and, thus, the boundary is represented by the level set v = 0. A closed-form solution exists for this problem for v, and v = 0 can be written as an implicit relationship between x and y. Details of this calculation appear in Appendix A.5.) Figure 4.3 compares the Neumann and Dirichlet methods for a variety of N. The plots show the mean .square error in the boundary points, ,n+l ,y?+l) = (x?,y?) + kvjn, (4.2) N-1 (4.3) at each iteration over a range of N taken from the set N G {20,50,100,200,500}. (4.4) Time steps k from the set k e {1,0.5,0.2,0.1,0.05,0.02,0.01,...} (4.5) Chapter 4. Computer Examples 41 iter = 1 iter = 2 iter = 4 -4 -2 0 2 4 -4 -2 0 2 4 -4 -2 0 2 4 iter = 8 iter = 15 iter = 30 -4 -2 0 2 4 -4 -2 0 2 4 -4 -2 0 2 4 Figure 4.1: Evolution of the ice-water boundary by the Neumann method. The black line represents T while the shaded data represents u. The initial boundary is y — 0. were tested for each N and the error data for the stable method having the largest k for each given N is reported in the figure. Using the Dirichlet method, the same k can be used for each choice of iV, whereas with the Neumann method, to maintain stability k must be decreased as N is increased, consistent with our analysis. As a result of this need to increase k, the number of iterations required to converge to the solution increases also. The maximum stable step sizes for each N are given in Table 4.1. We note here that even at (mean) steady state, there is oscillation present in the Neumann method not present in the Dirichlet method. Figure 4.4 shows this oscillation in detail for N = 10,20,50. We see a reason for this oscillation by looking at the linearised velocity for our generalised Stefan schemes. Using equation (2.28), we have vn = -|a|[l,o]c (4.6) Chapter 4. Computer Examples 42 iter = 1 iter = 2 iter = 4 -4 -2 0 2 4 -4 -2 0 2 4 -4 -2 0 2 4 iter = 8 iter = 15 iter = 30 -4 -2 0 2 4 -4 -2 0 2 4 -4 -2 0 2 4 Figure 4.2: Evolution of the ice-water boundary by the Dirichlet method. Again, the black line represents T while the shaded data represents u. Although difficult to make out, the initial boundary is again y = o. for the Neumann method, versus ^ = - - L [ l , l } c + ;(l + i)ft (4.7) for the Dirichlet method. The finite element method incurs high-order wave-modes in the solution, and although the magnitude of these wave-modes is small, the linear dependence on |a | of the Neumann method puts this method close to the limit of stability, resulting in oscillations. The Dirichlet method does not have linear |a|-dependence and hence does not suffer these oscillations. Finally, we observe that although the error tends to zero as N —> oo, we are ultimately limited by the accuracy of the F E M solve. The mean square error decays as 0(N~ll2), supposedly consistent with accuracy estimates for our finite element method. These estimates are future work. Accuracy as a function of N is shown in Figure 4.5. The accuracies of the two methods are comparable. However, for N = 500, while the Dirichlet method took approximately 3 hours Chapter 4. Computer Examples 43 iterations Figure 4.3: Plots of mean square error versus iteration number, showing the effect of N on the maximum size of k for each method. Solid lines represent the Neumann method; dashed lines represent the Dirichlet method. The oscillations for the Neumann method are discussed in the text. N maximum k Neumann method Dirichlet method 20 0.02 0.2 50 0.01 0.2 100 0.005 0.2 200 0.002 0.2 500 0.001 0.2 Table 4.1: Comparison of the Neumann and Dirichlet methods by the maximum stable step sizes k for each discretisation N. Chapter 4. Computer Examples 44 1 1 1 • i r — - N=10 -N=20 - - - N=50 i\ 1 , v . v ^ v : v . ~ : ••- /•v 7 — v y/ "" / i s I \ V r ' '• ' V ' i \ ' ' VAi; \/V! \; • \ - T - - - - i - ' i ' j ' ! ' ; • / , • > • • • • ' • , • V v v 'f •; ) i i i,' \ ' V V '? v v V 1 i 0->l 1 1 1 1 50 60 70 80 90 100 iterations Figure 4.4: Close view of the oscillations in the solution by the Neumann method. Unlike Figure 4.3, the horizontal axis is linear. to solve for the free boundary, the Neumann method took approximately 6 days running on the same machine. This observation is encouraging justification for using generalised Stefan velocities for solves where the governing P D E s might be more complicated. The work culminating in these computer examples constitutes our main body of study. The central idea, that making good velocity choices can lead to superior numerical methods, has been verified in this example, and'the way is now paved to apply these ideas in more general settings. This will be the course of the next and final chapter. Chapter 4. Computer Examples 45 Figure 4.5: Comparison of the accuracy of the Neumann and Dirichlet methods for increasing iV with the k in Table 4.1. Chapter 5. Continuing Work 46 Chapter 5 Continuing Work While one immediate goal of this project is to apply a generalised Stefan scheme to the porous media phase change problem presented in the introductory chap-ter, other work continues. The first part of this chapter discuss a level set method for free boundary adjustment that retains the good stiffness charac-teristics of the generalised Stefan velocity method. Level set methods have an advantage over the front tracking method employed thus far in that they do not require a priori knowledge of the free boundary topology. In Section 5.1, we treat a level set method with the flat boundary used in Chapter 2, then in Section 5.2, we consider both the front tracking and level set methods applied to more general boundaries. Section 5.3 concerns the well-posedness of boundary data near points where the free boundary contacts the physical boundary of the problem, an issue we have thus far avoided through the use of infinite and periodic domains. Finally, Section 5.4 summarises this work. Work presented in this chapter is motivated by the well-posedness of generalised Stefan problems and the workability of their numerical solutions. 5.1 Generalised Stefan velocities for a level set method Our analysis has thus far assumed that one coordinate of our free boundary can represented as a function of another coordinate, that is, for example, that in R 2, the boundary location is described by the points (x, r)(x)). This will not be true in general, such as in computations on bubble formation and Hele-Shaw flow, cf. [19, ch 7], where the free boundaries (surfaces, in IR3) enclose a region of the domain. Computing such free boundaries can be challenging, primarily because the topology of the enclosed regions can change, for example, when two bubbles merge, or a single bubble divides. In the case of steady-state solves, we must be concerned with making accurate initial guesses for the steady-state topology of the free boundary. In three (or more) dimensions, the front tracking technique used thus far for our trial method can be particularly difficult to implement stably, and so we present an alternative. A successful strategy for dealing with free boundary problems, particularly those in which the free boundary topology changes over the course of the solve, is the level set method of Osher and Sethian. Their method was first presented in [21] for addressing curvature-driven boundary evolution, yet has been further developed, whereupon Sethian's text [23] treats the extension of the level set Chapter 5. Continuing Work 4 7 idea to other techniques and applications. Instead of tracking the free boundary, the level set method solves for a function, </>(x, r), that typically represents the signed distance from any point x in the domain to the boundary. Hence, the level set cp — 0 represents the boundary location. <f> obeys an evolution equation that depends on the velocity of the boundary. Our goal in this section is to show that the evolution equation for <f> inherits the desirable stiffness properties obtained by making good velocity choices for the corresponding front tracking method. Osher and Sethian treat the time-dependent problem of front propagation in W1 as an initial value problem for a function 0 : [0, oo) x I " - f E. Again, the zero level set of <j> represents the (moving) free boundary. We develop an evolution equation for <p beginning with 0(x(*),t)=O, (5.1) the set x(t) representing the free boundary. The -^dependence of x represents the general tendency of the free boundary to change location with time. By the chain rule, & + V0(x ( t ) , t ) -x'(t) = 0, x e r (5.2) Now, the normal velocity of the boundary, vn = x'(i) • ii. The gradient of a function is normal to its level sets, so n = V 0 / | V 0 | . Thus, <j> obeys / <j>t+vn\V4>\ = 0, x 6 D , . \ 0(x,t = O)=/(x), X G U • ^ Note that this level set method is a capturing method, and we have extended its evolution equation to the entire domain D, and thus removed the explicit dependence of x on t to reflect this. vn may depend on any combination of local, global and intrinsic properties of the system: for our development, we take vn to be a generalised Stefan velocity of our trial method. This means that time t in (5.3) represents artificial iteration time, just as throughout our development of the front tracking trial method. Since <f> must be solved for all x, not just for those x on T, and that vn must also be known for all x. Extending vn to the entire domain from T may not be trivial, and we shall revisit this point at the end of this section. <f> is not unique for a particular T: there are many functions that share the same zero level set while satisfying the evolution equation (5.3). This means that any number of initial /(x) are permissible for the same free boundary problem. It is often desirable, however, to specify <f> as a signed distance function, where <j> is +1 times the distance from the free boundary on one side of T and —1 times the distance from the free boundary on the other side. With this framework, we suggest an algorithm for using the level set method with generalised Stefan velocities to solve the steady state free boundary prob-lem. 1. Choose an initial level set function. Its zero level set is the initial approx-imation to the free boundary. Chapter 5. Continuing Work 48 2. Solve the resulting fixed boundary P D E using all but one of the free bound-ary conditions. 3. The remaining condition will not be satisfied unless the position of the free boundary is accurate: the residual in the remaining condition is the boundary velocity. 4. Extend the boundary velocity to the entire domain, and step the level set equation (5.3) forward. If the adjustment of the level set function is small, stop; otherwise, adopt the new zero level set as the boundary and return to step 2. The level set method is touted for its ability to evolve a free boundary without knowledge of its location at every step. It may seem as though we nullify this particular advantage in that we must know the boundary location at each iter-ation in order to solve the fixed boundary P D E compute vn. However, unlike front tracking, we do not need to follow the individual boundary points from one iteration to the next. This algorithm suggests extra work in everywhere com-puting a function and its velocities for which we desire only a small subset of its values (ie near </> = 0), but we needn't pay too much attention to the accuracy of the solve for 4> far from its zero level set. Furthermore, as previously men-tioned, there may be great advantages higher dimensional applications, where the boundary parametrisation required for explicit front tracking of changing topologies is not trivial. In order to ascertain the stiffness of the hybrid method, as an example, let us return to the linearised system considered in Chapter 2, where a base boundary separates upper and lower regions, D+ and D~, respectively. The Laplacian applies to the dependent variables above and below the free boundary, T. Instead of perturbing the base boundary, we choose a level set function describing the base boundary and perturb that base level set function. We shall reuse much of the notation from Chapter 2. Let 4>(x,y) = <Po(x,y) +e<j>1(x,y), (5.4) where 4>o{x,y) = y, (5.5a) M*,V) = h(t)g(y)eiax. (5.56) We have suppressed the t-dependence of cf> for clarity. Although we shall discover that the perturbation amplitude of the eigenfunctions associated with V do not depend on y, we retain this freedom for the time being. We next determine how the perturbation in (j) affects the boundary location. Assume again that j] ~ 770 + em + (2-3) Now we can compute the boundary location, y = n(x), where (j)(x,n(x)) — 0. Again, we suppress the i-dependence of n for clarity. Computing term by term, <j>o{x,r](x)) = <t>0(x,T]o(x)) + ern(x)<f>0y(x, rjo (x)) + 0{e2), (5.6a) Chapter 5. Continuing Work 49 4>i(x,n(x)) = Mx,m(x)) + 0(e). (5.6b) So, 0 = (f>(x,ri(x)) = 4>0(x,r){x)) + e4>±(x,r)(x)) = 4>o(x,Vo(x)) + 4<f>i(x,ilo(x)) +m(x)<Poy(x,rio(x))] + 0(e2). (5.7) Considering our choices for (po and (pi, and matching terms, we have Vo = 0, (5.8a) 77! . = -h(t)g(0)eiax. (5.86) To within a change in sign, we have the same boundary perturbation we applied in the linear stability analysis of Chapter 2. Other choices for 0O which set y = 0 to as the base boundary are possible (for example, (po = ay, a = const), but this choice makes (po a signed distance function. Thus, the boundary velocity is given by the eigenvalues for v T G K c - h(t)g(0)vTGjlo = g(0)^ = Xh(t)g(0), (5.9) where the sign of h on the LHS has been reversed from that in equation (2.28) according to the change in 771 (a;) by the level set setting. We can repeat exactly the same calculation as that in Chapter 2 to obtain X = - ^ ^ - (5-10) The decay for h(t)g(0) differs only from that of the original linear analysis by a change in sign. Just as in Chapter 2, features of the problem, the linear operator, the global fluxes and the boundary conditions, appear in the numerator, and the affect of velocity choice appears in the denominator. Even in this new setting, A remains dependent on our generalised Stefan velocity choice. We require vn in order to return from this result to the level set evolution equation (5.3). Recalling that "»= % and drj/dt = eg(0)eiaxdh/dt, we have vn = -e\(a)h(t)g(Q)el0cx, so 0 evolves according to 4>t - e\(a)h(t)g(0)eiax\V(p\ = 0 (5.11) near y = 0. Numerically, (5.11) will be solved by a method explicit in time, and must obey a Courant-Friedrichs-Lewy (CFL) condition for stability of its P D E scheme. This condition depends directly on A, contained in the leading term to the spacial derivative, and so we see that level set equation does indeed inherit stiffness properties according to our velocity choice, v. If we substitute cp = cp0 + ecpi into (5.11), and collect the 0(e) terms (the 0(1) terms are zero), we find that g(y) = g(0) is the natural way to extend the Chapter 5. Continuing Work 50 velocity away from T for the linearised problem. However, this analysis assumes that the amplitude, h(t), decays uniformly throughout the domain. If, instead, we replace h(t)g(y) with the single amplitude h(t, y), we see that more generally, the stiffness of the level set solve is related to how the velocity is extended from the r to the entire domain. Such freedom may be more convenient numerically, especially in more complicated domains, and we expect the method of velocity extension to have an effect on the problem stiffness. Velocity extension schemes are presented both in the original paper [1], and in the text [23]. It is not known if such schemes will cause the stiffness of the level set equation to deviate significantly from the predictions for the associated velocity choice. For the time being, such study remains future work. 5.2 Extensions to general operators and base boundaries Before leaving the level set framework, we carry the analysis of Chapter 2 one step further, and consider the behaviour of generalised Stefan schemes applied to a larger set of boundaries and (elliptic) operators. Such extension is motivated by the ability of level set functions to describe boundaries that may not be representable as functions y = n(x). Before incorporating the level set method, we return to the basic front track-ing trial method, and consider the stability of perturbations to more general boundaries. We consider the base boundary location xo parametrised by s. (Our boundaries all have co-dimension 1, so dim(s) = dim(x) — 1.) Perturb the free boundary, x = x 0 ( s ) + ex 1 (s) , (5.12) Just as in Chapter 2, we seek perturbation solutions, u ~ u 0 + eui + (5.13) Again, we use Taylor series to approximate conditions at the perturbed bound-ary by conditions at the base boundary, so u(x) = u ( x 0 + e x i ) , = u(xo) + exi • Vu(xo), = u 0 (x 0 ) + e[ui(x 0 ) + x i • V u 0 ( x 0 ) ] + 0(e 2 ) , = u 0(s) + e[u1(s) + x 1 ( s ) - V u 0 ( s ) ] + O(e 2 ). (5.14) The last equality follows because all boundary information can be parametrised by s. Using the matrix representation, G , for the free boundary conditions, we assume that the base solution is satisfied, so Guo(xo) = b. To satisfy the 0(e) condition, we need the normal derivative, u„(x) = u „ ( x 0 + e x i ) , Chapter 5. Continuing Work 51 = u 0 „ ( x 0 ) + e[ui(x 0 ) + x i • V u 0 ( x 0 ) ] n + C(e 2 ) , = V u n ( x 0 ) • n + e [Vui (x 0 ) • h + x i • n A u 0 ( x 0 ) ] + 0{e2), = V u 0 ( s ) • h + e[Vui(s) • h + x1(s) • n Au 0 (s ) ] + C(e 2 ) . (5.15) (5.16) Next, we generalise the elliptic operator. In each region D± sharing the free boundary F, we solve f L u i = 0 in D± \ u i = ipa(s) on r We expect far-field conditions, such as global flux conditions, on the problem to be satisfied by the base solution, so u i —> 0 far from T. As long as L is linear and separableMn the boundary coordinates, s, this reduces the solution to a linear combination of solutions who's values and normal derivatives at the boundary are represented by i i i = •0 a (s)Kc. (5.17) K again gives us the relationship between the solutions near T, given by ipa a n d ipan- b i the case of the Laplacian with base boundary r)(x) = 0, we had K = — \a\ 0 0 \a 1 0 0 1 (2.18) but in general, we will have K ¥>i(<*) fm{a) 1 (5.18) {(fii(a)} represent the ratios of the normal derivatives the solutions Uj a ( s ) to their values at the base boundary. {ifi(a)} depend on the operator L and the base boundary location, and in effect describe the smoothness of the Dirichlet-to-Neumann map at T. Note that in this more general setting, {ipi(a)} may be complex, yielding complex A. However, the stability criterion illustrated by Figure 1.8 still applies. The appropriate perturbation, x i = h{t)g(s)ipa(s)n(s), (5.19) 1 Separability as we saw for the Laplacian near the base boundary rj(x) = 0 is not guaranteed by linearity. The behaviour of general L near T has not been addressed here, and may indeed be important for the nonlinear, singular examples of operators presented in Chapter 1. Chapter 5. Continuing Work 52 is now apparent. Sturm-Liouville theory states that the set {ipa} will be com-plete, so we are again performing a spectral stability analysis. This means that the C(e)-terms of (5.14) and (5.15) become, respectively, m(s) + h(t)g(s)ipa{s)h(s) • Vuo(s), (5.20) and V u i ( s ) • n + h(t)ipa(s)Au0(s). (5.21) We can now compute A as in Chapter 2, obtaining w J M i e 1 where, according to (5.20) and (5.21), we now have Po A u 0 ( s ) V u 0 ( s ) • n(s) (5.22) As long as each <Pi(a) > 0(1) in a, the derivative terms of K give a similar contribution to the stiffness as in the ice-water example. In Chapter 2, we applied our linear analysis to a variation on the Riabouchin-sky flow problem addressed by Garabedian. We made a simplifying modification to the differential equation, solving the Laplacian rather than the complete flow equation Uxx + u,rr ur = 0. (1.3a) r We are now in a position to consider the linearised free boundary dynamics of this more general linear equation. Let us assume, as in Chapter 2, that the base boundary is at y = 1, that we are solving in the upper half-plane, and that we can linearise the boundary conditions about this point. Equation (1.3a) is separable. We set u(x,r) = X(x)R(r), and choose a for the eigenvalue. The separated equations are X" + a2X = 0, (5.23a) r2R" -rR' -a2r2R = 0. (5.236) The first equation gives us the familiar orthonormal set of functions parallel to the boundary, {elax}. The second equation gives us the solution behaviour perpendicular to the boundary: solutions to this equation determine <p{a). m = 1 for this problem, so there is only one such Dirichlet-to-Neumann function. The solutions to (5.236) are yl\(\a\y) and yKi(\a\y) for y > 0, where h(£) and -Ki(£) are modified Bessel functions of order 1. We use the second solution, which decays for y —» oo. Using the series expansion for Xi (£ ) , cf. [4, ch 11], we compute the asymptotic expansion for V = f e £ « f * f (5.24, Chapter 5. Continuing Work 53 near y = 1, finding that tp(ct) ~ —|a|, just as for the case when the operator L is the Laplacian. Therefore, K ~ [—|a|,l] T . This explains why our Chapter 2 analysis of the problem using the Laplacian, so K = [—|a|,l] T , so accurately explains the success of Garabedian's method. Results regarding the regularity lost in the Dirichlet-to-Neumann map [29], as well as this example, suggest that in general, when L is a second-order elliptic differential operator, tpi(a) ~ |a | . However, the analysis needed to prove this conjecture is beyond the scope of this work. To apply this analysis to the level set method of Osher and Sethian, instead of perturbing the boundary, x(s), we again perturb the level set function, seeking how 0 = 0O + eh = 0 (5.25) affects the boundary perturbation, x i , above. Again, we expand 0 = 0 ( x o + e x i + 0(e 2 )) , = 0(x o ) + exi • V0(x o ) + 0(e 2 ) , = 0 o (x o ) + e[0i(x o) + x i • V0 o (x o ) ] + 0(e2), = 0 o(s) + e[0 1 (s )+xi(s ) -V0 o (s ) ] + (9(e2). (5.26) where we have again suppressed the i-dependence of 0 for clarity. 0o is a base solution, so 0o (s) = 0, and, to first order, 0i(s) = -x i ( s ) -V0o(s ) . (5.27) To isolate x i , we note that we wish to find x i normal to the base free boundary, so that x i = x\fi. Recalling that n = V0o/|V0o|, 0i(s) = - X I T ^ T • V 0 O = -x i |V0o(s ) | , (5.28) |V0o| so as a result of our perturbation to 0o, *' = - » V W s » - < 5-2 9 > Having identified how a perturbation in 0 affects x i , we use equations (5.14) and (5.15) to see that the (9(e) value and normal derivative conditions at xo become, respectively, V0o • V u 0 UI (s) - 0i (s) — — — , (5.30) V0o • V0o and V m t s J . n - ^ f s ) - ^ . (5.31) |V0o| Computing the condition on the normal derivative is, admittedly, some work, and we omit the details. We again write 0i(x) = /i(£)<j(r)^ Q(x), where yba is from a set of orthonormal functions, this time taken over the entire domain. Chapter 5. Continuing Work 54 Wi th r defined to be the coordinates orthogonal to s, g(r) will again represent how the velocity varies away from T. We can again identify the problem stiffness as in equation (2.29), with the modification, Mo = A u 0 V 0 u 0 V 0 O • V 0 O Applying this to the development in the previous section, where <po t)X — h(t)g(0)eiax, we use equation (5.19) to compute X T , X l = -h(t)g(0)eiaxh(s). (5.32) y and (5.33) In this case, n(s) = n(x,0) = (0,1), so this exactly represents the perturbation y = er)i(x) = -eh(t,y)eiax. V0 o (x o ) = V<j>o{x,0) = 1, so since u 0 = u0(y), we have / i = — [uo^y, u o y ] T . This gives the same A as that in the linear analysis, with the change in sign computed in the previous section. Just as with the argument in Chapter 2, to complete this discussion, we would have to verify the existence of base solutions, not only for u 0 , but also for 0o, the corresponding level set function. The conclusion, however, is the same, that there exist generalised Stefan velocities that make the solve i) workable, and ii) 0(1) in the mode number a. This reiterates the idea of the previous section, that regardless of the topology at hand, the level set method of Osher and Sethian can be used in combination with generalised Stefan velocities to form a numerical method for steady state free boundary problems with good stability properties. 5.3 Boundary contact points In our exposition, we have chosen to make predictions and run simulations on periodic domains, such that the free boundary never contacts the outer bound-ary of the problem. This is by no means physically realistic, and near points of boundary contact, we must supply some additional information to the problem to determine the locations and angles of contact between the free boundary and the outer boundary of the domain. In his discussion of the dam seepage prob-lem, Crank mentions the issue of boundary contact, but does not address it in general [8]. Determining contact angles for elliptic solves with our general class of boundary conditions does not appear to have been studied in the literature. We focus here on the well-posedness of free boundary problems near contact points, rather than on numerical schemes. Consider the half-plane of Figure 5.1 divided into two wedge-shaped regions by a free boundary, a straight line through the origin. This is a linearised view near a point of contact, 0, between the free boundary, T, and the physical boundary, dD+ U dD~, of the problem. Here, the point of contact of T is fixed, but the contact angle, <p, in the figure, is not. This suggests that our boundary Chapter 5. Continuing Work 55 8D~ O dD+ Figure 5.1: Free boundary near a point of contact with the physical boundary. data explicitly fixes O, and we expect the behaviour of the solutions u*1 near O to determine <j>. We will use the polar coordinates (r, 9) near O, with the convention at all surfaces that n = 8, so d/dn = (l/r)d/d6. As an example, we use the free boundary conditions for the ice-water prob-lem, for which u+ = u~ = 0 and dnu+ — adnu~ = 0 on T at steady state. Sup-pose we have Dirichlet conditions along the physical boundary. We know that at O, u = 0, so we choose the conditions u+(r, 0) = pr on dD+ and u~(r, n) = —qr on dD~, p, q constants. In polar coordinates, the general solution for A M = 0 is oo u(r,0) = ] T r n (An cos n6 + Bn sin n<9). (5.34) n=0 However, for our linearised treatment, we take only the n = 1 term (the n = 0 term is zero because u = 0 at O) so to solve the ice-water problem in these wedge-regions, we have u ± r(A± cos6 + B±sm.9), (5.35) and we seek the four coefficients A±,B±. The conditions on dD+ and dD~, where 0 = 0 and 0 = TT, respectively, give A + = p and A ~ = <7. At Y, we assume some 0 < 4> < TT, whereupon the conditions u + = 0 and u~ = 0 give B+ = —pcotcf) and B ~ = —qcot<f>. With A ^ B * determined, the final free boundary condition should determine the unknown contact angle, <f>. u\ — au~ = 0 gives —A+ sin <f> + B+ cos 0 + aA~ sin 0 — aB~ cos <A = 0, (5.36) which simplifies to ^ = 0. (5.37) sin <p This is a surprising result, for it does not specify <j> (except to say that sin</> ^ 0), but, rather, specifies the condition p = aq on the boundary data along Chapter 5. Continuing Work 56 dD+ L)dD~, such that all three conditions along T, wherever it may be, can be satisfied. Strictly speaking, even with p = aq, the problem remains ill-posed, for the contact angle of the free boundary remains unknown. As a solution to this difficulty, we suggest applying the mixed boundary data u+-u+=pr ondD+, ,g 3 g . u~+au~ = —qr ondD~. These are radiation, or negative feedback, conditions, the same as those con-sidered on T in Chapter 2. We use the same ansatz (5.35) for the solutions for u1*1 in terms of A±,B±, and applying the revised boundary conditions at dD+ ,dD~, we have A+r-B+^=pr, ( 5 3 g ) —A r — aB = —qr, while the conditions at V again give r{A+ cos 0 + B+ sin 0) = 0, r(A~ cos 0 + B~ sin 0) = 0, (5.40) —A+ sin 0 + B+ cos 0 + aA~ sin 0 — aB~ cos 0 = 0. Comparing powers of r in (5.39) gives A+ = p,B+ = 0,A~ = q,B~ = 0 and with conditions (5.40), the problem seems overdetermined. However, we see that the third condition on T recovers the requirement p = aq, and the first two conditions agree that cos 0 = 0, and thus determine the contact angle, 0 = n/2. Choosing radiation conditions at the outer boundary seems physically rea-sonable, but also turns out to be necessary to give a well-posed problem. How-ever, it is not clear that these are the only conditions on c l D * giving well-posed problems. A more general treatment of this topic with different conditions on T is yet to be done, and remains outside the scope of this work. Yet, this initial investigation indicates that boundary conditions near free boundary contact points cannot be chosen arbitrarily, and a system for setting such boundary conditions must be considered in order to develop algorithms for solving more general steady state free boundary problems. 5.4 Summary This work considers the convergence of free boundary solves by the trial method. In a system of m dependent variables sharing a free boundary, the trial method assumes a fixed approximation to the boundary, solves the problem using m of the boundary conditions, and uses the residual in the remaining boundary condition to adjust the free boundary for the succeeding fixed-boundary solve. This procedure iterates until the residual is sufficiently small. For its motivation from a change-of-phase model, the residual condition, the choice of which is somewhat arbitrary, is called a generalised Stefan velocity. As solving the fixed-boundary problem is expensive, we explore a simple strategy for maximising the convergence rate while maintaining stability. Chapter 5. Continuing Work 57 We consider the general class of free boundary problems, having second-order elliptic operators in each domain, with mixed conditions at the free boundary. As a model for this class of problems, we study systems having boundary condi-tions linear in Dirichlet and Neumann data at the free boundary, and satisfying Au = 0 for each of m dependent variables sharing the free boundary. A linear analysis of such systems reveals that given certain restrictions on global conditions of the problem, as well as local well-posedness, maximising convergence of the trial method without forsaking stability is possible. In choos-ing the residual condition, or velocity, using a boundary condition containing Dirichlet terms yields a numerical method that has a convergence rate, hence stability properties, independent of the discretisation used for the free boundary. This is our main result. This result is confirmed in a computer experiment for a model of a water-ice change-of-phase system, which computes the boundary between the two phases. The test program uses a commercial finite element code to solve the fixed-boundary problem at each iteration. Hence, little original code is required to implement the generalised Stefan method. We close this study by considering a mild extension to elliptic systems with more general operators. We also comment on a hybrid scheme that uses gen-eralised velocities to capture the free boundary using a level set method. This future work can allow the computation of free boundaries with unknown topolo-gies, while retaining the good convergence characteristics of our generalised Ste-fan velocity choices. Bibliography 58 Bibliography [1] D. Adalsteinsson and J . A . Sethian. The fast construction of extension velocities in level set methods. Journal of Computational Physics, 148:2-22, 1999. [2] Robert A . Adams. Sobolev Spaces. Academic Press, New York, 1975. [3] H . W . Al t and L . A . Caffarelli. Existence and regularity for a minimum problem with free boundary. Journal fur die Reine und angewandte Math-ematik, 325:105-144, 1981. [4] George B . Arfken and Hans J . Weber. Mathematical Methods for Physicists. Harcourt/Academic Press, San Diego, 5th edition, 2001. [5] Wil l iam E . Boyce and Richard C. DiPrima. Elementary Differential Equa-tions and Boundary Value Problems. Wiley & Sons, Toronto, 5th edition, 1992. [6] L . Bridge, R. Bredean, M . J . Ward, and B . R. Wetton. The analysis of a two-phase zone with condensation in a porous medium. Journal of Engi-neering Mathematics, 45:247-268, 2003. [7] Y . C. Chang, T. Y . Hou, B . Merriman, and S. Osher. A level set formula-tion of eulerian interface capturing methods for incompressible fluid flows. Journal of Computational Physics, 124:449-464, 1996. [8] John Crank. Free and Moving Boundary Problems. Clarendon, Oxford, 1984. [9] C. Cuvelier and R. M . S. M . Schulkes. Some numerical methods for the computation of capillary free boundaries governed by the Navier-Stokes equations. SIAM Review, 32:355-423, 1990. [10] John R. Dormand. Numerical Methods for Differential Equations. C R C , New York, 1996. [11] M . Flucher and M . Rumpf. Bernoulli's free-boundary problem, qualitative theory and numerical approximation. Journal fur die Reine und ange-wandte Mathematik, 486:165-204, 1997. [12] Andrew Fowler. Mathematical Models in the Applied Sciences. Cambridge University Press, Cambridge, 1997. Bibliography 59 I. A . Frigaard and C. Nouar. Nonlinear stability of poiseuille flow of a Bingham fluid: theoretical results and comparison with phenomenological criteria. Journal of Non-Newtonian Fluid Mechanics, 100:127-149, 2001. P. R. Garabedian. The mathematical theory of three-dimensional cavities and jets. Bulletin of the American Mathematical Society, 62:219-235,1956. B . Gustafson, H.-O. Kreiss, and J . Oliger. Time Dependent Problems and Difference Methods. Wiley & Sons, New York, 1995. Richard W . Hamming. Numerical Methods for Scientists and Engineers. Mc-Graw Hi l l , New York, 1973. Comsol Inc. h t tp : / /www . comso l . com. Accessed June 2003. Ka r i T. Karkkainen and Timo Tiihonen. Free surfaces: shape sensitiv-ity analysis and numerical methods. International Journal for Numerical Methods in Engineering, 44:1079-1098, 1999. J . Ockendon, S. Howison, A . Lacey, and A . Movchan. Applied Partial Differential Equations. Oxford, New York, 1999. J . T . Oden and J . N . Reddy. An Introduction to the Mathematical Theory of Finite Elements. Wiley & Sons, Toronto, 1976. S. Osher and J . A . Sethian. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations. Journal of Com-putational Physics, 100:209-228, 1988. W . H . Press, S. A . Teukolsky, W . T. Vetterling, and B . P. Flannery. Nu-merical Recipes in C. Cambridge, New York, 2nd edition, 1992. J . A . Sethian. Level Set Methods and Fast Marching Methods. Cambridge, New York, 2nd edition, 1999. J . A . Sethian and J . Strain. Crystal growth and dendritic solidification. Journal of Computational Physics, 98:231-253, 1992. N . Shamsundar. Comparison of numerical methods for diffusion prob-lems with moving boundaries. In Moving Boundary Problems (Gatlinburg, 1977), pages 165-185, New York, 1978. Academic Press. J . Simon. Differentiation with respect to the domain in boundary value problems. Numerical Functional Analysis and Optimization, 2:649-687, 1980. J . Sokolowski and J . P. Zolensio. Introduction to Shape Optimization, Shape Sensitivity Analysis. Springer-Verlag, New York, 1992. Gilbert Strang. Linear Algebra and Its Applications. Harcourt, Brace, Jovanovich, San Diego, 3rd edition, 1988. Bibliography 60 [29] John Sylvester and Gunther Uhlmann. The Dirichlet to Neumann map and applications. In Inverse Problems in Partial Differential Equations, (Artaca, CA, Jul 28-Aug 4, 1989), Philadelphia, 1990. S I A M . [30] Roger Temam. Navier-Stokes Equations. North-Holland, Amsterdam, 1977. [31] Timo Tiihonen. Shape optimization and trial methods for free bound-ary problems. RAIRO, Mathematical Modelling and Numerical Analysis, 31:805-825, 1997. [32] Timo Tiihonen. Fixed point methods for internal free boundary problems. Numerical Functional Analysis and Optimization, 19:399-413, 1998. [33] Timo Tiihonen and Jari Jarvinen. On fixed point (trial) methods for free boundary problems. In Free Boundary Problems in Continuum Mechanics (Novosibirsk, 1991), International Series of Numerical Mathematics, 106, pages 339-350, Basel, 1992. Birkhauser. [34] K . S. Udell. Heat transfer in porous media heated from above with evapora-tion, condensation, and capillary effects. Journal of Heat Transfer, 105:485-492, August 1983. Appendix A. Calculation Details 61 Appendix A Calculation Details A . l El l ipt ic i ty of the porous media flow problem We consider the derivation of the problem of two-phase flow in a porous medium introduced in the examples in Chapter 1. By studying the physics of the prob-lem, our aim is to provide enough detail to indicate that the problem is elliptic. In the vapour only zone, we can describe the flow using the two variables pressure Pv(x,t) and temperature T(x,t). The molar density C will be deter-mined by the ideal gas law c = i < A 1 ' where R is the ideal gas constant. Conservation of heat and mass are given below: A T = 0 (A.2a) V - ( C u ) = 0, (A.26) where we assume here in (A.2a) the heat is transported primarily through con-ductance in the solid phase. The average vapour velocity u is given by Darcy's law: n=-Kv{VPv + (0,pg)) (A.3) where Kv is a constant containing vapour viscosity and the permeability of the porous medium, g is gravitational acceleration, and p is the mass density. p = MC, where M = 18 is the molar mass of water vapour. Combining (A.26, A.3 , A . l ) we obtain *•* W o ^ H - a (A.4) RT V V RT Since the highest order term above is K„P, RT APV (A.5) we can recognise (A.4, A.2a) as a (nonlinear) elliptic system for Pv and T. In the two-phase zone we make a saturation assumption, that is that the vapour pressure Pv is a given function Psat{T) that can be experimentally de-termined. Fitted Psat approximations can be found in [6] for example. The Appendix A. Calculation Details 62 saturation pressure is an increasing function. In the two phase zone, the vari-ables of interest are the temperature T(x,y) and the saturation s(x,y), the fraction of pore space occupied by liquid water. Gradients in vapour pressure drive a flow of vapour and gradients in the liquid pressure drive a flow of liquid: ut = - s 3 ^ ( V P ; + (0,p / 9 )) (A.6) uv = -(l-s)3Kv(VPv + (0,pvg)) (A.7) where ui is the liquid velocity, Ki is a constant, Pi is the liquid pressure, pi is the liquid density (taken constant), uv is the vapour velocity, Pv = Psat (T) is the vapour pressure, and pv is the vapour density that can be determined from T and Psat[T) as above. The terms s 3 in (A.6) and (1 - s ) 3 in (A.7) are empirical relative permeability terms: it is easier to move liquid water and harder to move vapour when s is large. The liquid pressure and the vapour pressure are related as follows: Pi=Pv- Pc (A.8) where Pc(s) is the capillary pressure, taken to be a known (decreasing) function of s. We can write the conservation of liquid, vapour and thermal energy as follows: V(piui) = 7 (A.9a) V(p„u„) = - 7 (A.96) A T = -Krl (A.10) where 7 is the local condensation rate and KT is a constant, the ratio of the heat of condensation to the conductivity of the sand pack. Equation (A.10) can be seen to be elliptic in T using (A.96) for 7 by following back through (A.7) and recalling that Pv = Psat(T) is increasing with T. Equations (A.9a, A.96) can be combined to eliminate 7 giving V(piut+ pvuv) = 0. ( A . l l ) Considering ui from (A.6) using the capillary relationship for Pi (A.8) with P'c{s) < 0 it is possible to interpret ( A . l l ) as an elliptic equation for s, although one that is degenerate when s = 0. Thus we have an elliptic system for T and Pv above the free surface T and an elliptic system for T and s below the free surface. We will show below that at the interface, five mixed Dirichlet-Neumann conditions are specified, thus that this problem also fits into our general framework. On the free boundary T we have the following conditions: s = 0 (A.12o) [T] = 0 (A.126) [P„] = 0 (A.12c) (pvuv)+-n = (pvuv + pm)-. • h " (A.12d) ~8T~ (jr = KT{piu{)- (A.12e) on Appendix A. Calculation Details 63 where subscripts + denote limiting values from the vapour region and — from the two-phase region, and [•] denotes the difference in quantities from one side of T to the other. Here, (A.126) and (A.12c) denote continuity of temperature and pressure, respectively, (A.12d) is the conservation of total water, and (A.12e) is a thermal balance, allowing the presence of a condensation or evaporation front at r. Note that condition (A.12a) forces the system in the two-phase zone to be degenerate at the free surface. At physical boundaries we specify two conditions: a given temperature and a no-flux mass condition. Finally, the total water content in the closed system is specified to make this a well-posed problem. A.2 Calculation of A Recall from page 25 that we wish to isolate A in the expression for the velocity, v T G K c + / i v T G / x 0 = ^ = \h. (2.28) dt Define A to be the m+l-by-m matrix with unit vector columns independent of v. A can be interpreted as the linear combinations of interface conditions used in the P D E solve at each iteration. Thus, A T G K c + hATG}lo = 0. (A.13) From A and v construct the matrix M = AT (2.30) Wi th this construction, we can combine equations (2.28) and (A.13) to obtain M G K c + / i M G f i n = [A/i ,0 ,0] T = Xheu (A.14) where e i = [1 ,0 , . . . , 0 ] T . M is full-rank, so its inverse exists, and G K c + hGjlo = A / i M - ' e i . (A.15) Finally, we multiply by the left nullspace of G K = w T , so / i w T G £ 0 = A / i w T v . (A.16) h ^ 0, so we divide through by h and rearrange to give w' M 1 e 1 Expression (2.31) can be simplified further if we can assume something about A , the method used to solve the fixed-boundary P D E . If the columns of A are Appendix A. Calculation Details 64 orthonormal to v , we can row-reduce M so that it is an orthonormal matrix, Q. That is, Q = E M , where E is a row-echelon matrix. The first row of Q is v T and row reduction only affects M by making its rows A T orthonormal to each other. Therefore, " 1 0 ••• 0 " 0 E = 0 where X depends on how A is row-reduced, so M - 1 e ! = Q T E e ! = v. This gives A = (A.17) (A.18) (2.29) the decay rate representing the stiffness of the problem presented in Chapter 2. A . 3 The Fourier integral of the cubic spline Recall from page 33 that we wish to calculate V(n) = ^ fj y(x)e~inxdx, (3.7) where we recall that _ 1 1 66 y { j l ) ~ 2n~ (in)3 (Ax)3 N-1 r X j + l (w% - 1 ) £ w% \ 3=0 Jxi e~mxdx, (A.19) since y"'(x) is constant on (XJ,XJ+I). We compute the integral, as well as substituting Ax = 2n/N and Xj = j2ir/N: y(n) = 1 I 66 2TT \inj (Ax)3 N-1 (W% - 1) w i=o N -in 3b f N 87m \™J TC-1)SWwe"in^(e"in^-l), n j . (A.20) 36 f N_ 87m V'"""' 3=0 N-1 we recall that 6 = TC-lKW^-l)^^'^ 3=0 cos(2im/N) - 1 cos(2im/N) + 2 ' (3.6) Appendix A. Calculation Details 65 and compute TO - iXWy - 1) = 2 ( l - cos . (A.21) Since W^WZ-™3 = 1, the sum is N, and we can simplify the trigonometric expressions using 1 — cos(2#) = 2sin 2(0) to obtain , 3 sin4(7rn/AQ y ( n j ~ cos(27rn/iV) + 2 (TITI /TV) 4 ' 1 ' } This is the Fourier transform of a single discretised mode, W^j*. If, more generally, we have J V - l 3=0 where [yj] is arbitrary data, then the Fourier transform of the cubic spline fit to these points is the sum of the spline fits of each of the modes y(n), so for general [yj], -i ^ 3 sin4(7rn/A0 cos(2irn/N) + 2 (irn/Np as stated in the main text. A.4 The Fourier transform of periodic functions In Chapter 3, page 35, we discuss smoothing by a 2/T-periodic Gaussian window, and as such use the Fourier transform of this window. Let us compute the Fourier transform of more general functions, those of the form oo hn(x)= f(x-2nr). (A.23) r= — oo Our Gaussian window is, of course, a special case of this form. If f(x-2irr) G L 2 [0 , 2ir] for all r 6 Z , and the sum (A.23) converges uniformly for x £ [0,2ir], then ^ { n ) = h fj ( £ f { x ~27Tr)) e~inxdx> -I oo - 2 J T = — J f(x-2Trr)e-inxdx. (A.24) Appendix A. Calculation Details 66 Substitute u = x — 2irr, and recognise that e i n 2 7 T T = 1 for r, n £ Z so f M = IT E / f(u)e-inuein2"rdu, 2 7 7 r=^oo ^ ( 1 - r ) = ^ f(u)e-inudu, = f(n), (A.25) the Fourier transform of f(x) over all space. However, unlike the Fourier trans-form of f(x), x € ( — 0 0 , 0 0 ) , / is to be evaluated only for integer arguments. The Fourier transform of our 27T-periodic Gaussian window follows directly. A.5 Series solution to model ice-water problem The model ice-water problem computed numerically in the examples can be solved analytically by scaling u~(x,y) by a, the ratio of thermal conductivities. We set v(x, y) — u+(x, y) in D+ and v(x, y) = au~(x, y) in D~. This guarantees the continuity of the flux — Vv everywhere in the domain, and determining the location of T amounts to solving for the level set v(x,y) = 0. W i t h this scaling, we have the problem Av = 0 in D v(x,l) =u+{x,l) o n y = l , (A.26) v(x, — 1) = au~(x, — 1) on y = — 1 where D = D+ U D~ is our domain — TT < x < TT, - 1 < y < 1, 27r-periodic in x, and we recall that for our specific problem, a = 10, u~(x, —1) = —1, and ( ^ — - < V < — -+(^ ) = { I e l s 2 e-hL-2 • (4-1) We can solve (A.26) by separation of variables, or by inspection, to obtain the two solutions exp(±ny) exp(inx),n € Z for the homogeneous problem, and ay + b,a,b = const for the non-homogeneous problem. Thus, v(x, y) is given by the series solution 0 0 v(x, y) = ay + b + ^ {aneny + bne-ny)einx. (A.27) n= — 0 0 Solving the problem for v(x,y), then, amounts to determining the coefficients a,b,{an},{bn}. We solve the zeroth-order problem (that for a, b) using the boundary data VQ = —10 at y = — 1 and the minimum of v\ at y = 1, v\ = 1/2. Intuition alone suggests this will work a priori, but the v(x,y) we calculate does turn out Appendix A. Calculation Details 67 to solve the problem. Wi th this boundary data, the inhomogeneous solution is a = 11/2 ,6= - 9 / 2 . Looking now at the homogeneous terms, we note that the factor (aneny + bne~ny) in the sum is the Fourier transform of v(x,y) with respect to x at location y. Thus, comparing this factor to the Fourier transform of the ho-mogeneous data at y = ± 1 will recover {an},{bn}. The homogeneous data at y = —1 is 0, as is its Fourier transform. The homogeneous data at y = 1 is 1 for —TT/2 < x < TT/2 and 0 elsewhere, so its Fourier transform is -!- f ' 2 e-inxdx = — sin ( . (A.28) 2TT J_n/2 nn V 2 I So, matching data at y = ± 1 , ane-n + bnen = 0, (A.29a) anen + bne-n = — s i n ( ^ ) . (A.296) 7rn \ 2 / Solving, we have - s i n ( ^ ) , (A.30a) 17T V 2 / 2 sinh(n) mre~2n : 2 sinh(n) n?r """ V 2 6„ = - _  6 2. n, s i n f g g ) , (A.306) and a„e J _ s , n / m r x s i n h K , + l)] riTT V 2 / sinhf2n) v ; nn V 2 / sinh(2n) Note that this factor is even in n, as should be expected since the boundary data is even. We can use this property to represent the solution as a sum of cosines, . , 11 9 ^ 2 . (rms sinh[n(j/+1)1 . . , . n r , s n=l Observe that this solution does indeed satisfy the boundary conditions, and our choice for the zeroth-order boundary conditions is reasonable. Using this solution, and the derivative vy(x, y), we can use Newton's method to numerically solve for the u-location of the boundary at each fixed x by finding the root of v(x,y) = 0. This solution is an infinite series, yet we can approximate it numerically using.only a finite number of terms. Fortunately, the free boundary location is in y < 1, where the size of the summand decays as 0(e~n/n). When comparing this exact solution to that obtained by a boundary discretised by N points, N terms of this series are more than enough to approximate the exact solution such that the error of the numerical trial method solve far exceeds that of our truncated series approximation of the exact solution. Moreover, the Newton's method solve would generally be fraught with instability due to Appendix A. Calculation Details 68 poor initial root approximations. However, it is straightforward to implement here, as we know the approximate boundary shape from the numerical interface solves, from which Newton's method applied to our exact solution converges rapidly and stably. 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items