Industrial Mathematics MATH 680
Popular in Course
Popular in Mathematics (M)
This 132 page Class Notes was uploaded by Michale Kuhlman on Monday September 28, 2015. The Class Notes belongs to MATH 680 at George Mason University taught by Daniel Anderson in Fall. Since its upload, it has received 28 views. For similar materials see /class/215002/math-680-george-mason-university in Mathematics (M) at George Mason University.
Reviews for Industrial Mathematics
Report this Material
What is Karma?
Karma is the currency of StudySoup.
Date Created: 09/28/15
Notes for Industrial Mathematics Math 680 Contents DM Anderson Department of Mathematical Sciences George Mason University Fairfax VA 22030 USA August 2 1 2008 1 Introduction to Asymptotic and Perturbation Methods 11 Convergent vs Asymptotic Series Landau Order Symbols Asymptotic Equivalence 12 13 14 Asymptotic Sequence 15 Asymptotic Series Expansion 16 Uniform Expansions D Perturbation Methods for Algebraic Equations 3 Nondimensionalization 4 Method of Matched Asymptotic Expansions 41 Opening Example 411 412 413 414 42 MMAE General Discussion 421 422 423 424 425 426 Outer Solution lnner Solution Matching Uniform Solution Outer Solution Stretching Transformations lnner Expansion Overlap lnterval Matching Conditions Composite7 or Uniform7 Expansion 43 Where is the Boundary Layer 1 25 26 27 27 August 21 2008 Industrial Math Notes DM Anderson 2 5 Multi scale Methods 33 51 Two time Method 36 52 Three time Method 38 53 Nonlinear Problems 42 54 Exercises 47 6 Basic PDE Theory 49 7 Green7s Functions Methods for ODEs 51 71 Comments on Distribution Theory 52 711 Distributions and Differential Equations 53 72 Introductory Example 55 8 Heat Transfer 59 81 Introductory Concepts 59 82 Conservation Law for Heat 59 83 Boundary Conditions 60 84 Uniqueness of the Solution 61 85 Solution Methods 62 851 Similarity Solution 1D Heat Transfer From Heated Boundary 62 852 Similarity Solution 1D Heat Transfer Localized Heat Source 64 853 Green7s Function solution for the 1D Heat Equation 66 86 Homogenization Methods 71 9 Solidi cation of a Pure Material 74 91 Basic Concepts and De nitions 75 911 Interfacial Conditions 76 92 1D Freezing From a Cold Boundary 77 93 1D Freezing Into an Undercooled Melt 80 931 Attachment Kinetics Effects 81 94 Morphological Instability 83 941 Curved Interfaces 84 95 Morphological Instability of a Planar Interface Growing Into an Undercooled Melt 86 951 Base State Solution 87 952 Linear Stability Analysis 88 96 Homogenization of the Stefan Problem 92 10Fluid Dynamics 96 101 Lagrangian and EuIerian Descriptions 96 102 Conservation of Mass 98 August 21 2008 Industrial Math Notes DM Anderson 3 103 The Stress Tensor 99 1031 Newtonian Fluid 99 1032 Non Newtonian Fluids 100 104 Momentum Equation 100 105 Vorticity 102 106 Boundary Conditions 103 107 Recap Navier Stokes Equations 105 108 Nondimensionalization 105 109 Some Simple Flows 106 1091 2D Parallel Flows 107 1010Approximate Solutions to NS Equations 109 10101 General Lubrication Theory 110 1011Axisymmetric Gravity Current 113 1012Unsteady Parallel Flow 120 101210scillatory Flat Plate 120 10122lmpulsively Started Flat Plate 124 Appendix A 131 Appendix B 131 Abstract These notes are based on class notes from a number of different courses taken at Northwestern University D20 7 asymptotics D11 7 Green s functions Solidi cation Theory Fluid Mechanics D26D27 1 Introduction to Asymptotic and Perturbation Meth ods Good references for this material are Hinch 11715 21 22 23 Bender and Orszag p 89E 118E and Fowler 31733 11 Convergent vs Asymptotic Series see Bender and Orszag p 118 Suppose we have the function f and the series 2 DEFINITION 11 The series converges to f at a red z x0 if there erists ari iriteger N such that 7 221 lt E for h gt N DEFINITION L2 The series is asymptotic t0 f as x 7 0 iffor red N arid 6 we have lf 7 251 lt E as x 7 0 August 21 2008 Industrial Math Notes DM Anderson 4 Another way to state this is Convergent ESQ fnz a 0 as N a 00 for xed x Asymptotic ESQ fnx ltlt fNz as x a 0 for xed N Note that not all asymptotic series are convergent Convergent series is an absolute concept in that we can determine whether or not a series is convergent without knowing to what it converges However an asymptotic series must be asymptotic to something EXAMPLE 11 The following example shows important di erences between asymptotic and convergent or divergentl series See also Fowler p 37738 De ne the function Fz zeZZOO gdt 11 We are interested in the behavior of for large 2 Integration by parts gives Z 67 00 00 67 Z 00 67 7 2e 772 71 th 1ize A tTdt 12 Integration by parts again gives Z equot 00 00 2equot 1 Z 00 2equot 1ize 7721 71 t3 dt17zez t3 dt 13 Repeating this gives 1 2 3 71 1n 71 n z 00 e t 1772n 7171n12eZ Wdi 14 Remainder Note that 71 n1 132 7g 2 7 15 n0 2 since this series diverges for all 2 since eventually n gt z for any 2 However N 71 711 in is asymptotic to as 2 a 00 for ped N 16 n0 2 That is the remainder goes to zero faster than the last term kept For epample suppose z 10 Summing the rst 9 terms gives 091582 and summing the rst 10 terms gives 091546 suggesting that 9 terms gives accuracy to the third decimal place August 21 2008 Industrial Math Notes DM Anderson 5 It turns out that does have a convergent power series representation It is 2 3 77 17 2 F z i 71 i 2 2e 7 nzz 22 33 where the Euler Masheroni constant 7 06772156649 eg see Abramowitz and Stegunl It is interesting to note that in order to achieve three decimal place accuracy for z 10 using the convergent series we would need to use 40 terms compared to 9from the divergent series The divergent series is an asymptotic series for as 2 a 00 Consider the following example Note that Z 18 is an N term asymptotic expansion for as 2 a 00 if 7 2712 anznl a 0 faster than aN12N 1 a 0 That is7 for an asymptotic series we require that the terms thrown away approach zero faster than the last term kept as 2 a 00 In terms of the above expansion we require that a a Zlltlt aszaoo 19 N z z 1 12 Landau Order Symbols Landau order symbols are used to characterize the size7 of different quantities in particular limits or on particular domains The following are a few de nitions and examples DEFINITION 13 Big O f Ogx for z E D if there epists a constant k such that lfl lt hlgl for allz E D Given a domain D we must be able to nd a k that works for all z in that domain A related de nition is DEFINITION LA f Ogx as x a 0 if there epists a constant k and a neighborhood Nx0 near 0 such that lfl lt hlgl for all z E Nx0 Here we nd a constant k and a neighborhood around do such that the inequality holds for all z in that neighborhood DEFINITION L5 Little o f ogz as x a 0 iffor any 6 gt 0 there epists N0 such that lfl lt 6lgl for all z E Nx0 August 21 2008 Industrial Math Notes DM Anderson 6 This says that no matter how small you Choose 6 you can always nd a neighborhood around do such that the inequality holds DEFINITION 16 fp7t 0gx for z E D if there epists kt such that lfl lt ktlgl for all z E D This de nition is the same as Def 13 except that now the parameter k depends on another parameter or variable in the problem DEFINITION L7 fp7t 0gx as x a do if there epists kt and Np0t such that lfl lt ktlgl for allz E Nx0t Here we note that if h and N are independent of t then f is Big 07 of g uniformly in t Simply put7 f 0g means that fg is bounded in the appropriate limit or domain f 0g means that fg a 0 in the appropriate limit or domain see for example7 Kevorkian and Cole7 p 1 Often in this context we have fx7 6 and think of z xed as e a 0 EXAMPLE L2 The following epamples demonstrate the above de nitions of0 and o o sinx 01 for all real x Here given z E R choose h 2 say so that lsinpl lt 2l1l for all z E R o x2 0z as x a 0 is technically true although perhaps not as informative as possible In this case we can choose h 1 and Np0 0 as a neighborhood with radius less than 1 Then p2 lt 1m for all z E Np0 0 A more informative statement is necct x2 0p as x a 0 Here for any 6 gt 0 we can nd a neighborhood nearz 0 where le lt That is if6 001 then le lt 001M for lt 001 Ifrl 0001 then le lt 0001lxl for lt 0001 etc o x 0x as x a 0 for all or gt 1 This follows from a similar argument applied in the previous case e m 0z as x a 00 for all n Here we have that le wl lt SH ll for any 6 gt 0 Note that we can think of this as e w1p pne m a 0 as x a 00 for all n o tcost 01 for0 lt t lt T Here we can choose h gt T so that ltcostl lt o tcost 31 01 for allt gt 0 Here for any h we choose independent oft we can always ndt such that ltcostl gt August 21 2008 Industrial Math Notes DM Anderson 7 o tcost 0t for all t gt 0 Here choose h 11 for epample tcosx 01 for 0 lt t lt T andfor all x tcosx 0t fort gt 0 and for all x e a6 06 as e a 0 for all n where a is a cced positive constant Again we note here that faEVE a 0 as E a 0 o e w6 06 as e a 0for alln providedx is bounded awayfrom zero eg z 2 f 2 0 eme 31 06 as E a 0 for all n ifp ge0 Note that ifz 0 then klenl lt1 as E a 0 Based on the previous encamples we say e w6 is not uniformly 06 as e a 0 for 0 S x S 1 0 However e m6 is uniformly 06 for 0 lt f S x S 1 o e w6 01 ase a 0for0 S x S 1 Here we can choose h 2 say so le mel lt 211 13 Asymptotic Equivalence We write fx g as x are if fga 1 as x are This reads f is asymptotic to g as x approaches x0 if the ratio of f to g approaches 1 as p approaches x0 EXAMPLE L3 The following are emamples offunctions that are asymptotically equivalent a x14x asxaoo since 14gt1 aszaoo b sinhz ew as x a 00 since sinhz em 7 e w2 a ew2 as x a 00 c f g as x a no if fx g1 01 as x a 00 Here f gets arbitrarily close to g as x a 0 14 Asymptotic Sequence DEFINITION 18 see also Kevorkian and Cole p is an asymptotic sequence as x a no if pn own as x a 0 for all n EXAMPLE LA The following are emamples of asymptotic sequences a is an asymptotic sequence as x a 0 since xn 0x as x a 0 That is pn x a 0 as x a 0 b is an asymptotic sequence as x a 00 c 6 is an asymptotic sequence as e a 0 d is an asymptotic sequence as x a 00 if An gt An for all n eg An 7104 August 21 2008 Industrial Math Notes DM Anderson 8 15 Asymptotic Series Expansion DEFINITION 19 see also Bender and Orszag p 89 or Hinch section 22 p 20 N fx N Z fnwnw as x a 07 n1 is an asymptotic series eccpansion with respect to the asymptotic sequence if N Z 0 F7nl39 n1 That is the stu that gets thrown away is smaller than the last term kept So N Wit 7 1 fnwzl own Consider the asymptotic expansion N 9570 Z ant9 OWN as 95 m 9507 to S t S It17 n1 N Zaath 110 n1 Here the coef cients are given by z t 7 for t z W hm f 7 gt H lt w gt39 111 maze pn For example7 consider the asymptotic expansion xriolwrilwlagltpg7 as pare 112 Here we have i 0 a0 7 lim fz07 113a mama al hm m hm M7 11310 mama p1 14mg p1 a2 hm W7 1130 mama p2 etc All of these limits must exist in order for fp to have an asymptotic expansion as x a 0 Note that in general both the numerator and denominator approach zero so that these are inde nite limits Another way to look at this is to directly solve for the coef cients aj and take the limit z a 0 For example7 aof ailt71azlt 2m7 114 August 21 2008 Industrial Math Notes DM Anderson 9 so that 10 hm lfz7a1ltp17a2ltp27l 115 14mg implies a0 flt0 Similarly solving for al and taking the limit gives a hm 7 116 mam lt71 implies that al hmx a x0fz 7 a0ltp1 These results follow since pgp1 pgp1 a 0 as x a 0 since n is an asymptotic sequence In this way when fx and pnx are given the coef cients 17 are uniquely determined So in the following situation where N f N 20n n67 13917 with 6 ltlt 1 given f6 and ltp6 the coef cients can are unique However consider now the following situation where the coef cients can an6 are also functions of 6 Then in the asymptotic expansion N f6 N glauewe 118 it is not true that the coef cients an6 are unique when given f6 and pn6 See also Hinch p 22 This situation will come up for example in our discussion of multi scale methods OPERATIONS ON ASYMPTOTIC SERIES EXPANSIONS We can multiply divide add subtract and integrate term by term but only under certain conditions can we differentiate see also Bender and Orszag p 1257128 For example given N f96 N Z awnm 119 n1 we can write N f W N Z aw xm 120 n1 only if f exists and is continuous and ii f has an asymptotic series expansion with respect to the asymptotic sequence August 21 2008 Industrial Math Notes DM Anderson 10 16 Uniform Expansions Consider the asymptotic expansion fx7e Zanz6 121 n0 This expansion is uniform in 171 if the expansion is valid for all z E cub That is7 EiN u i a z a 0 as e a 0 12 If this condition holds for z E cub for N independent of x then the asymptotic expansion is uniform For example7 consider ooen f76 7 as 60 123 n0 n This expansion breaks down when x 06 and so the expansion is not uniform for z E 07 1 It is uniform for z 6 71 where f gt 0 for E a 0 2 Perturbation Methods for Algebraic Equations Consider the algebraic equation 67712 l 2m107 21 with 6 ltlt 1 Since this is a quadratic equation we can write down its exact solution mm 7 i 22 If we expand 1 7 E in a Taylor series about 6 0 we nd that 1 1 1 m1 N ii 7 g5 7 E62 7 H Regular root7 233 m 2 116 162 Sin ular root 2 3b N 7 7 7 7 H 39 2 e 2 8 16 g In general for algebraic or other problems we do not have the exact solution available and so we would like to develop a method that will allow us to write down the asymptotic series expansion directly from the original equation In our rst attempt to do so7 we make the following ansatz or assumed form of the solution mm0m16m262 24 Then7 m2 N 7713 2771077116 27710771262 25 August 21 2008 Industrial Math Notes DM Anderson 11 We plug these expansions into equation 21 6 mg 2771077116 2m0m262 2 me mle 771262 1 07 26 and then rearrange by collecting powers of E 27710 1 E 27711 mg 62 2mg 277107711 0 27 Since 17 67 627 are linearly independent each coef cient rnust vanish identically7 so 1 m0 757 28a 1 1 m2 7 28c Therefore7 the assumed form of the solution gives mwiiieeiiez7 29 which is regular root identi ed in the expansion of the exact solution We have not identi ed the singular root with this particular approach since our assurned form of the solution was that of the regular root but not the singular root If we look at the singular root7 it suggests that the correct assurned form of the asymptotic expansion is mEm0mlenn 210 e However7 it would be nice if we could somehow deduce this by some method that did not depend on the availability of the exact solution Suppose we introduce a change of variables 71 m 67 211 which can also be thought of as a stretching77 if E is small and 04 gt 0 Plugging this change of variables into the original equation 21 gives 1720 2 for i E n 26 n 1 70 212 Term 1 Term 2 Term 3 Now7 depending on how we choose 04 we get different balances7 as e a 0 An appropriate balance is actually a dominant balance in which the terms that balance in the limit as E a 0 are also the dorninant7 or largest terrns7 in that limit In the above exarnple7 there are three terms and so there are three possible two terrn balances We examine all cases here and assess each one to identify a dominant balance August 21 2008 Industrial Math Notes DM Anderson 12 Balance Term 1 and Term 2 In this case7 to balance the terms as e 7 0 we require that 1 7 2o 704 which implies Oz 1 This constitutes a balance7 of the rst two terms It is also a dominant balance since the rst two terms are 06 1 ie large as e 7 0 while the third term is 01 This appears to be a suitable scaling to pursue7 however7 we shall examine rst the other two possible balances Balance Term 1 and Term 3 In this case7 to balance the terms as e 7 0 we require that 04 12 In this case7 however7 it is the second term that dominates in the limit 6 7 0 so that this choice does not constitute a dominant balance Balance Term 2 and Term 3 In this case7 to balance the terms as e 7 0 we require 04 0 While this technically leads to a dominant balance since Term 1 in this case is 06 we have already examined this limit in our original attempt at solving the problem So7 in agreement with our expectation from the exact solution we take 04 1 so that the rescaled problem to solve is n22n60 213 We seek a solution to this problem of the form nno6n162n2 214 Plugging this into equation 213 gives ng2non162no6n160 215 At 01 we have 713 2710 07 216 which has solution no 0 and no 72 At 06 we have 2non1 2n1 1 0 217 The solution to this problem depends on the value of no If no 0 then n1 712 If no 72 then n1 12 At 062 we have 2710712 nf 2712 0 218 When no 0 and n1 712 we have no 718 When no 72 and n1 12 we have This rescaling has led to the solution 1 1 n N 0 7 i6 7 ez n 721616 2 8 7 219a 2 21 August 21 2008 Industrial Math Notes DM Anderson 13 In terms of the original variable in ne we have 1 1 220 m 2 86 a 2 1 1 77 7 7 22010 W 6286 These results agree with the expansion of the exact solution in 23 Note we have identi ed not only the singular root in this case but also have recovered the regular root we had derived earlier The reason behind this is that the leading order problem in the latter case was a second degree polynornial leading to two solutions while that in the former case was a rst degree polynornial leading to only one solution In fact7 the singular nature of this problem follows from the loss of the highest power of m in the original equation as e a 0 EXAMPLE 21 This example demonstrates a more general approach to singular algebraic problems see problem 13 in Hinch Find the rst two terms in an asymptotic eccpansion for 6 ltlt 1 for all four roots of the equation 6p4iz27x20 221 We note that without rescaling the limit 6 a 0 gives 7z2 7 z 2 0 or z 72 and z 1 the leading order parts of only two roots are revealed Note that this simple calculation has implicitly assumed that the value ofx is 01 or more speci cally that 6x4 01 as e a 0 We shall nd shortly that the remaining two roots do not have this property and consequently have not been revealed so far We shall restart this problem by rst introducing a rescaling X p 6 222 Plugging this into equation gives 614D X4 i e ZD XZ i X 2 0 223 In order to choose or we seek a dominant balance of terms in the above equation Note that we shall only consider balances that involve the rst term the 4th degree term Balancing the rst and second terms requires 1 7 4o 7204 or or 12 This leads to a dominant balance which we shall pursue Note that balancing term 1 and term 3 leads to or 13 but that term 2 dominates in this case no dominant balance Balancing term 1 and term 4 leads to or 14 but again term 2 dominates so again this choice does not lead to a dominant balance We take or 12 which leads to the rescaled equation X4 i X2 i e X 2e 0 224 August 21 2008 Industrial Math Notes DM Anderson 14 The leading order problem in this rescaled problem is X2X2 7 1 0 This shows the four roots X i100 We already know that the two zeros that appear here must correspond to the two regular roots we identi ed in our rst calculation Therefore our leading order estimates for the four roots are given by am As we have already seen one cannot identify a solution if the assumed form of the solution is not correct In the present case we do not know in advance the appropriate scaling for the correction terms for the above four roots To address this situation we can let the problem tell us what the expansion should be by starting with a more general form for the solution That is we start with x 606x0 616x1 626x2 226 where 606 gtgt 616 gtgt 626 as e a 0 Here 606 is an asymptotic sequence and for regular roots we have 606 1 It is also assumed that x0 x1 x2 etc are 01 as e a 0 We shall make use of this approach for each root individually in this problem making use of the information we already know about the leading order terms in the expansions First Root We seek an expansion of the form where 61 ltlt 1 but is otherwise still to be determined We plug this expansion into equation 2 and retain only the minimum number of terms required to determine 61 and x1 We have 6726114726l l27726112 0 0 mwmn HQ The nal form of this equation tells that we should choose 61 E and x1 7163 Therefore the rst root is approximated by x N 72 7 1663 Second Root We seek an expansion of the form where 61 ltlt 1 In this case we nd l o 6726117611 0 August 21 2008 Industrial Math Notes DM Anderson 15 So we choose 61 E and 1 13 Therefore the second root is approccimated by z N 163 Third and Fourth Roots In this case we work with the rescaled version of the original equation X4 7 X2 7 elZX 26 0 and so take XX061X1 231 Plugging this in and throwing out unnecessary terms giues 07 7 7 261XOX1 7 E XO 7 661X1 26 70 Edamining this equation suggests that we should take 61 6 Then 4X8 7 2X0X1 7 X0 07 233 which if X0 i1 for roots 3 and 4 then X1 12 for both roots Therefore the third and fourth roots are approccimated by X N 1 ei and X N 71 ei or in terms of the original scaling z N i and z N 767 EXAMPLE 22 Another problem to consider is from Hinch Z5ab p 11 Find the rst two terms in asymptotic eccpansions for the three roots of 6953 z 2 i ez 1 07 234 for the two cases indicated by the i signs 3 Nondimensionalization Consider the following model for a mass on a spring with damping Newton7s law says F Ma where F are the forces applied to the mass M and a is the acceleration of the mass Using Hooke7s law and a linear damping term we nd F 7KY 7 BdYdT where K is a spring constant and B is a damping coef cient K has dimensions of mass per unit time squared and B has dimensions of mass per unit time Y measures the displacement of the mass from an equilibrium position and has dimensions of length The governing equations for this system are dZY dY Y0 07 3110 dY IO Wo M7 31c August 21 2008 Industrial Math Notes DM Anderson 16 where IO is an impulse dimensions of mass times length per unit time and can be interpreted as 0 A FdT M V0 i V0 32 These initial conditions correspond to an initially undisplaced mass that is given an initial velocity at time T 0 We make this system dimensionless by introducing a length scale l and a time scale 739 and dimensionless variables y and t Y lly7 T wt 33 Introducing these new variables into equation 31 gives Iggggaggumy 0 34 210 07 34b 0 34c We now ask the question How do we choose l and 739 Below we outline a number of different choices that could be made CHOICE A Suppose that we anticipate that the damping term is going to be a relatively small effect Then an appropriate choice of time scale can be obtained by balancing the rst and third equations in the ODE ie the time scale associated with undamped oscillations Choosing KW 35 gives Tl ltgtl2 36 An appropriate length scale may be the maximum displacement of the spring7 however7 this quantity is not known in advance A suitable choice for the length scale7 then7 may be to make the initial condition for velocity reduce to dydt0 1 This choice gives ll 1 m 37 so that our length scale is m L 38 August 21 2008 Industrial Math Notes DM Anderson 17 Then our governing equations can be written as dzy dy 7 2i 0 39 dtz edt 11 7 21 110 07 39b dy 7 0 1 39 dtlt 7 C where we have de ned a new7 dimensionless7 parameter E e 310 At this stage we have not made any assumptions about the size of the parameter 67 this choice of nondimensionalization may be appropriate for physical systems in which 6 ltlt 1 ie when the time scale for damping is much longer than the time scale associated with undamped oscillations CHOICE B Suppose we are interested in a system in which the mass is small7 in some sense In this case we may nd it appropriate to choose 739 based on a balance of the second and third terms in equation 31 This leads to the time scale 7 7 311 Then if we choose the length scale l based on a balance of terms in the second initial condition we have 10BK l 7 312 H M lt gt Then the governing equations are dzy dy 7 7 0 313 edtz dt 11 7 110 07 314 dy am 7 17 315 where MK 6 E 316 CHOICE C Another choice is B IO M f7 lll W7 317 August 21 2008 Industrial Math Notes DM Anderson which gives dzy dy E dLEdLy 7 110 07 dy 1 where iMK e 313 319 320 321 CHOICE D An additional choice in which the length scale is independent of the mass M B IO m7 mig which gives dzy 7 0 Edit dt y 7 dy 1 where MK 6 E 7 32 This may be useful if one wishes to study the effects of small but changing M 4 Method of Matched Asymptotic Expansions See exercise 51 in Hinch7 pp 527677 Bender and Orszag7 Chapter 97 pp 4197430 41 Opening Example Consider the boundary value problem dzy dy 7 07 110 07 111 1 This problem has exact solution 671 7 eime Wm m 322 323 324 325 326 41a 41b 41c 42 Examining a graph of this solution shows that there is a rapid transition that occurs near z 0 the so called inner region and a region of moderate change the so called outer region This is a singular problem a simple expansion fails7 as we see below August 21 2008 Industrial Math Notes DM Anderson 19 411 Outer Solution Suppose that y 6 y0x 6y1x 62y2z 43 The 01 problem is given by 7 0 44 dx y0 This is a rst order problem that has solution y0x Aoe m 45 There is only one unknown constant to determine in this solution but we have two boundary conditions in the original problem We cannot satisfy both conditions with this solution Based on the exact solution7 we can anticipate that this result represents the solution in the outer region and so we apply the boundary condition at z 1 This gives y0z el w The behavior we see here is a signature of a singular problem we lost the highest deriva tive when we took the limit 6 a 0 Let us continue with the expansion despite this to see how the solution progresses At 06 we nd iyl iiii07 46 subject to y11 0 Therefore y1x 0 It follows that 0 for all n 2 1 Note that y e w is actually an exact solution to the original problem 41 but does not satisfy both boundary conditions One ofthe problems is that we have assumed a series in which each term is independent of 6 Therefore7 there is no way that y0x can describe the region of rapid change in 0 lt z lt 6 The solution we have written down describes the outer region So as an outer expansion we can use y 6 y0x 6y1x 62y2z 47 412 Inner Solution We now need to nd an inner eapansi0n Recall that the exact solution had a term like e we The idea is that we need to rescale or enlarge the region near z 0 To do so we introduce an inner variable X me This is also sometimes called a stretching transformation The rapid variation of the solution near z 0 in terms of the original variable x now looks like moderate variation with respect to the variable X The solution is smoother and in some sense more well behaved in terms of the variable X August 21 2008 Industrial Math Notes DM Anderson 20 One may ask How can this new variable be predicted if the exact solution is not available The idea is to rescale the variable by introducing X 3 48 60 and transform the original equation in terms of the new variable X This gives dzy dy 1720 70 7 1 7 0 49 E dX2EE ch y Term 1 Term 2 Term 3 We seek a dominant balance that involves Term 1 and at least one other term in the equation Balancing Term 1 and Term 2 as e 7 0 implies that 04 1 In this case Terms 1 and 2 are the dominant terms as e 7 0 Balancing Term 1 and Term 3 suggests that 04 12 However7 in this scenario7 Term 2 is actually dominant as e 7 0 and so this balance is no good Consequently7 we choose the former and take X xe With this choice we can identify an inner problem We let y Y denote the solution in the inner region The inner problem is given by dY 7 1 7 Y 0 410 dX2ltegtdXe lt gt subject to the boundary condition Y0 0 Note that we do not apply the boundary condition y1 1 in the inner region since z 1 corresponds to X 7 00 as E 7 0 The inner and outer solutions will ultimately be patched7 together using a matching procedure We solve equation 410 by introducing a standard power series YX Y0XEY1X62Y2X 411 At 01 we nd 07 412a Y00 0 41210 which has solution Y0X B0lt17 5X 413 where B0 is a constant to be determined At 06 we nd d2Y1 le dYO 7 777573 414 dXZ dX dX 0 07 2 Y10 0 414b August 21 2008 Industrial Math Notes DM Anderson 21 which has solution Y1X B1 1 7 5X 7 BOX 415 Summarizing7 we outer and inner expansions youte x 617m transcendentally small terms7 416a 14W 7 3017 5X 4 43117 5X 7 BOX 4 08 416b Note that in the inner expansion we have an undetermined constant at each order in 6 This occurs because we are solving a second order problem for the inner solution but only have one boundary condition to apply there These constants will be determined by matching The idea for matching is that the two expansions hopefully connect together smoothly over some transition region Assuming the existence of an overlap region7 the idea is that these two results are asymptotic expansions of the same function so they must match7 or agree 413 Matching There are a number of principles that can be used to match inner and outer asymptotic expansions We describe one method in this example and mention other approaches later Our approach here is as follows We rst write the outer expansion in terms of the inner variable and collect powers of E that is7 expand in E for xed X 6x 1 yam el m el EX e 17 EX EEZXZ 063X3 417 Note that this is equivalent to expanding the outer solution near z 0 We denote 41231124 e 0z 418a 413192410 e eHX 0z2 418b Similarly7 denote 5635151 301 6 06 41 563431250 3017 Ex e 311 7 e X 7 BOX 08 41 Now7 we require for xed z in the overlap region as e a 07 X ze a 00 that the two asymptotic expansions represent the same function If we rst apply this procedure at 01 we require that 1 0 1 0 1 gljgo ymftgr 7 751 This translates to the condition 130 5 0z 7 B0 06 5 7 0 421 August 21 2008 Industrial Math Notes DM Anderson 22 This condition tells us that B0 e and that these approximations require that z lt 1 from outer solution expansion and e X ltlt 1 inner solution The latter translates to z gt 6 Therefore7 the matching interval for 01 matching is given by 6 lt z lt 1 This con rms that the assumed matching interval does exist See also Hinch7 p 57 and Bender and Orszag7 p 430 We apply the same procedure7 but now to 06 Here we require 41m 441244 7 aggro o lt422gt This translates to the condition lim 6 676X 0z2 7 B0 7 6B1 7 BOX 0627 e X7ee X 0 423 This condition tells us that B0 e and B1 0 note we could have skipped the matching procedure at 01 The matching interval is determined by noting that the above expansions require in the outer expansion that 2 ltlt 67 which is equivalent to z lt 612 and in the inner expansion e X ltlt 67 which is equivalent to m gt el lnel Finally7 we de ne 4351120 E gim yiquot l Xhm Yig yjoklye7eeXOz2 424 These results are used to identify a uniform solution as described below Note that other matching principles7 can be used For example7 Van Dyke7s Match ing Principle is one such procedure see Hinch7 section 5157 p 58 Also methods using intermediate variables are possible 414 Uniform Solution We de ne a uniformly valid solution by yuniformz youte39r Yinne39r 7 ymatclu 617m 617 67X 676X 7 617 6X7 f 1 EH lter inner match 5H 7 el X 425 We can compare this solution with the exact solution given by equation 42 671 7 6715 7 7X yemact 671 7 67176 N 55 m i 5 7 for small 6 In this case7 the boundary layer and matching procedure has missed an expo nentially small term 6 15 see also Bender and Orszag7 p 430 A comparison of the outer solution7 the inner solution7 the uniform solution and the exact solution is shown in Figure 1 This gure indicates the regions for which the inner solution and the outer solution agree well with the exact solution The uniform solution is clearly very accurate These gures correspond to E 01 August 21 2008 Industrial Math Notes DM Anderson Method of Matched Asymptotic Expansions s 01 l l l l l l l O exact solution 7 uniform approx outer solution 7 inner solution Figure 1 This gure compares the exact solution7 the outer solution7 the inner solution and the uniform approximation based on the method of matched asymptotic expansions August 21 2008 Industrial Math Notes DM Anderson 24 42 MMAE General Discussion The following discussion follows from Reiss D20 notes p 41745 Consider the function u76 67 427 This function shows boundary layer behavior near z 0 as e a 0 For xed z 74 0 we have 1355 0 428 while for z 0 1355 1 429 We see that the result depends on the order in which the limits are taken That is7 0 liamjuw7 6 74 lirr iirnJuz7 6 1 430 Therefore u7 6 does not converge uniformly in x as e a 0 For small 6 and small x ux7 6 changes rapidly from u m 0 to u 1 this is boundary layer behavior Another example is e u z E 7 431 lt gt i E lt gt foreZOandzZO One way to resolve this type of nonuniformity is to stretch out77 the region of nonuniform convergence In these examples the nonuniformity is near z 0 but this is not generally the case Here is how the method of matched asymptotic expansions works Suppose we seek a solution ux7 6 to an ordinary differential equation7 for example 421 Outer Solution Suppose ux7 6 has an asymptotic expansion u7 6 u0x u16 u2x62 432 where E a 0 and z 01 as e a 0 In general this expansion is not uniformly valid as E a 0 the expansion is singular at one or more values of For simplicity7 suppose the singularity is at z 0 So the outer expansion is an asymptotic expansion ofuz7 6 as e a 0 for z bounded away from x 0 An example of such a function and expansion is up x6 lt17 mgt5 i063 433 This expansion is valid only away from x 0 August 21 2008 Industrial Math Notes DM Anderson 25 422 Stretching Transformations We stretch out the region of nonuniformity by introducing an inner variable or boundary layer variable r 5 456 434 where gt6 is a function to be determined7 by nding an appropriate balance of terms7 but for now has the properties that gte a 0 as E a 0 and 0 0 Examples include xe7 xelZ7 etc So for xed 7 5 a 00 as E a 0 The region of rapid variation gets smoothed out If we follow our previous example7 E 1 u76 6 1 435 If we let xe7 we write ue 7 6 E U 7 6 1 1 This last form clearly converges as E a 0 uniformly in g to 1 1 In the case where we are solving an ODE rather than looking at an explicitly de ned function we use the stretching transformation to nd a rescaled ODE 4 2 3 Inner Expansion Here we seek an asymptotic expansion of U 767 ie the inner solution7 with xed g as e a 07 U e N 105 6U1g 436 This series is assumed to be uniformly valid in g for the interval 0 3 00 4 2 4 Overlap Interval A fundamental assumption of the MMAE is that there exists an overlap region where both the outer and inner expansions are valid asymptotic expansions ofu7 6 That is7 the outer expansion is assumed to be valid for z gt x06 and the inner expansion is assumed to be valid for z lt x16 and we assume i x06 S x16 as e a 07 ii x06 and x1e a 0 as E a 07 iii z06 6 and z16 e a 00 as E a 0 The rst point says that an overlap interval exists this not always true and some cases require one or more intermediate expansions in addition to the inner and outer solutions The second point tells us that the overlap window shrinks and moves to the singular point as E a 0 The third point tells us that for values of z in the overlap region7 x e a 00 as E a 0 The region of nonuniform convergence is stretched out to 00 August 217 2008 Industrial Math Notes DM Anderson 26 425 Matching Conditions Based on the assumptions given above7 since the outer and inner expansions are asymptotic expansions of the same function u7 6 they must match There are at least two ways to accomplish this STANDARD METHOD In the rst method we write the outer expansion in terms of the inner variable ab Wt 6 uo 6U1 437 We then re expand this expression in E for xed g as e a 0 This gives u 7400 uglt0 lte 7 67410 4 1105 438 So 2 496 7 U496 0 439 0 HF Outer in Terms of Inner lnner Expansion as E a 0 for z in the overlap region So 71n Un for 71 07172 as e a 0 ie g a 00 So the matching conditions are lim 7 Un 0 for 71 071727 440 4700 METHOD OF INTERMEDIATE LIMITS see Reiss notes p 48 and 49 and Hinch7 section 5147 p 57 Recall the outer variable is z and the inner variable is g zZJ6 De ne an intermediate variable 77 77 441 where Me a 0 as e a 0 and Md756 a 00 as E a 0 That is7 Me goes to zero slower than gt6 goes to zero as e a 0 As a result ltlt 95 ltltg 95 77 7 77 7146 456 as E a 0 For example7 if xe7 then 77 The intermediate limit has 77 xed and E a 0 So in order to match7 write everything inner and outer expansions in terms of the 442 intermediate variable7 777 x 77 and let 6 a 0 Note that z E77 a 0 as E a 0 with 77 xed and g ZE77 6 a 00 as E a 0 with 77 xed August 21 2008 Industrial Math Notes DM Anderson 27 426 Composite or Uniform Expansion Here we put the above pieces together to obtain a solution that is uniformly valid in x Basically7 we sum the inner expansion and the outer expansion and subtract the overlap or matching part So u7 6 Outer Solution lnner Solution 7 Matching Solution 443 An alternative is Van Dyke7s matching principle in which we equate the outer limit of the inner solution expressed in terms of the outer variables m xed as e a 0 and the inner limit of the outer solution expressed in terms of the inner variable xed as e a 0 43 Where is the Boundary Layer See Bender and Orszag7 p 42274267 Hinch7 Section 5177 Fowler7 Sections 41743 Consider the boundary value problem 6y axy bzy 07 444a y0 A7 444b y1 B7 444c where 6 ltlt 17 a and bx are given functions For an outer solution assume y 6 y0x 6y1x 62y2z 445 At 01 we nd awn6 bzyo 0 446 The solution to this problem is y0z coexp i m dt 447 For simplicity we shall assume that a 31 0 for z E 071 Since this is a singular problem note that the leading order out problem is only a rst order ODE while the original equation is second order we expect a boundary layer How do we determine the locations of the boundary layers Our approach here will be to test different locations and see under what conditions a boundary layer there appears possible First try z 0 Here we introduce the inner variable X r This gives 62Y Y b6XY 0 448 August 21 2008 Industrial Math Notes DM Anderson 28 We assume that 105X 7 449a b 06X 4 44910 For a balance we choose 5 6 assuming that 10 is nite This gives Y 1EXY EbEXY 0 450 Seek a solution of this rescaled equation of the form YX7e Y0X6Y1X 451 The 01 inner problem in this case is Y0 10Y0 07 452 where we note that the coef cient on the rst derivative term has been expanded and is now a constant The general solution is Y0X A0 Boe ale 453 We apply the boundary condition Y0 A and nd Y0X A 7 301 7 5 40 454 Now7 if 10 gt 0 then this is a valid solution that is possible to match to the outer solution note that the exponential decays as we approach the outer region7 X 7 00 and a boundary layer is possible here On the other hand7 if 10 lt 0 the exponential grows as X 7 00 and this solution cannot be matched unless B0 0 and there is no boundary layer So if 10 gt 0 then a boundary layer on the left at z 0 is possible if 10 lt 0 then a boundary layer on the left is not possible Next7 consider 1 1 Here we introduce the inner variable X 1 7 This gives 6 11 7 5X my 7 52 5 Assuming that 11 74 0 and 11 is nite we get a dominant balance if 5 6 So Y b176XY 0 455 Y 7 117 6XY 6b1 EXY 07 456 with boundary condition YX 0 B note 1 1 corresponds to X 0 August 21 2008 Industrial Math Notes DM Anderson 29 The 01 inner problem in this case is Yg ia1Y0 0 457a Y00 B 457b This has solution Y0X 4 B 4 B01 4 5 10 458 where B0 is an arbitrary constant Now if a1 gt 0 the inner solution blows up as X a 00 so no boundary layer is possible If a1 lt 0 the inner solution can be matched as X a 00 so a boundary layer is possible So if a1 gt 0 there is no boundary layer on the right if a1 lt 0 boundary layer on the right is possible EXAMPLE 41 Consider the case where ar gt 0 so that based on the above arguments there is a boundary layer on the left at z 0 As we have seen the outer solution is A1 dt 459 where ywt1 B satis es the boundary condition on the right Also we have seen that youtx N B exp the inner solution is Y X N A 4 B0 1 4 WM 460 where X ze and A satis es the boundary condition on the left To do the matching we write the outer solution in terms of the inner variable and eccpand in 6 We get ywt1 B exp A dt Bexp 1 6 0 a dt O X 461 So matching giues Xgig y m X 4 mm 4 o 462 This implies 1 W 1 A 4 B0 4 B exp0 mdt 4 0 463 Then we can identify ymatch B exp 1 dt 464 0 at The resulting uniform eccpansion to leading order in 6 is A1 dt A 4 B exp 01 dtgt WWW 465 N B exp August 21 2008 Industrial Math Notes DM Anderson 30 EXAMPLE 42 Consider the case where ap lt 0 so that based on the above arguments there is a boundary layer on the right at z 1 The outer solution is 10 dt 7 466 where th0 A satis es the boundary condition on the left The inner solution is ywtp A exp Y X N B 4 B0 14 WW 467 where X 1 7 p6 and B satis es the boundary condition on the right The matching solution is ymatch A exp 4 dt 468 The resulting uniform eccpansion to leading order in 6 is w b t 1 b t 7 th B 7 Aexp 7 th e 117m5 469 0 at 0 at EXAMPLE 43 See Bender and Orszag p 431435 Use the method of matched asymptotic eccpansions to solve the problem N Aexp 6y 1 xy y 07 47021 410 17 470b y1 1 4700 Here ap 1 z gt 0 for z E 071 so based on the previous analyses we eppect a boundary layer on the left near z 0 Outer Solution Assume a solution of the form y0z 6y1z 62y2z 471 Plugging this eccpansion into the governing equation gives at 01 Hm61 07 472a 4101 1 472b This has solution 2 7 473 W 11 y lt gt At 06 we nd 1 9044 41 743 741 x 47421 4111 0 474b August 21 2008 Industrial Math Notes DM Anderson 31 Guessing a solution of the form yl ll1 z3 bl1 leads to the solution 2 1 y1x W 7m 475 At 062 we nd 1 7 7 77i y y 91 148 14837 1121 0 476b 476a Guessing a solution of the form yg lg1 z5 b21 z3 021 leads to the solution 6 1 1 W W W W39 7 Inner Solution We let X 956 and transform the equation to the new variable and nd 6 1 EY 31 6XY Y 0 478 To obtain a dominant balance we choose 6 6 so Y 1 6XY EY 07 4798 Y0 1 479b Assume a solution of the form YX Y0X 6Y1X 62Y2X 480 The leading order inner problem is given by Y0 YO 7 07 4818 Y00 1 481b This has solution Y0X 40417405 X7 482 where A0 is an arbitrary constant The 06 inner problem is given by Y1 Yl 7Y0 7 XYO 740 7 1 7 40eX X17 A0e X7 483a Y10 0 483b The general form of this solution is Y1 A1 B1e XClXe XD1X2e XE1X Plugging this into the di erential equation and determining eoe eients giues Y1X 7 411 7 5X 7 40X 4 40 71X25X7 484 August 21 2008 Industrial Math Notes DM Anderson 32 where A1 is an arbitrary constant Similarly we nd that Y2X 7 X2 7 2X 17 A0 X le X 7 X2 2X 7 41 X 7 xzaX 4 421 7 5X 485 Matching solutions and matching regions39 At 01 we eccpand the outer solution 4108 N 21784824 486 The inner solution has the form Y0X N A0 17 A0e X N 407 487 provided e X is su ciently small Then y N 20lt677 48821 Y N 40405X7E7EX 488b so A0 2 and we require that z lt 1 from the outer solution e X ltlt 1 from the inner solution so that these together give the matching region for the 01 match 6 lt z lt 1 489 At 06 we eccpand the outer solution 4118 N 21738717x73787 490 The inner solution has the form Y1X N 4117 5X 7 40X 4 40 71X25X N 41 7 40X 491 Then y N 217 n 6 06276727 49221 Y N A0 4 A1 7 40Xe 0e X7 627 X627 X262 49210 Matching requires A0 2 as before and A1 32 We require that 2 ltlt 6 from the outer solution and e X ltlt 6 from the inner solution so that the matching region for the 06 match is 6 1nd lt z lt 612 493 At 062 we eccpand the outer solution 4128 N 61757173x7170x 494 August 21 2008 Industrial Math Notes DM Anderson 33 The inner solution has the form Y2X N X2 7 2X 17 A02X 7 X2 7 A1X A2 495 Then 2 3 11 21 2 3 2 2 3 21 ixz 6 7 Z6 O6 7E 76 7x 7 496a N A0 A1 7 AOXE A2 7 AlX 7 7 62 Oe X7 637 63X7 63X27 ESXS 496b Matching requires A0 2 A1 32 as before and A2 214 We require that 23 ltlt 62 from outer solution and e X ltlt 62 from inner solution so that the matching region for the 062 match is el lnel lt z lt 623 497 Note that from the inner solution we actually have e X ltlt 62 but this is just 7iE ltlt ln 62 2 lne Notice that the matching region is getting smaller at higher orders In fact for some small values of 6 this region fails to eccist since el lnel gt 623 Finally if we de ne the matching functions 99mm A0 7 AM Ao z7 498a ylnatch A1 7 Aw 7 2A0 498b alum A27 4980 we can write down a uniform appropimation to 062 yam yea em 622206 YoX 6Y1X E2Y2X ygnatch Eyvlnatch Ezyvznatchl 499 EXERCISE 41 Use the method of matched asymptotic eccpansions to solve the problem ey ay by 07 4100a y0 047 4100b y1 B7 4100c where a7b gt 0 Use the above analyses to anticipate the boundary layer location 5 Multiscale Methods Recall the damped harmonic oscillator problem d22dz 0 51 dtz Edt 95 7 7 39a 950 07 51b dx August 21 2008 Industrial Math Notes DM Anderson 34 where E is a small parameter 6 ltlt 1 In this section we will solve this problem using multiscale methods First7 however7 let us recall the results obtained by applying regular perturbation meth ods We assume that the solution can be expressed in the form t 6 z0t 61t 62z2t 52 The 01 problem is given by dzz Lzo E 20 0 07 53a z00 07 53b dig 7 0 1 53 dt lt gt lt c where L E dZd2 1 is a linear differential operator The general solution to this problem is z0t A0 cost Bo sin t Applying the initial conditions gives z0t sint 54 The 06 problem is given by dx Lxl 7 7261707 55a x10 07 55b d l gm 7 0 55c The general solution to this problem is z1t itsint A1 cost B1 sin t Subject to the initial conditions we nd z1t itsin t 56 Our solution7 then7 based on this regular perturbation expansion is zt sintietsint 57 As we have observed earlier7 the problem is that t is not bounded Therefore7 the expansion breaks down when t 06 1 The expansion is not uniformly valid for all t In order to resolve this problem we can consider the exact solution available in this problem but not generally available in other problems andor ii use physical reasoning to identify the real long time behavior of the system Thinking about the latter7 in our problem there are two time scales One is an 01 time scale associated with the oscillations and the second is a long time scale 06 1 related to the slow decay of the envelope of the oscillations In this particular case7 let us examine the exact solution to gain further insight August 21 2008 Industrial Math Notes DM Anderson 35 EXACT SOLUTION TO EQUATION 51 Let xt 6quot This implies m2 26m1 0 Then7 m1 7e i Z39V1 7 62 This leads to the general solution t 6 A cos V1 7 621 B sin V1 7 6 58 Applying the initial conditions gives the exact solution Mt 6 5 sin V1 7 621 59 V177 T This solution shows the presence of a time scale at related to the exponential decay amplitude of the oscillations We can also see the 01 time scale associated with the oscillations Note that we can recover the result obtained by the regular perturbation expansion by taking at small and expanding the exact result in a Taylor Series As we have learned already7 this is valid for limited time only t lt 016 Note also that if we expand V1 7 62 we nd that the oscillatory term has the form 1 sinV1762tsinltt7 562tgt 510 Here we can see that many more times scales can be observed in addition to the 01 time scale t we have 6217 64t7 etc How can we try to predict these time scales without knowledge of the exact solution Consider again the equation dzz dx tQEE i 0 511 There is clearly a balance of terms 1 and 3 as e 7 0 giving undarnped oscillations o However7 as noted above7 we expect term 2 to become important eventually7 so we may ask What balance includes term 27 To address this7 rescale time and let 739 eat Here we anticipate that t will be large so that 739 will be 01 This gives 6277 261M957 z 0 Balancing the rst two terms gives 20 1 047 or 04 1 This does not lead to a dominant balance in the previous sense 7 narnely7 the third term in this case dorninates Balancing the second and third terrns gives 04 71 which is also no good The problem is that 01 time scale balance has not gone away 7 there is just a second time scale sornehow superirnposed in the problem August 21 2008 Industrial Math Notes DM Anderson 36 A solution to this problem is to introduce simultaneously two time scales into the problem These will be treated as mathematically independent variables Speci cally7 we treat 739 eat and t as independent variables 512 The time derivatives in the governing equations transform as d 6 d7 6 a a a axafaa a 513 so that d3 7 63 61 514 dt 7 6t 6 677 39 a dzz 3 a 32x Zadzz 7 kQE ataT 6 7 Note that z is now treated as a function of two variables7 t and 739 and so we use partial differentiation Plugging these expressions into the original equation gives m 26 6277 26 96 6 z 0 515 Term 1 Term 2 Term 3 where subscript indicates partial differentiation with respect to the speci ed variable The balance that ends up working in this case is a balance between term 260 in Term 1 and the term 26 in Term 2 This argument does not always work the detailed explanation is that the time scale 739 is chosen to bring in 739 derivatives into the expansion at the point where secular terms 7 eg 6t sint 7 arise 51 Twotime Method First we solve equation 51 with a two time method and then later apply a three time method We introduce a slow time 739 6t but retain the 01 time variable t dx 7 3a 3a 5 16 dt 7 at 6677 39 a dzz 32x 32x 262 267 E 7 7 516b dt2 3t 3th 372 7 Plugging these expressions into the original equation gives at 26th ezz 26 z 6 z 07 517 We seek a solution to the above transformed equation in the form t 6 x0t7 739 61t7 739 622t7739 7 518 August 21 2008 Industrial Math Notes DM Anderson 37 where we require that the xjt7739 are bounded Enforcing that the 7 are bounded is one of the key steps in the method7 as we shall see We plug this expansion in and obtain 012 6 11 20t7gt 62 21 21t7 077 08 26 mt 6 1 07 08 0 6x1 622 063 07 519 Notice that in the initial conditions we have 50070 61070 622070 07 5205 007 0 E x1t070 20707 0 1 520b At 01 we nd Lzo E zott x0 07 521a 50070 07 5215 50070 1 5210 The solution to this problem is given by x0t7739 AT cost BT sin t7 522 where the initial conditions imply that A0 0 and B0 1 Note that the difference between this result and the 01 result from the straightforward perturbation expansion with a single time variable see equation 54 is that the coef cients now depend on the new7 slow7 time scale The functional forms of AT and BT are chosen at the next order in the expansion so that secular terms terms like tsin t in the solution expansion can be avoided ie so 1 is bounded At 06 we nd Lxl 72 7 2x077 523a 51070 07 5235 107 0 70707 0 523c lnserting the 01 solution 522 leads to the equation Lxl 2 fl739 AT sint 7 2 BT BT cost 524 Now7 to avoid secular terms in the solution z1t7739 we must choose AT and BT such that the coef cients of sint and cost in the above equation vanish Therefore7 we require that A A 07 525a A0 07 5255 3712 3th2 371 372 3722 August 21 2008 Industrial Math Notes DM Anderson 39 Plugging these expressions into the original equation gives 96 z 269651 96 62 29 96m 2le 263 96 9672 6 967m 0533 We seek a solution to the above transformed equation in the form t 6 x0t7 73917 7392 61t7739177392 1L7 534 where now each Q is assumed to depend on three independent variables and are required to be bounded At 01 we nd Lzo 07 535a 95007 07 0 07 5355 5007 07 0 1 5350 The solution to this problem is given by x0t7 7177392 A739177392 cost B739L7 7392 sin t7 536 where the initial conditions imply that A07 0 0 and B07 0 1 At 06 we nd Lxl 72mm 7 27 537a 95107070 07 5375 5107070 407107 07 0 5370 lnserting the 01 solution 536 leads to the equation Lzl 2 A71739177392 A73917T2 sint 7 2 B71 71772 BT177392 cost 538 Now7 to avoid secular terms terms proportional to sint and cost lead to terms like tsint and tcost in the particular solution for 1 in the solution 1t7739177392 we must choose A739177392 and B739177392 such that the coef cients of sint and cost in the above equation vanish Consequently7 we require that AT1 A 0 and B71 B 0 which imply that A739177392 17392e n7 B739177392 b7392e 7 539 subject to the conditions 10 0 and b0 1 Given the above choices7 and noting that x07107 07 0 A71070 07 it follows that 1 satis es Lxl 07 540a 95107 07 0 07 5405 5107 07 0 0 5400 August 21 2008 Industrial Math Notes DM Anderson 40 This equation has solution 1 0 there are no driving terms At 062 we nd that 1amp2 2 hm ht 2th 07171 20n7 5413 520 07 0 07 5415 5207 07 0 75171 0070 7 50720 07 0 5410 Simplifying the right hand side of the above equation gives ng 72 7AT2 sint BT2 cos t 7 A7171 cost B7171 sin t 7 2 ATI cost BT1 sin t7 cost 72BT2 7 A7171 7 2An sint 2A72 7 E7171 7 2Bn 542 To avoid secular terms in 2 we choose 72B72 7 Ann 7 2A71 07 242 7 BM 7 231 0 5435 Using equations 5397 these translate to 72b a 07 544a 2d b 0 5445 These can be combined to give 4b b 07 which has general solution b7392 61 cos73922 g sinT22 The initial condition b0 1 implies 61 1 Since a 2b we have 17392 7sinT22 62 COST22 Then 10 0 implies g 0 So a 7 sin739227 b COST22 545 Therefore7 x0t7 73917 7392 6 71 sintcos73922 7 cos tsin73922 546 Using the trigonometric identities in equation A1 gives 0t7739177392 6 71 sint 7 739227 5 sint i 62t2 547 Compare this with the expansion of the exact solution as described in equations 59 and 5 10 It is informative to examine the plots ofthe various approximate solutions to this problem and compare with the exact solution Speci cally7 we plot a hierarchy of solutions sir11 7 sint 7 etsint 7 6 5 sin t7 6 5 sint 7 62t27 548a regular expansion regular correction twotime three time in Figure 2 August 21 2008 Industrial Math Notes DM Anderson Multiscale Method 2 01 l O exact solution 7 leading order regular regular correction 7 multiscaletwo time multiscalethree time Figure 2 This gure shows a comparison of solutions based on different rnulti scale methods and the exact solution of equation 51 The plots in this gure have used 6 01 Clearly the regular perturbation expansion results fail to capture the long time behavior while both rnulti scale results agree very well with the exact solution 41 August 21 2008 Industrial Math Notes DM Anderson 42 53 Nonlinear Problems Up to this point we have considered linear problems For nonlinear problems exact solutions are typically not available and our ability to nd approximate solutions by way of multi scale methods is therefore of even greater importance Consider the nonlinear ODE dzu du 2 dtz EltEgt u 0 54921 u0 0 5495 du 70 l 549 d 0 Note that in the context of the massspring system this problem does not make physical sense since it would correspond to a damping term of the form FD 76d2udt2 which is always negative even when dudt lt 0 The oscillations do not decay as one would normally expect However it will demonstrate certain points that need to be under consideration in nonlinear multi scale problems We will address this problem with a two time multi scale method We assume that the time scales of interest are t and 739 at As we shall see shortly this is not a correct assumption The important lesson to be learned here is how to identify a better assumption With time scales t and 739 treated as independent we have du 6n 6n 5 50 7 7 7 7 dt 6t 6677 a dzu 321i 321i 321i 7 7 2 7 2 5505 dt2 T at Eon37 E 672 Then the governing equation 549 can be written um 261m 621177 E 261ml 6211 u 0 551 or um u e 2 u 62 u 21W 6 0 552 subject to the boundary conditions u0 0 0 and ut0 0 eu0 0 1 The 01 problem is given by Luo 0 553a u00 0 0 5535 u0t0 0 1 5530 Notice that in this case since the nonlinear term is multiplied by E the leading order problem is linear We will nd that all the subsequent problems are linear as well that is the nonlinear August 21 2008 Industrial Math Notes DM Anderson 43 problem is to be approximated by solving a sequence of linear problems The solution7 as before7 is given by uot7 739 A0739 cost B0739 sin t7 554 where A00 0 and B00 1 The 06 problem is given by LuL 72110 7 uof7 555a u107 0 07 5555 u1t070 7uOT07 0 7A60 555c In order to identify the solution ul we need to simplify the right hand side terms in this expression We nd LuL 72 7A6 sint B6 cost 7 7A0 sint Bo cos t27 2 A6 sint 7 B6 cos t 7 A6 sinzt 7 2A0B0 sintcost B6 cos2 t 7 1 1 1 1 2A6sint 7 2B6 cost 7 A6 7 E cos 2t AOBO sin2t 7 B6 i cos 2t7 1 1 2A6 sint 7 2B6 cost 7 A6 B3 A6 7 B3 cos 2t AOBO sin 2t7 556 where we have used trigonometric identities in equations A1 and A2 Now7 to avoid secular terms like tsint and tcost in the solution ul we require that A6 0 and B6 0 This suggests that the coef cients A0 and B0 do not actually depend on the time scale that we had assumed that they did With this result and application of the initial conditions we nd A0 0 and B0 1 so the solution at 01 takes the form u0t7739 u0t sint 557 While this is a clue that there is a problem7 let us continue by nding the solution uh which satis es 1 1 LuL 7i 7 Ecos 2t7 558a 11107 0 07 558b u1t070 0 558c The particular solution has the form ulp co 01 cos 2t Plugging this into the governing equation shows that co 712 and 01 16 Then the general solution is 1 1 u1t7739 7i 8 cos 2t A1739 cost B1739 sin t 559 The boundary condition 11107 0 0 implies that A10 13 and 11107 0 0 implies that B10 0 The functions A1739 and B1739 can be determined at the next order August 21 2008 Industrial Math Notes DM Anderson 44 The 062 problem is LUZ 2u1t7 2u0tu1t 077 21101107 7 5603 1120 0 07 560b uZ 07 0 7u17070 7A10 560c We simplify the right hand side by rst plugging in the known formulas for no and ul 1 Luz 72 7Af1L sintBicost costlt7 sin2t7A1 sintB1 cost 2 2Af1L sint 7 2BL cost g costsin 2t 2A1 costsint 7 2B1 cos2 t 561 Before we can identify the form of the solution uz we write all of the terms on the right hand side in terms of sinpt and cospt where p is constant Speci cally7 note7 for example that from equations A2 we have 1 1 1 costsin 2t i sin2t 7 t sin2t t i sint E sin 3t 562 2 1 1 Luz 2A Lsint7BjLcostnLglt sint sin3tgtAlsin2t7B11cos2t7 1 1 sint 2Af1L cost 7B 7 B1 7 B1 cos 2t A1 sin 2t g sin 3t 563 Lead to secular terms Do not lead to secular terms The above equation identi es the terms on the right hand side that do and do not lead to secular terms in the solution u2 In order to avoid secular terms we choose 2Af1L 0 and B1 0 note that the latter satis es B10 0 derived earlier Recalling the boundary condition A10 13 we nd that A16 J 7 564 We could proceed and solve for uz but our primary interest in the 062 problem was to determine the unknown functions appearing in the 06 solution7 which we have now done Let us examine the consequence of these choices in the resulting solution u uo eul We have t 1 1 t 11 2t02 565 us1n e 777739 cos 77 7cos 6 3 6 2 6 We can see that this expansion is nonuniform due to the term 6739 cost That is7 the expansion breaks down when t 06 2 Although our objective of nding a uniform solution has failed7 it suggests a solution to the problem the need for a dz erent time scale 621i August 21 2008 Industrial Math Notes DM Anderson 45 Let us know retry the problem with the time scale It and 739 621 The equation 549 expands in these two new variables as um 2627 641177 E 262mm 6411 u 07 um u 6113 2627 263mm 064 7 07 As before7 the 01 problem is given by Luo 7 07 568a 1100 0 07 568b u0t07 0 1 568c The solution7 again7 is given by uot7 739 A0739 cost B0739 sin t7 569 where A00 0 and B00 1 The 06 problem is given by LuL 7u0t2 570a 1110 0 07 5705 u1t07 0 0 5700 Substituting in using the 01 solution 569 we nd that LuL 7 A3 cos2 t 2A0B0 cos tsint B3 sin2 t 7 1 1 7 7 33 7 5 cos 2t Ag 7 33 7 AOBOsin 2t 571 The particular solution takes the form 1 ui 7013 B3 1 01 cos 2t d1 sin 2t 572 Note that none of these terms are secular terms they are bounded as t 7 00 In order to satisfy the ODE we require 7401 01 B3 7 A2 and 74d1 d1 7A0B0 so that 01 7B3 7 A6 and d1 AOBOB Then the general solution to the 06 problem is 1 2 2 1 2 2 1 7 A0 B0 7 8B0 7 A0 cos 2t gAOBO s1n 2t A1739 cost B1739 sint 573 U1 The initial conditions can be used to determine conditions on the functions A1 and B1 August 21 2008 Industrial Math Notes DM Anderson 46 The 062 problem is given by Luz 72u0tu1t 7 2110 E RHSZ 574a u20 0 07 5745 11207 0 71104070 574c Our primary interest here is to obtain conditions for the functions A07 and B07 appearing in uo These conditions will be determined by requiring that terms in RHSZ do not lead to secular terms in u2 Consequently7 we shall identify terms in RHSZ that contribute to terms proportional to sint and cost and group all other terms such as sinpt and cospt where p 31 1 as non resonant terms In this calculation we will not actually nd uz and so we do not need to keep track of the speci c form of these non resonant terrns RHS 1 2 722 7A0 sint Bo cost lt Bg 7 A6 sin 2t gAOBO cos 2t 7 A1 sint B1 cos t 7A6 sint B6 cos t7 1 2 sintsin 2t 7 A0Bg 7 713 sintcos 2t 7gA6B0 7 cos t sin 2t 1 2 B0Bg 7 713 cos tcos 2t E71033 A0A1 sinZt 7 A0B1 AlBO cos tsint B0B1 cos2 t7 7 A6 sint 7 B6 cos t7 1 1 78710093 7 A6 cost 7 cos 3t 7 gAgBO sin7t sin 3t 7 1 1 8B0Bg 7 A6 sint sin 3t gAOBg cost cos 3t 7 A6 sint B6 cos t7 non resonant terrns7 1 1 1 1 7 gm 7A6 7 g71350 7 85053 7 713 7 cost 33 7 8710093 7 Ag 7 gAOBg non resonant terrns 575 Therefore7 to avoid secular terms in the solution uz we require that A0 and B0 satisfy 1 A6 8377137 B3 576a 1 B6 78A0A 33 576b which are subject to the conditions A00 0 and B00 1 We can solve these equations by multiplying the rst equation by AO7 multiplying the second equation by Bo and adding the results together to obtain A0A6B0B6 0 577 August 21 2008 Industrial Math Notes DM Anderson 47 This implies that Ag Bg 07 or A3 B3 is constant By the initial conditions we must then have Ag B3 1 Consequently7 the equations for A0 and B0 reduce to 1 Ag 8307 578a 1 Ba i Ao 578b Combining these gives A3 A036 so that A0 is given by 1 l 1 A0739 acos 8739 bs1n 8739 579 Applying the initial condition A00 0 implies that a 0 Then B0 6AJ bcos 76 and initial condition B00 1 implies that b 1 Therefore7 A0739 sin 73 580a 1 B0 739 cos 8739 580b Finally7 the leading order solution u0t7 739 can be written as l 1 1 l u0t7739 sin 8739 costcos 8739s1nt7 1 sin ltt 581 Figure 3 shows this solution compared with a numerical solution to the original problem 5 49 and also with the solution obtained using a multiscale method with timescale 6t 54 Exercises EXERCISE 51 See Hinch eccercise 71 p 120 see Bender and Orszag 17 5527554 for solution Obtain an asymptotic approrimation for z to 01 which is valid fort 06 1 to the solution of dzz dx 3 i 0 5 82 dt2 6 dt 95 7 7 39 a 950 17 582b dx 0 0 582 dt 0 Note that this problem can be related to a massspring system with nonlinear damping ie FD 7Bdzdt3 F5 in and Mdedtz FD F5 August 21 2008 Industrial Math Notes DM Anderson 7 numerical solution 7 7 solution based on at 7 solution based on e2t i u H l 100 150 Figure 3 This gure shows a comparison of solutions based on different multi scale methods and a direct numerical solution of equation 549 The plots in this gure have used 6 02 The solution based on timescale 739 at is clearly nonuniform7 with amplitude growing with t The solution based on 621 is uniform in the sense of its amplitude but the phase eventually can get completely out of phase 48 August 21 2008 Industrial Math Notes DM Anderson 49 EXERCISE 52 See Hineh erereise 72 p 120 Obtain equations for the drift in amplitude and phase of the solutions to d2 d d7 6x2 71 1 6hz ECOS t 583 Assume h 01 as e 7 0 Hint write the leading order solution in the form 0 R0739 sint 60739 where 739 6t 6 Basic PDE Theory See Zauderer27 Evanssg Kevorkian and Cole 7 p 330 Consider the linear 2 Order PDE Lu E Aum 2Bumy Cuyy Dum Euy Fu G 61 principle part We can rewrite the principle part as the product of two rst order differential operators and some rst order terms Aum 2Bumy Cuyy 6m 7 may 6m 7 wgdy u 3me 7 wldng 3211i7 62 where Wz Wllv 63a wle NONE Multiplying the rst of these by oil and combining it with the second equation and then repeating the same with W2 shows that oil and Mg satisfy sz 2Bw o 0 64 This implies B 1 u 7ZiZVB27AO 65 So the nature of the 2nd order part of the differential operator Lu changes with the sign of B2 7 AC If B2 7 AC gt 0 then oil and Mg are both real and unequal Here Lu 6m 7 may 6m 7 wgdy u rst order stuff 66 leads to wave like motion and this equation is of hyperbolic type For example7 if A 17 C 717 B D E F G 0 then we have a standard wave equation Lu um 7 uyy 0 67 August 21 2008 Industrial Math Notes DM Anderson 50 Anotherexample isA17 B7327 027 DEFG0so um 7 311m 2uyy 6m 7 26y6m 7 3yu 0 68 This equation has solutions of the form u cof y 0192z Kevorkian and Cole p 331 point out that real characteristic curves exist A point disturbance solution of Lu 0 in uences a restricted portion of the domain7 bounded by characteristic curves Two initial conditions are required on a space like boundary One boundary condition is assigned on a time like boundary If B2 7 AG lt 0 then wl and 102 are both complex conjugates For example7 wlg it gives 61 726261 my u um um Vzu 0 69 This is Laplace7s equation D E F G 0 More generally7 umuyy Dum EuyFu G 610 If D E F 0 we have um uyy G Poisson Equation or if D E G 0 we have um uyy 7Fu Helmholtz Equation These equations are of Elliptic Type Kevorkian and Cole p 331 state No real characteristics A point disturbance solution of Lu 0 in uences the entire space One boundary condition is prescribed on a closed boundary77 If B2 7 AG 0 then wl 102 7BA Then Lu 6m 7 wldy2u Dum Euy Fu G 611 This equation is of Parabolic Type Examples include the heat equation E 717 B F G 0 7uy um 07 612 the heat equation with transport D E 717 B F G 0 7uy 7 um um 0 613 These become the more familiar equations if one replaces y with t A nal example is ut7Dum umFuG7 614 which includes the time derivative term7 a transport term7 a diffusion term and two source terms Kevorkian and Cole p 331 state that this is the dividing case between hyperbolic and elliptic type equations Here real characteristics coalesce Half the space is in uenced by a point disturbance One initial condition on a characteristic arc or one boundary condition is prescribed7 August 21 2008 Industrial Math Notes DM Anderson 51 7 Green s Functions Methods for ODEs References for Green7s Function Problems Boundary Value Problems of Mathematical Physics I IL lvar Stakgold Green7s Functions and Boundary Value Problems7 lvar Stakgold Methods of Techniques of Applied Mathematics77 Friedman Methods of Math ematical Physics 1 Courant and Hilbert A First Course in PDE7 Weinberger Advanced Calculus for Applications7 Hildebrand See also Hildebrand Methods of Applied Mathe matics 7 pp 2287229 Here we consider a boundary value problem of the form d du u g qzu 7 7Fx7 71a du a1u l 0 atza7 71b du 0421i 62 0 at z b 71c A reformulation of this problem involves a Green7s Function G The function G is given by 7 G1z zlt GWl 7 GM 95 gt 72 This function has the following four properties 0 G1 satis es G1 0 when x lt zz39 and G2 satis es G2 0 when x gt g o G satis es homogeneous boundary conditions at z a and z b ie a1G1a 616G16za 0 and 042G203 o G is continuous at z that is G1 G2 We call this the continuity condition 0 The derivative of G has a discontinuity of magnitude 71p at the point z that is7 G 2 7 Gi 71p We call this the jump condition Another way to de ne the Greens Function G is as the solution to G 76x 7 g 73a dG 041G 61 0 at z 1 73b dG 042G 62 0 at z b7 73C where 6Q 7 g is the Dirac Delta Function see below for more details Then7 the solution to equation 71 is given by 17 W a Gltxl F d 74 August 21 2008 Industrial Math Notes DM Anderson 52 A heuristic de nition of the Dirac Delta Function is 6ltze gt000 75 such that b 7 1 if a lt g lt b a 6m 706 T 0 otherwise 7 76 or b 7 45 if a lt g lt b a pawa 7 06m 7 0 otherwise 7 7397 holds The latter is typically known as the sifting property Technically7 this statement is not rnathernatically correct but as written it has great utility These ideas can be put on a more solid foundation using the theory of distributions see Section 71 71 Comments on Distribution Theory DEFINITION 71 A distribution is any continuous linear functional which takes a function as input and returns a number or a rule tx de ned on 08 the space of in nitely di er entiable functions that vanish outside some nite interval which we write as We gt tltzgtltgtdz 78 ltt96ltp96gt awash 79 for Md 6 08 DEFINITION 72 If the rule ltt7 pgt is generally by a locally integrable function tz then it is called a regular distribution otherwise it is called a singular distribution In the latter case tz is generalized or syrnbolic DEFINITION 73 The 6 distribution is given by the rule W907 M96 0 710 This is a singular distribution and cannot be interpreted as a classical integral since 6p is not locally integrable August 21 2008 Industrial Math Notes DM Anderson 53 DEFINITION 74 Two functions t1 and t2z locally integrable or generalized if t7 p is a singular distribution are distributionally equal whenever they generate the same distribution That is if ltt17 M96 ltt2967 M96 f0r 11W 6 08 711 then t1 t2 712 is a distributional equality which is di erent than the usual equality Heuristieally this is like being equal almost everywhere For exarnple7 dZG W 6 97 73913 is a distributional equality Another example of a distributional equality is H z To see this note that l H ltzgt ltzgtdz sayammo 7 l Hzltp 96d 7 0 Winch e we 00 714 and also that L 6ltpdx ltp0 715 In general7 distribution theory will allow us to interpret things that have no Classical interpretation 711 Distributions and Differential Equations lf tz is such that t7 p f7 p for all p 6 08 then we say that tx is a solution to the distributional differential equation t 716 This is a distributional equality an nth order ODE for the unknown distribution t o If tx is nitirnes differentiable then it is a Classical strong solution of equation 7 16 o If tz is locally integrable but not nitirnes differentiable then it is a weak solution of equation 716 0 Then lt95t7ltPgt lt07ltpgt0 ltt795ltPgtltt7 95lt gt7 719 for all p 6 08 So we d like to nd tx such that ltt7 mow 0 Try tx c so ltc ltwgt gt cltwgt dz memo o 720 for all p 6 08 So tx c is a classical solution Also try the Heaviside function tx Then ltHltwgtgt lHltzgtltz gt dzA ltz gt dzmom 721 for all p 6 08 So tx is a weak solution since is locally integrable but not n 1 times di erentiable So tx c1 c2H7 722 is a solution make up of a classical solution and a weak solution x 0 Then ltx2t 7ltpgt ltPgt t7 952ltPgt7 723 for all p 6 08 As before tz e and tz work Now try tz Then lt6ltzgtltz2wgt 6ltzgtltz2wdz 29677 W So tx 6a is a distributional solution since it is not locally integrable So tx 01 02H 036z7 725 0 724 m0 is a solution made up of a classical solution a weak solution and a distributional solution We eventually want to consider a problem like dZG W 7595797 726a G0l 0 726b Gll 0 726c Then ltGH7ltpgt 7 gt7 G7ltp 727 More generally this is lt G7ltPgt ltG7 Wgt 77707 728 with suitable boundary conditions on Note that to solve for G we can use dZGdx2 0 if z gt g or if z lt We expect G to look like a delta function7 G to look like a Heaviside function with the jump at z 7 and G to be continuous So to obtain G7 we solve for G in z lt g and in z gt g and then apply continuity on G and a jump condition on G This is shown in the following example 72 Introductory Example This problem is from Stakgold7 Vol I Here we consider the static de ection of a string by an applied load We let T be the unde ected tension in the string and f the load per unit depth a given applied load Assuming small de ections we have the boundary value problem 2 7T fx7 0 lt z ltl 729a u0 07 72 U 07 729c August 21 2008 Industrial Math Notes DM Anderson 56 where and u are continuous on OJ ie is continuously differentiable ls this BVP well posed Well posed in the Hadamard sense means the solution exists7 ii the solution is unique7 iii the solution depends continuously on the data To solve this problem we introduce a Green7s Function Problem 75375 0ltxltl7 7302 G0l 0 730b Gm 0 7300 In this problem7 the Greens Function can be thought of as the solution to the case in which an extremely high load is applied at z Here G represents a de ection which is continuous but whose derivative dGdx is not continuous at z We solve the Greens Function problem in the following way Consider cm W 0 zy g 731 Which has general solution 7 Ax B if 0 S x lt Gmg Oz D ifg lt z 1 732 This is subject to the two standard7 boundary conditions G0l 07 733a GUM 7 07 733b which implies that B 0 and D iCl so 7 Ax if0 zlt Gala Czil ifgltzgz 39 73934 The conditions that determine A and C are the continuity and jump conditions described above The continuity condition gives G l G l 7 735 so Ag C 7 l The jump condition can be derived by integrating the ODE from to 5 That is7 dZG 7 i 6 i d g dxz 7 ox E 7 71 736 dx 6 dx 6 August 21 2008 Industrial Math Notes DM Anderson 57 This condition requires that C 7 A 71 Combining this result and the previous one gives 0 7 l and A 1 7 l Therefore7 the Greens Function for this problem is 7 17 lx if0 zlt Gala 14 9505 ifg lt z 1 39 73937 This is the Greens Function solution to this problem We now show how the solution to the original problem is obtained from this solution We can think of the solution in terms of the following superposition N WC Gzl if iA i 738 in which the Greens function is multiplied by a weighting function f and a horizontal spacing In the limit N 7 007 A51 7 0 we have 1 um 7 0 Gzl f d 739 A more rigorous method to derive this solution is obtained by rst de ning Lu dziu 7fx 7 40a 7 dz 7 7 39 dZG LG W 76x 7 g 740b Now consider ltG7Lugt 7 ltu7LGgt7 741 where 1 lt1 41gt E 0 pqdz 742 Using the above equations we have 1 dzu l dZG ltG7Lugt7ltu7LGgt 4 AGwdx7Ouwdx l l 0 Glt7fltzgtgtdz7 0 ultzgtlt76ltz7 gtgtdz l 7 0 Gltziggtfltzgtdz 144 743 On the other hand7 we can integrate by parts to nd dziudx 7 ludZiG dz 0 dz 1 421 4 an 0 dx dx dx ltG7Lugt 4 ltu7LGgt Ala dz du l l du dG dx G 7 odzdx 7 7 0 0 0 744 August 21 2008 Industrial Math Notes DM Anderson 58 Therefore7 l we 0 Gzl f9cd 745 01 l um 0 G lf d 746 Note that in this problem Gzl G lx can see this from the solution7 for example that is7 they are symmetric In general this is not true unless the problem is self adjoint That is7 we have Gxl G l where indicates the adjoint problem 6 is a solution to the adjoint problem lfthe problem is self adjoint G Gzl and hence Gxl G lz The solution to the BVP in terms of its Green7s Function is l we 0 Gzl f d 747 We now consider the well posedness of the BVP We rst examine the issue of existence Consider 4a 748 Now um lt17gt f d lt17gtzf d 749 Then a I lt1 7 f d jig 17 molme Am l 179mg we 7 0m flt gtd l lt1 e f d 750 lt17gtxfx7lt17gtx 7Tf07lt17jgt f950f96 751 Therefore our expression in equation 747 satis es the differential equation Also note that and um 0lG0l f d 0 7522 W 01Gll f d 0 75 A solution7 namely the one in equation 7477 exists We can further test the uniqueness of the solution and the continuous dependence of the solution on the data August 21 2008 Industrial Math Notes DM Anderson 59 8 Heat Transfer These notes are partially based on Solidi cation77 notes at NU by Worster and by Davis Also see See Fowler10 Ch 5 and Carslaw and Jaeger5 For further discussion on similarity solutions of heat transfer and related problems see Barenblattz 81 Introductory Concepts The transfer of heat in materials can occur due to o Conduction of heat heat passes through without motion of the substance 0 Convection of heat heat transfer is due to motion of the substance 7 eg in a uid 0 Radiation of heat eg electromagnetic radiation In this section we focus exclusively on conductive transfer of heat If different parts of the material are at different ternperatures7 heat ows from hotter regions to cooler regions Here we shall develop a mathematical description of this process First7 we introduce some terrninology o Ti 7 t Temperature is a function of space i and time t UNITS C or 0 Op Speci c Heat speci c heat capacity per unit mass is the amount of heat energy required to raise the temperature of a unit mass of material by 1 C without change of phase UNITS calg C This can be a function of temperature Note that 1 cal 4184 J Nm unit of energy 0 H Enthalpy is the heat energy per unit mass of material UNITS calg The derivative of H with respect to temperature is the speci c heat Op dHdT o The conductive heat ux f per unit area is proportional to the negative temperature gradient f 7kVT where k is the heat conductivity This is a constitutive model of heat conduction known as the Fourier Law of Heat Conduction UNITS of k cal 0 s cm UNITS of if calcrnz s The heat ux vector j is in the opposite direction of the temperature gradient so for example in one dirnension7 if dTdz lt 0 then q gt 0 so heat ows from hot to cold 82 Conservation Law for Heat Consider a xed but arbitrary control volume V of the material Let 7 be the outward unit normal to this volume and let 6V denote the boundary of the volume Conservation of heat requires that rate of change of heat ux of heat through the surface 81 August 21 2008 Industrial Math Notes DM Anderson 60 We assume there are no internal sources of heat such as a phase boundary This general principle can be expressed as d 7 HdV i i dS 82 dt V p 6V q For simplicity we assume that the material is incompressible so that the density p is constant Then we can write 3H idv ivch 83 V p at V q Here we have used the divergence theorem and have assumed that j is continuously differ entiable Now if we take H HT we have 63 61061 84 at dTat Paw 39 V Now7 since the control volume is arbitrary we require that so that 6T pCpEVtj dV 0 85 6T 0 7 7V a 86 p p at 7 which is a local form of the heat conservation equation With Fourier7s Law if 7kVT we obtain 6T pOpE V kVT 87 If the thermal conductivity k is constant we get the standard heat equation 6T 7 2 at HV T7 88 where H kpCp is the thermal diffusivity UNITS of H cmZs 83 Boundary Conditions Generally we are interested in solving equation 88 on some speci ed domain subject to conditions at the boundary of the domain and initial conditions For a closed domain we need a boundary condition on the surface S of the domain and an initial condition The following are common types of boundary conditions arising in the context of heat transfer 0 T ff on S prescribed temperature on the boundary Dirichlet Condition August 21 2008 Industrial Math Notes DM Anderson 61 o kVT 92 on S prescribed temperature gradient heat ux in the normal direction Neumann Condition 0 7kVT hT 7 Too on S mixed condition7 heat ux is proportional to the difference between the surface temperature and the surrounding temperature or far eld temperature outside the domain Here h is the heat transfer coef cient UNITS of h cal 0 s cmZ This boundary condition is sometimes called Newton7s Law of Cooling or sometimes a radiation boundary condition 84 Uniqueness of the Solution Here we examine the uniqueness of the solution of the heat equation 88 with associated boundary conditions and initial conditions The problem can be stated as LT WT in D amp VT igcriiuy onS7 am T m where D is a spatial domain xed in time and S is its boundary Now suppose that T1 is a solution to 89 ls it the only solution To answer this question7 suppose that T2 is also a solution to 89 De ne T T1 7 T2 Then 7 satis es 61 KV277 in D7 6t VT 71 777 on S7 810 7270 07 in D Now multiply the PDE in 810 by T and use the results 67 7 1672 TE 7 i at 7 811a TVZT v7v771v7127 811b to obtain 1672 5 2 H r7v771vr l 12 Now integrate this equation over the domain D Daizt dV KD lV397V7 lVTlZldV 813 August 21 2008 Industrial Math Notes DM Anderson 62 Since the domain D is xed in time and using the divergence theorem we have lg72dV HTVT dSin lVTlZdV 2dt D s W D ihTk h 2 2 77 TdSin v7 dvgo 814 k D D Now we de ne J TZdV20 815 D however7 1dJ ia 7 816 and Jt 0 0 since 70 0 Consequently7 Jt S 0 for t gt 0 Therefore7 the only way that these conditions can be maintained is for J E 07 or T 0 for t 2 0 Therefore7 T1f7t T2f7t so that solutions to equation 89 are unique This tells us that once we have a solution we do not need to worry about nding others 85 Solution Methods There are many situations de ned by geometries7 boundary condition combinations7 initial conditions7 etc for which one can seek solutions to this problem Some methods in clude Fourier Transforms7 LaPlace Transforms7 Separation of Variables7 Green7s Functions7 Similarity Solutions7 and Numerical Methods For now we will consider a one dimensional problem in which a similarity solution is possible 851 Similarity Solution 1D Heat Transfer From Heated Boundary Consider an in nite slab of material initially at uniform temperature To0 in the region 2 gt 0 above a xed boundary at z 0 whose temperature is xed at T3 for t 2 0 For t gt 0 the temperature in the far eld z a 00 is Too With these boundary conditions we solve 6T VT 8 17 7 H7 3t 322 As a rst step to solving this problem we examine a nondimensionalization We introduce T i T T 700 818 TB 7 TOO lt a t t 77 818b Tl 2 i 818c E August 21 2008 Industrial Math Notes DM Anderson 63 where 739 and L are time and length scales to be determined The heat equation becomes dT MT NT 7 7 819 dt L2 32 2 Hence7 we choose 739 LPH Then in dimensionless form we have 6T 62T W W7 82021 T 0 as 2 7 007 820b T 1 at 2 07 820c T 2 7t 0 0 820d Notice that we have made the system dimensionless without choosing a speci c length scale That is7 there is not an external length scale in the problem ie no domain size or boundary length This suggests the possibility of a similarity solution If we take L 2 and notice that 739 22H so m E7 2 1 821 t This tells us that there is really only one independent variable in the problem namely the combination t22 There are actually a number of ways to argue that this is an important ratio In this particular case we shall go back to the original dimensional system to develop the similarity solution A little hindsight and experience suggests that we introduce a similarity variable of the form n 7 822 and seek a solution of the form T T77 We transform the derivatives as 3T 7 den7 dT 2 i ndT 823 at duet dnzmE 2td777 39 a 6T dT 677 dT 1 7 if 7 823b 32 dn 32 dn 2M7 611 61 6L1 1 1 LdiT 8230 322 T dn 32 32 T dn dn 2 2M 7 4m dnz39 39 We plug these transformations into the PDE and nd that the problem reduces to dZT dT 7 2 7 0 824 dnz 77 dn 7 a T77 0 TB7 824b T77700 Too 824c August 217 2008 Industrial Math Notes DM Anderson 64 Note that this is an ODE reduced from the original PDE with boundary conditions There is a singularity at t 0 note that T 0 at t 0 coincides with the far eld condition T 0 as 77 a 00 We can solve this ODE by introducing an integrating factor e712 Upon integrating twice we nd that 7 T77 00 e 52ds 01 825 0 The boundary condition T0 TB implies cl T3 The boundary condition T77 oo To0 implies To0 00 e 52ds T3 826 0 Note that the error function is de ned as 2 m erfz WA e s2ds7 827 which has properties erf0 0 and erfoo 1 Then7 co 2Too 7 TB so that T07 TB 2 Too 7 TB Aquot 5 9st TB Too 7 mam 828 W This solution is plotted in Figure 4 We see that T is a monotonic function of z and that there is a nonuniforrnity at z 0 for t near zero Here the isotherms lines of constant temperature are horizontal Suppose we wish to know the location7 27 of a speci c isotherrn T T where To0 lt T lt TB Then we use the equation T TB Too 7 TBerf777 829 to determine 77 and then use 77 z2 to nd that the location of the isotherm T T is given by 2 MM 830 Since 77 is constant any particular isotherrn advances like 852 Similarity Solution 1D Heat Transfer Localized Heat Source For a one dimensional heat transfer problem ut Hum with localized source of initial heat E5 per unit area at z zo and boundary conditions u a 0 as a 007 Barenblattz pp 64776 identi es the similarity solution Q 7 7 2 4 0z t e m 70gt 7 831 7 F 7 2 met where Q ESpCp Note 30 de Q EXERCISE 81 Derive equation 831 The stability of such a similarity solution is discussed in Barenblattz7 pp 2097219 August 21 2008 Industrial Math Notes DM Anderson Temperature Temperature Pro le for 1 D Heat Equation r r r r r Figure 4 This gure shows the error function solution to the one dimensional heat equation Here we have taken TB 1 and T00 0 Values of at are 00017 0017 017 027 05 and 10 Note that there is a nonuniformity at z 0 for t near zero 65 August 21 2008 Industrial Math Notes DM Anderson 66 853 Green7s Function solution for the 1D Heat Equation Here we describe the Greens Function solution for the 1D heat equation in an unbounded spatial domain Consider the dimensionless heat conduction problem ut um q7t for lt 007 t gt 07 M967 0 f967 832 for u ux7t Let Lu 7U um and note that the adjoint operator is Du U um Note that for an appropriately de ned inner product 739 A lt 7 gt dzdt 833 0 7A the adjoint operator is de ned by Luw u7 LW Note vLu 7 uLW 07ut um 7 111 um7 7WD Wm 7 Wm 7uvt 7 um 7 yum 834 De ne 6 E 616 8 835 and I3 um 7 uu m 181067 836 so vLu 7 uLW 76 Here 61 and 6 are unit vectors in the z and t directions If we integrate this over 739 A d dt dR 837 0 U gt z Rlt gt lt gt we nd vLu7uL1 dR 7 W dR7 I3 dS 838 R R 373 where 73 is the outward unit normal to 372 Note that R represents the domain 7A lt z lt A and 0 lt t lt 739 So 7 A 7 A 7 P dS P 7 tt0dz P tt7dx7 878 7A A 0 13 zzAdt 0 P awryAda A d A d 7 A 7ltuvgtto z Ammo 8 Tuvm 7 vummAdt T 7uvm 7 vummAdt 839 0 0 August 21 2008 Industrial Math Notes DM Anderson 67 So we have the identity 739 A A 739 vLu 7 uLW drdt 7 1w dx 7 uni 7 vuwfA dt 840 0 7A 7A 0 Wed like to solve 832 using a free space Green7s Function The Greens Function for this problem z7t 7n satis es 1 07 tlt 777 841 1371m 7 07 as 7 00 The associated adjoint Green7s Function ltIgtx7t 777 satis es ltIgt 07 tgt777 842 ltIgt7ltIgt 7 07 as M700 If in our identity 840 we have 739 gt 7771i ltIgtz7t 777 and i ltIgt7t 777 we obtain 739 A 0 A ltIgtltztigxngt6ltz 7 6W 7 n 7 ltIgtltztisngt6ltz 7 we 7 7m dzdt 7 o 843 if we take A 7 00 Note that ltIgt 0 if t 739 gt 77 and ltIgt 0 at t 0 lt 77 Therefore n n 7 lt1gt 777 777 0 or 18067106777 lt1gt 777z7t 844 Next in our identity 840 we have 739 gt 777 u u7t our solution to equation 832 and i ltIgtx7t 777 In this case the delta function res in the domain of integration Then OT mm ngti7qltz m 7 we tgtlt76ltz 7 6W 7 77 dzdt 7 mA 1 ux7tltIgtx7t 777 ilkOT m 4mm x H 7A dt 845 0forA7oo Note that ltIgt 0 for t gt 77 and so the upper limit in the time integral 739 can be replaced by 77 Also7 we consider the limit A 7 00 This equation reduces to n A A 7 0 AltIgtltzt5ngtqltztgtdzdtman 7 fltzgtlt1gtltzoisngtdz 846 7A OF 14677 gig A 7 A Aflt1gt70 777dz0 LAqz7tltIgtx7t 7ndzdt 847 August 21 2008 Industrial Math Notes DM Anderson 68 Renarning the variables A t A um 3930 Aims 70l967td 0 Aqlt 7ngtlt1gt nln d dnl 848 Finally7 noting that ltIgtx7tl 777 mlnt we have u7 t A t A 3330 A f lt1gt967tl 70d 0 A qlt 7 77lt1gt967 tlmd dnl 849 This is a solution to problem 832 involving the data f and q integrated over the appropriate domains However7 to use this solution we need the free space Green7s Function ltIgtz7 tl n Recall that z7tl 7n satis es equation 841 First consider 7 m a 0 as a 007 850 b 0 for t lt 0 Note that the last condition here says that nothing happens until the 6 function res We solve this problem using a Fourier Transform where 43 Ooe m 7td7 851 and um Leiwmw7tdw 852 Multiply equation 850 by 6 and integrate A e iwa tdz e iwa mdx 5z5tdz7 it 57iwm m i000 7 m7 weiiwmd 6097 0 OO m 455W 7 7 we iwdz my 853 0 so it 7w2 3 6t 854 This equation has the general solution A Ae w2t t lt 0 August 21 2008 Industrial Math Notes DM Anderson 69 where A and B are to be determined First note that b 0 for t lt 0 so A 0 Next if we integrate equation 854 across t 0 we have e tdt wzje dt76tdt 436 wZla dt 856 The integral term on the right hand side vanishes in the limit 6 a 0 since i5 is locally integrable Therefore 0 7 150 1 This irnplies 7W2 7 33 Be i 0 1 857 so B 1 so A 7 0 t lt 0 iwzt 7 gt 0 mm 858 where Ht is the Heaviside function H 0 if t lt 0 H 1 if t gt 0 Therefore 1 00 139qu 7w2t Mat e e dw for t gt 0 859 7T 700 The following calculation shows that we can rst convert it to an ODE and then solve the ODE Note 69 i 00 2 7 welmmeiw tdW7 a i ioo Zr 00 mm 3 7W2 1 iooe Elte tgtmdw 1 0 6 2t 77 mm 7 7w d 47ft 700 ar an lt6 gt wl u T I woo 0 1 7i elwzeiwzt 7 Zieuumeiu tdw 77 2t myoo 700 2w 0 1 0 2t 77 7 1 7 d 860 2t 2w co 6 6 w l as So b satis es the differential equation da z a was 861 August 21 2008 Industrial Math Notes DM Anderson 70 where t appears only as a parameter This ODE is subject to the initial condition 1 0 2 7139 0t 7 d 7 862 66 gt 2Wme wt lt gt W 2 17rt Equation 861 can be solved by integrating do i pdz 6 7 2t 7 2 ln 7Ec7 863 or 68 0562640 864 Subjecting this to 862 leads to the result 1 x t e m24tl7 865 7 2x7Tt for t gt 0 This is the free space Green7s function for the one dimensional heat equation with a source at the origin Now we really want the solution to the problem in equation 841 Note that the solution we seek is just a shifted version of the solution 5 that we have obtained in 865 Therefore7 1 m 68 n eltw gt2lt4lttquotgtgtHltt e n 866 2 7T6 7 77 This is the function we use in 00 t 00 MM own 0gtd 0 qlt ngtlt1gtltzt ngtd dn 867 which solves our original equation 832 Compare this to the solution to the localized heat source solution 831 of the previous section ie here use qp7t 0 and fp 6z0 What is the long time behavior of this solution for general fx compared to the localized heat source problern EXAMPLE 81 Suppose we want to solve ut um7 t gt 07 lt 007 868a up7 0 6p 7 0 868b The solution to this represents the thermal eld associated with a localized initial disturbance Our solution 867 takes the form 1 mar4615 1 00 u z t 6 7 z 7 700 g 027Tt 2x7Tt A sketch of this function at di erent ualues oft shows the slumping of a delta function initial heat source We have an initial spike which attens out as t increases e 024t7 tgt0 869 August 21 2008 Industrial Math Notes DM Anderson 71 EXAMPLE 82 Suppose we want to solve ut um7 t gt 07 lt 007 870a u7 0 870b Here our solution 867 takes the form 00 1 2 1 0 2 7t H ems416 7 W 4td 871 M96 700 92 7f 6 2 Tm 7006 6 If we let r x 7 so dr id xZE and note that g 0 corresponds to r and 700 corresponds to r 00 Then 1 m B 2 7 7 4t r d 2m f 00 5 T 1 0 2 7 r d W A e r7 gamutm 872 u7t where the error function and complementary error functions are de ned in A5 This so lution has the properties I as x a 00 uz7t a 0 as x a 700 u7t a 1 at z 0 u07t 12 4 for pedp as t a 00 u7t a 12 and 5 for cced x as t a 0 u a 0 ifp gt 0 andu a 1 ifp lt 0 If we sketch this solution initially it is a down step and as time increases the step di uses or smooths out 86 Homogenization Methods See for example7 Hinch7 section 74 Here we use ideas of multi scale methods to derive an approximate model from exact models such as the heat conduction model we7ve just studied In the heat conduction problems so far we have taken h and p017 constant so a hpCp is constant Homogenization methods can be applied to situations in which these coef cients may vary rapidly in space or time 7 for example7 a rapid oscillation about a slower77 variation of the coef cient One can imagine the situation of thin alternating layers whose thermal properties differ in each layer To describe the method7 we shall consider heat conduction in one dimension in which the thermal conductivity h varies rapidly in space so that we may write h where E is a small parameter characterizing the length scale of the variations An example is h 1 Acos27Tze which has rapidly oscillating thermal conductivity about mean value 1 So suppose h and for now take p017 constant We shall assume that is a periodic function Hinch describes an approach for the more general case The heat equation in this case is pCpg g 873 August 21 2008 Industrial Math Notes DM Anderson 72 The dif culty with solving this equation lies in the non constant coef cient k Our object in this case will be to use homogenization techniques to derive an effective7 heat equation in which a constant7 effective thermal conductivity replaces the variable one Such an effective7 heat equation will not capture the details of the thermal conductivity structure but will capture the more global features of the problem Homogenization techniques are a type of multi scale method In this case we de ne a new variable Z 26 and treat Z and 2 as independent variables As in multi scale methods7 the idea behind this is that Z and z are associated with vastly different length scales as E a 0 Here we focus on just the heat equation for now ignore any boundary or initial condi tions We introduce the transformation 6 6 16 a a EE7 874 so the governing equation becomes 6T 9 1 9 9T 13T pOpE 875 Now we expand the temperature in an asymptotic series T N T027Z7t 6T127Z7t 876 At 06 2 we have 3 3T0 0 7 E 877 This equation can be integrated once to nd kZ6T03Z 1027 t where 10 is an arbitrary function of z and t lntegrating again gives T027Z7t a0z7tOZalz7t 878 Now in order for To to be a bounded function of its arguments 7 in particular Z 7 we require that 1027 t 0 Therefore the leading order temperature is independent of the fast variable Z AtOe 1 we nd i 3 6T0 6T1 a aT39O 0 7 E a 7 880 August 21 2008 Industrial Math Notes DM Anderson 73 6T0 6T1 k Z 7 7 62 62gt where be is independent of Z Solving for dTlQZ gives 6T1 b0lt27t 7 6T0 b0zt 881 7 7 7 882 92 kltZgt 62 Now if we integrate this expression over one period we suppose the period of the function kZ is Zp and that the functions T07 T17 etc are also periodic functions with the same period 7 we can examine the functions we identify for these quantities to assure that they indeed have these properties so that Zp dZ 9T0 TZZ 7T Z0 b 77 7 883 1lt pgt1lt gt 00 Z W ltgt zero by periodicity Consequently we nd that 6T0 b t k l 17 884 w gt lt gt 62 lt gt where we have introduced the notation ltgt 7 1 Zde 885 i Zp 0 39 39 We can go back and solve for T1 to nd that 6T0 Z ltk71gt71 Tb t 7 71dZ 886 1 127 62 0 7 where b1zt is an arbitrary constant One can con rm that this function is periodic with period Zp At 01 we have 9T0 6 9T0 9T1 6 9T1 9T2 If we use equation 881 we have 6T0 abo i a 6T1 6T2 independent of Z We now integrate over one period and nd that the resulting expression on the right hand side vanishes and we are left with the equation 6T0 POPE ltk71gt 1 62TO 322 39 889 August 21 2008 Industrial Math Notes DM Anderson 74 This is an effective7 heat equation in which the constant effective thermal conductivity is the harmonic average of the actual thermal conductivity function In this case we can de ne an effective thermal diffusivity by 1 K8 7 1 17 890 ff poplt gt where 1 Zp dZ k 1 17 7 891 lt gt ZP 0 Z lt gt Recall our solution to the one dimensional heat conduction problem using the similarity solution7 T TB Too 7 TBerf77 where 77 892 2 2V Kt Suppose now that we consider a similar one dimensional problem in which the thermal conductivity has the form kZ koll bsin27rZ7 where Z 26 Based on the above homogenization calculation we can use equation 889 to write down a homogenized solution to this layered7 problem To TB Too 7 TBerf775ff where neff 893 2 WW7 where Heff is given by equation 890 That is7 we can apply our previous results for uniform thermal properties if we just use the effective thermal diffusivity rather than the original constant value Again we emphasize in this particular case that it is the harmonic mean rather than perhaps the more familiar average7 of the arithmetic average The table shows values of the effective thermal diffusivity as a function of b Note that the arithmetic mean of the thermal diffusivity function is no while the harmonic mean of the thermal diffusivity which in this case is the effective thermal diffusivity can be signi cantly smaller than no especially as b approaches 1 so that the minimum thermal diffusivity value approaches zero This behavior can be compared to a number of other physical systems including transport in layered porous media7 series and parallel circuits7 etc 9 Solidi cation of 21 Pure Material These notes are based on lecture notes given at Northwestern University by Davis and Worster Much of these notes have since appeared in two publications Worster7 2000 and Davis7 2001 Refer to the book by Davis 2001 for topics on solidi cation of a pure material as well as binary alloys Another good reference is the chapter by Worster in 7Perspectives7 2000 August 21 2008 Industrial Math Notes DM Anderson 75 1 b 1 Kmam 0 1 Kmin 0 1 Keff 0 1 01 11 09 0995 05 15 05 0866 09 19 01 0436 0999 1999 0001 00447 Table 1 This table shows the maximum and minimum values of nno 1b sin27rZ where no kOpCp and the effective thermal diffusivity the harmonic average as a function of b Note that the arithmetic average of this function is no for all values of b 91 Basic Concepts and De nitions DEFINITION 91 Melting Point The melting point ofa pure substance is the temperature TM below which the substance is solid and above which the substance is liquid once it is in a completely stable thermodynamic equilibrium In general7 the melting point depends on the pressure P so that TM TMP this is the Clausiusiclapeyron relation but for the present purposes we shall take TM to be constant DEFINITION 92 Latent Heat The amount of heat energy required to change a unit mass of solid into liquid at T TM We denote the latent heat by L and note that is can be expressed in units of calg ln solidi cation problems7 we generally have both liquid and solid phases present at the same time Additionally7 the position of the interface separating the solid and liquid phases is not typically known in advance it must be determined as part of the solution to the problem Consequently7 solidi cation problems typically fall into the category of free boundary problems For instance7 in the region occupied by the solid phase7 the temperature satis es 6T 95017567 V lkSVTsh 91 and in the region occupied by the liquid phase7 the temperature satis es 6T 0ij V lkLVTLl7 92 where p57 CPS and k5 are the density7 speci c heat and thermal conductivity of the solid phase and pL7 CPL and kL are the density7 speci c heat and thermal conductivity of the liquid phase As was the case for heat conduction7 these equations are subject to boundary conditions and initial conditions Boundary conditions include those at the boundaries ofthe container August 21 2008 Industrial Math Notes DM Anderson 76 holding the substance eg a constant temperature condition or a zero heat ux condition but also would include boundary conditions at the interface between the solid and liquid phases1 911 Interfacial Conditions In the simplest case7 one boundary condition applied at the interface between the solid and liquid phases is that the temperature of this interface remains at the equilibrium melting temperature T TM397 that is TL TM and TS at the solid liquid interface We shall see later that there are important modi cations of this condition due to curvature of the interface7 interface attachment kinetics and other effects The consideration of the conservation of heat at a moving interface between solid and liquid phases provides another condition The basic statement of this condition is that the amount of latent heat generated conductive heat ux away from the interface 94 That is7 if we have a solid liquid interface moving with normal velocity Vn where 7 is a unit normal vector of the interface pointing into the liquid the amount of solid generated per unit area in time St is Vn t 95 This requires the removal of an amount of heat per unit area ngVn t 96 This must be balanced by the conductive heat ux away from the interface it 72 7 as m at 97 The above balance statement then requires that pSLVn it n 7 qs 71 98 This says that the jump in heat ux across the interface is balanced by the generation of latent heat due to phase transformation note that if Vn 0 the interface does not move and there is no phase transformation If we assume that the heat uxes are given by f5 7kSVT5 and j L ikLVTL then this condition becomes psLVn 5st 73 7 kLVTL 1The boundary condition 7count7 for such problems will be made more explicit when we consider specific solutions in simple geometries August 21 2008 Industrial Math Notes DM Anderson 77 Note that the interface velocity can be written in terms of a time derivative of the position of the interface see details below Technically7 there are additional modi cations of this condition for a curved interface but these effects are typically neglected Equations 91 and 927 interfacial conditions 93 and 99 together with appropriate boundary conditions at domain boundaries and initial conditions constitutes a full system of equations to solve for the unknown solid and liquid temperatures which are both functions of space and time7 in general and the unknown solid liquid interface position which is a function of time7 in general In the sections below we consider some simple geometries in which we work out these details 92 1D Freezing From a Cold Boundary Consider a semi in nite region 2 gt 0 containing a uid at a uniform temperature TL TM bounded by a solid wall at z 0 Suppose for t 2 0 the temperature of the wall is held at a xed cold temperature of TB so that a layer of solid phase begins to solidify at this boundary Below we write down the governing equations for this solidi cation problem Let 2 at denote the position of the solid liquid interface at time t The equation governing the temperature in the solid region 0 lt z lt at is 6T5 62T5 7 910 at S 622 l The equation governing the temperature in the liquid region 2 gt at is dTL 62TL at KL 622 911 The boundary and interfacial conditions are T5 TB7 at z 07 912a TS TL TM397 at Z 1lt7 da 6T5 dTL L7 k 7 7 k 7 t t 912 PS dt 5 62 L 627 a Z 007 0 TL TM7 as 2 a 00 912d Note that in this geometry Vn dadt lnitially7 there is no solid phase so 10 0 and the liquid is at uniform temperature TM We shall make the simplifying assumption here that the temperature in the liquid region remains uniformly at TL TM Note that TL TM satis es the governing equation in the liquid region as well as the temperature boundary condition TL TM at z 0 and z a 00 Before solving this problem we shall write it in dimensionless form Let TL TB TM TB6L7 Ts TB TM T191957 t Tlt7 z le7 a Zla 913 August 21 2008 Industrial Math Notes DM Anderson 78 where 1 and 65 are dimensionless temperatures in the liquid and solid phases t and 2 are dimensionless time and space variables and a is the dimensionless interface position The quantities 739 and Z are time and length scales to be determined Our problem becomes 6L 1 in 2 gt a t 914a 66 H 739 626 7 a 325 in 0 lt 2 lt a t 914b with boundary conditions 65 0 at 2 0 915a 65 1 at 2 a t 915b da kg 7 TB 605 w at HW We Now choose Z sh Then k57TM 7 TB 7 1 7 i 916 pSLZ2 L S where AT TM 7 TB is known as the undercoolmg and S is the Stefan number The Stefan number measures the latent heat relative to the heat content of the material Dropping primes we nd that our problem simpli es to 1 1 and 665 6265 7 E 622 1n 0 lt 2 lt at 917 subject to conditions 65 0 at 2 0 918a 65 1 at 2 at 918b d 619 SET 67 at 2 at 918c As was the case in section 851 where one dimensional heat conduction from a boundary was considered we observe here that there is no external length scale in the problem This suggests an approach where we seek a solution that depends only on the combination Therefore we introduce a similarity variable i 2 77 i The factor of 2 is chosen for convenience later in the calculation It would not be wrong to introduce 77 or 77 2212 however it will turn out that the above choice is the most 919 natural one As before we nd 6 1 d 62 1 d2 6 77 d 777 777 7777 920 62 2xfd777 622 4td7727 6t 2td777 August 217 2008 Industrial Math Notes DM Anderson 79 so that the equation in the solid region becomes d g 2 d265 7n f quot7772 Note that if we substitute these transformations into the heat ux boundary condition we nd that 0 921 Sda i 1 deg dt 7 2xEd777 This suggest that at has the form at n at2f 922 at 2AxE 923 where A is a constant to be determined Therefore7 we must solve dds d205 a 277 dnz 07 on 0 lt 77 lt A7 924a 05 07 at 77 07 924b 05 17 at 77 A7 924c dds 7 2AS t A 924d dn 7 a 77 Note that we have a second order ordinary differential equation with two boundary conditions and an additional condition that will determine the position of the interface ie determine the value of A The ODE has solution erf77 6 7 925 577 erfA 7 which satis es the rst two boundary conditions Noting that this solution has d6 2 1 57 A 775 926 E end we nd that the third boundary condition requires that S 1 WA exp2erf7 927 which we interpret as an equation for the unknown A Analysis of this equation reveals that A is a monotone increasing function of S l7 and in particular7 there is a solution for all S gt 0 Recognizing that the value of S can be controlled by the undercooling or the wall ternperature7 we see that S increases with decreasing undercooling So7 since at 2Axf7 increasing undercooling rneans decreasing S and increasing A that is7 growth is faster with larger undercooling Conversely7 if one were able to control the latent heat7 larger latent heat increases S and hence decreases A large latent heat means that a lot of heat must be conducted away from the interface in order to advance the interface 7 slower growth As was recognized in Section 851 the similarity solution fails at t 0 However7 we expect that the similarity solution here is approached by any solution of the initial value problem August 21 2008 Industrial Math Notes DM Anderson 80 93 1D Freezing Into an Undercooled Melt In this situation7 the transformation from liquid to solid is accomplished by 7undercooling7 the liquid phase to a temperature To0 lt TM In practice undercooling often occurs and solidi cation of the liquid phase begins by a nucleation event The geometry that we consider here is the same as that for one dimensional freezing from a cold boundary That is7 the solid phase is located in a region 0 lt z lt at and the liquid phase is located in a region 2 gt at In the present case7 we now assume that the solid phase is isothermal with T5 TM and that the far eld temperature in the liquid phase T00 is less than TM Taking into account an isothermal solid region7 our mathematical model of this problem is then dTL 62TL W KL 622 7 in the liquid region 2 lt 1t7 928a TL a 007 as 2 a 007 928b TL TM7 at z 1t7 928c da dTL ngE i ikLE at z i 1t7 928d We again seek a similarity solution where 2 7 t 2AM t T t TOO AT6 929 77 2M7 a KL7 L27 777 where AT TM 7 To0 gt 0 and 6 satis es d26 dinz i Qng 07 907nm 07 977A 17 9300 2A8 7777 A As before7 S LCPAT Solving the equation for 077 and applying 0A 1 and 600 0 gives 977 931 Enforcing the heat ux boundary condition gives 8 1 WA expA2erfcA 932 We interpret equation 932 as giving the factor A AS in the equation for interface position August 21 2008 Industrial Math Notes DM Anderson 81 Investigation of this result shows that A increases monotonically with S 1 for S lt 1 However one also observes that A 7 00 as S 7 1 There is no similarity solution of this form ie equation 932 has no solution when S S 1 An observation we can take from this is that the interface is predicted to move in nitely fast as S approaches 1 from above This nonphysical prediction suggests the need for more physics in the model One such piece of physics that is important at high interface speeds is called 7attachment kinetics7 931 Attachment Kinetics Effects The following notes are from Worster7s Solidi cation notes Solidi cation is a reaction whose rate depends on the level of disequilibrium undercooling in this case driving it 0 When the interface temperature TI is equal to the equilibrium melting tem perature TM a solid and liquid interface is in dynamic equilibrium with atoms continually attaching and detaching from the interface at equal mean rates 0 When TI lt TM atoms prefer to attach to the interface rather than detach lf TM 7 TI is not too large the interface advances at a rate V that increases with TM 7 TI 0 As TI decreases further the attachment process slows down atoms become 77sluggish but not quite as slowly as the detachment rates slow down 7 giving solidi cation at an overall slower rate We7ll use a linear law 1 V TM 7 T1 933 where u is a kinetics parameter which we take to be constant This equation can be inter preted as giving the interface temperature T1 in terms of the interface speed TI TM 711V That is indicates that the interface temperature may deviate from the equilibrium melt ing temperature due to attachment kinetic effects This is also sometimes called kinetic undercooling With this new effect of kinetic undercooling we shall revisit the problem of one dimensional freezing into an undercooled melt With the same geometry as before we now have dTL 62TL W W7 TL 7 00 as 2 7 oo 935b in the liquid region 2 lt at 935a August 21 2008 Industrial Math Notes DM Anderson 82 TL TI7 at z 1t7 935c da 6TL nga 7 7amp3 at z 7 1t7 935d TM 7 T1 at z at 935e The added feature here is that the interface temperature is now directly coupled to the interface dynamics Note that the interface temperature is now part of the solution rather than an input to the problem We shall nondimensionalize this problem by introducing t Tlt7 z Xlz7 a Xla7 T To0 TM 7 Too6 936 These give dropping primes 2 g g in the liquid region 2 lt 1t7 937a 6 7 07 as 2 7 007 937b 6 617 at z 1t7 937c p L X da k T 7 T00 36 Ska la 7 L llX gt57 at Z W7 93 Wu 7 91 at z at 93 where 61 TI 7 TogTM 7 Too From the PDE7 we choose XV739 KL From the kinetic boundary condition we choose XlT TM 7 TogM7 so that 2 L KLM X 7T TW 938 l l M M m TM 7 M lt gt The governing equations then are 2 it g in the liquid region 2 lt 1t7 939a 6 7 07 as 2 7 007 939b 6 617 at z 1t7 939c da 66 S7 77 t t 939d dt 62 a z alt gt lt gt d d 17 617 at z at 939e Notice that unlike the problem of freezing into an undercooled melt in the absence of kinetic effects where we only obtained a single relationship between length and time scales suggest ing a similarity solution approach7 here our nondimensionalization has given speci c length and time scales related to physical parameters in the problem August 21 2008 Industrial Math Notes DM Anderson 83 Let us seek a constant speed front solution to these governing equations That is seek a solution of the form 9 92 4 Vt at Vt 940 where V dadt is a constant to be determined With the above form for 0 if we let C27Vtwe have 2 2 g g E a 941 at dC 622 dCZ These reduce the governing equations to 2 7Vj7 gig in the liquid region C gt 0 942a l9 7 0 as C 7 oo 942b l9 01 at C 0 942c 0 SV 77 t 0 942d 6C a c lt gt V 17 01 at C 0 942e The above ODE and boundary conditions 190 01 and 600 0 imply that 90 015 943 Two unknowns 01 and V remain to be determined Substituting this result into the heat ux boundary condition yields d6 dC so that 01 S Then the remaining interface condition requires V 1 7 S This solution is valid for S lt 1 so that V gt 0 and the liquid solidi es Summarizing our results shows that for 8 1 lt 1 a similarity solution is possible with 5V C 0 V91 944 at For 8 1 gt 1 a constant speed solution is possible where attachment kinetic effects are important For the case S 1 Umantsev23 where at tZB 94 Morphological Instability Observations show that planar solid liquid interfaces can become morphologically unstable which means unstable to shape perturbations for example see Kurz and Fisher From a practical point of view morphological instabilities can in uence the properties of the solidi ed material 7 this is especially true for alloy solidi cation in which diffusion of chemical species becomes important 7 so it is valuable to try to understand the mechanisms for morphological instability August 21 2008 Industrial Math Notes DM Anderson 84 A Mechanism for Morphological Instability Consider the solidi cation of a planar inter face into an undercooled melt for the case of an isothermal solid phase Here the driving force for growth is the temperature gradient dTdz in one dimension in the liquid phase The temperature near the interface is relatively 7hot7 while the temperature in the far eld is relatively 7cold7 undercooled The ow of heat is away from the interface Now suppose that we 7perturb7 the interface by making the solid bulge out slightly into the liquid region Further suppose that the interface remains isothermal and that the tem perature in the far eld is unchanged In this case the isotherms in front of the bulge become compressed and the magnitude of the temperature gradient is increased There fore the heat ow away from the interface is increased and consequently the local growth rate is increased Hence the perturbed 7bump7 tends to grow faster and become a bigger bump Thus undercooling is a destabilizing mechanism This scenario despite its simpli cations explains one possible means of amplifying an interface perturbation and causing a morphological instability Note that we have previously argued that the interface temperature may in fact depend on the normal velocity of the interface and so our assumption that the interface remain isother mal may not be correct Also note that we have not incorporated the effects of interface curvature and interfacial energy into our formulation Before addressing the morphological instability problem in more detail we next turn to a description of curved interfaces so that those effects as well as ones previously discussed can be incorporated into a solidi cation model 941 Curved Interfaces The effect of curvature and surface energy on the solid liquid interface temperature is gener ally known as the Gibbs Thomson effect The following conceptual argument gives a general understanding of this effect Consider an interface between solid and liquid phases that instead of being planar varies sinusoidally thus making 7hills7 and 7valleys7 on the solidiliquid interface A solid atom at the top of a hill has relatively fewer 7neighbors7 than does a solid atom at the bottom of a valley Where atoms have fewer neighbors they are 7less bound7 to the interface and accordingly the interface temperature will be lower there one needs to cool the atom to a greater degree in order to get it to attach to the solid On the other hand where solid atoms have more neighbors they are 7more bound7 to the interface and the interface temperature is correspondingly higher there one needs to cool the atom to a lesser degree in order to get it to attach to the solid Therefore since the 7hills7 need to be colder to solidify while the 7valleys7 do not need to be as cold to solidify the tendency is for the hills to melt back and the valleys to freeze up Thus the Gibbs Thomson effect is stabilizing More formally we write the Gibbs Thomson condition on the interface temperature as T1 TM i W 945 August 21 2008 Industrial Math Notes DM Anderson 85 where TM is the equilibrium melting temperature7 P TM yng where y is the surface energy per unit area units of calcm27 L is the latent heat units calg and p5 is the solid density units gcm37 and is a unit normal vector to the interface pointing towards the liquid phase Here7 V is twice the mean curvature of the interface Note here that for the purposes of describing the Gibbs Thomson undercooling effect we are neglecting kinetic undercooling In two dimensions with the interface parameterized by z hz7 t the unit normal is ihm 1 7 7 7 946 1 h so that 6 h h v A 7 7m 7 947 ax MT15gt wag2 gt To simplify this slightly7 assume that the interface de ection is small so that hi ltlt 1 Then TI TM th 948 From this form it can be seen that the interface temperature TI is less than TM where the interface is concave down7 or hm lt 07 ie a 7hill7 while TI gt TM where the interface is concave up7 or hm gt 07 ie a 7valley Again we see that 7hills7 melt back and 7valleys7 freeze up the Gibbs Thomson effect is stabilizing The curvature and surface energy effect also modi es the heat ux boundary condition and a correction to our original form can be written as psLVn lt1 i psiLv ksVTS i IWTL 949 where Vn is the normal velocity of the interface For references see Wollkind and Maurerz6 or Umantsev and Davis Finally7 note that our previous attachment kinetics effect led us to T1 TM 7 Therefore7 points on the interface moving faster ie a forming 7hill7 have correspondingly lower interface temperature than do points on the interface moving slower ie a 7valley Therefore7 kinetic undercooling is a stabilizing effect 7 the hills tend to melt back and the valleys tend to freeze up Putting the effects of attachment kinetics together with the Gibbs Thomson effect gives us the following expression for the interface temperature TI TM i W a i W 951 August 21 2008 Industrial Math Notes DM Anderson 86 Summarizing our situation for solidi cation into an undercooled melt we have a driving force due to thermal undercooling which is destabilizing However we have now the effects of kinetic undercooling which is stabilizing and the Gibbs Thomson effect which is also stabilizing That is we have competing mechanisms that combine in some way to in uence the ultimate morphological stability or instability of the solid liquid interface Which effect wins In order to examine this question we shall use a linear stability theory 95 Morphological Instability of a Planar Interface Growing Into an Undercooled Melt See also Davis6 pp 35741 Let us revisit again the problem of solidi cation into an undercooled melt The solid liquid interface is de ned by z a t where we note speci cally that we allow the interface shape to be non planar De ne a unit normal vector to the interface f1 which points into the a A liquid phase Denote the normal velocity of the solid liquid interface by Vn where Vn Vn We de ne the interface velocity vector by T7 dadtg where I is a unit vector in the z direction For simplicity we shall assume that the solid is isothermal The governing equations for this problem are 66 HLVZTL for z gt ax t 952a TL 7 Too as 2 7 oo 952b TL TM 7 TV 7 Ma at z ax t 952c ngVn 7kLVTL at z at 952d We shall take the approach in which we identify a simple 7base state7 solution correspond ing to a planar solid liquid interface Then we shall 7perturb7 this solution and examine how these small perturbations grow or decay in time As we did before we shall nondimensionalize this problem by writing t Tt z X 2 a Xa T To0 ATQ 953 where AT TM 7 Too Also as before we choose X 2 X AT L l Q 7 954 Tl Tl M based on balacing terms in the PDE and also balancing the temperature and kinetics terms in the interface temperature condition Further we have A 7 70171 n 7 1a 12 955 August 21 2008 Industrial Math Notes DM Anderson 87 Vn 7W7 956 i A A at v 7 atkni 1ag127 957 A am 1 Vn 7 6m 7 1a 12gt01 lt lt1ag12gt7 958 where subscript t7 z and 2 indicate partial differentiation with respect to that variable Our resulting dimensionless system is dropping primes 66 626 626 a y 7 for z gt 1z7t7 959a 6 a 07 as 2 a 007 959b L L 7 6 7 1 P0lt1 lag2 1 air27 at z i 1z7t7 959c 1 at 7E 61 7 116m7 at z 1x7 t7 959d where dimensionless parameters To and S are de ned by P L P 7 S 7 960 0 ile CPAT 951 Base State Solution For the case of a planar solid liquid interface so 1amp7 t at and where 6 is independent of the horizontal coordinate 7 equations 959 reduce to equations 939 given in the section on interface attachment kinetics Therefore7 a base state solution to equations 959 with a planar solid liquid interface corresponds exactly to the solution identi ed earlier Denoting the base state with subscript 7B7 63 Sexp7VC7 961a aB V22 961b V 1787 9610 where C z 7 Vt Since this solution corresponds to an interface moving at constant speed7 we can view it as a solution that is steady in a frame that moves with the interface velocity Therefore7 it will be convenient for the stability analysis that follows to introduce a change of variables corresponding to a change of reference frame Speci cally7 we introduce z L 962 E t7 963 c z 7 Vt 964 August 21 2008 Industrial Math Notes DM Anderson 88 so that 6 6 6i 6 6t 6 6C 7 6 a n 57E37 93965 6 6 6i 6 6t 6 6C 6 7 n if if 7777 966 62 6s 62 6t 62 6C 62 6C 6 6 6i 6 6t 6 6C 7 6 6 a n EEEampE E V C 93967 Additionally7 we de ne hx7t a7t 7 Vt so that hm am7 ht at 7 V and in the base state hB 0 This gives us dropping bars 66 66 626 626 7 7 V7 7 7 f h t 968 at a g 627 or C gt 7 7 a 6 7 07 as C 7 007 968b hmm ht i V i 6 7 1 Tom 7 W7 at C 7 hz7t7 968c ht V 7 7 94 7 1191 at g 7 W707 968d 952 Linear Stability Analysis The linear stability analsis proceeds as follows Recall that the base state solution has 63 Sexp7VC7 hB 0 and V 1 7 S We introduce perturbations amt 93lt0 zltt 969a Mm hBh zt 96 where 6 and h are perturbations of the base state solution We shall be interested in only in nitesirnally srnall perturbations and so we substitute the above expressions for 6 and h into equations 968 and linearize in the prirned quantities Later7 we shall consider speci c forms for the perturbations 6 and h but for now we just linearize the equations Let us examine rst the PDE Substituting 6 63 6 into equation 968a gives 66 d6B 66 626 d26B 626 7 7 7 7 2 970 6t dc 6 6x dc 6C2 However7 by de nition of the base state we know that d6B d26B 7 TC 7 CK 7 971 so that our linearized PDE becomes 66 66 626 626 7 972 E7 a egg2ec239 August 21 2008 Industrial Math Notes DM Anderson 89 Note that this equation contains only terms that are linear in primed quantities Note also that the original PDE was itself linear Next7 the far eld boundary condition 6 7 0 becomes 03 6 7 0 Since we know the base state solution has the property that 63 7 0 then we must have 6 7 0 Now consider the interfacial temperature condition 968c Substituting into this equa tion gives h h V 636 1P0 m 7 1 h 32 1 19127 at C 95713 973 Unlike the previous two equations7 this is a highly nonlinear equation evaluated at C h and the process of linearizing it is considerably more involved First7 we note that the quadratic terms h will be neglected This leads to the result 63 6 1 Fobm 7 h 7 V7 at C h x7t 974 Recall that the base state solution satis es 63 17 V7 at C 0 975 It is tempting to therefore cancel these same terms appearing in equation 974 but this would not be correct Note that equation 975 says that 03 C 0 1 7 V However7 it does not say that 03C h 1 7 V and it is 03C h 7 and not 03C 0 that appears in equation 974 The problem is that 03C h is not yet linearized about the base state solution Nor is 0 C h ll However7 this can be easily accomplished through the use of a Taylor Series expansions7 63Ch 93ltohlt0 976 and 6 Ch 9 lt0hquot7 lt07 9 lt0 977 Note that we only retain terms that are linear in primed quantities Therefore7 we can rewrite equation 974 as 93K 0 h 0 97c 0 1 Poh m i h i V 978 Now7 using equation 975 we obtain the result h dTK 0 97c 0 roh m i h 979 August 21 2008 Industrial Math Notes DM Anderson 90 Again we note that each of these terms remaining is linear in a primed quantity Finally we consider the heat balance equation 968d The procedure to linearize this equation is similar to that for the previous interface condition Noting that the corresponding base state interfacial condition is 1 dQB V 7 igidC 7 0 980 and using Taylor Expansion on the appropriate quantities we nd that the linearized heat ux boundary condition is dZQB 66 0 7 dcz C M iSht h C0 981 Again we emphasize that the resulting linearized equation should contain only terms that are linear in the primed quantities Therefore our disturbance equations are 96 919 920 920 7 7 7 982 at a 6952 6 a 0 a 0 as C a 00 982b h 0 Fobm 7 h on C 0 982c dZQB 30 Sht dcz 67C on g 0 982d 982e Linear disturbance equations should contain only term that are linear in the primed distur bance quantities It should be clear that 0 0 and h 0 is a solution to these equations Before we proceed we note that since 03 Sexp7V we have dQB TClt0 7V5 983a dflt0 V25 983b We solve the disturbance equations by seeking a 7normal mode7 solution of the form 0z t expUt zaz cc 984a hz t flexpat zaz cc 984b where 039 is the growth rate which may be complex 04 is the horizontal wave number and cc indicates complex conjugate The quantities and Ii are to be determined and we seek 039 0a S P0 to see how different modes grow or decay in time Note that since this is a linear problem we can superpose these modes August 21 2008 Industrial Math Notes DM Anderson 91 Substituting the normal mode solutions 984 into our disturbance equations and sim plifying gives U 7 V0 624 7 0202 a 07 as g a 007 985b fins lt 0 iPoozZi i U137 9850 75w vzs 4lt0 The ODE for has the general solution exp7p We require Rep 2 0 in order for a 0 as C a 00 The ODE itself further requires that 039 Vp p2 7 042 986 This can be solved to give p ltVi 987 but again since Rep 2 0 we have p WNW 988 Therefore7 exp 07 satis es the ODE and the far eld boundary condition The remaining boundary conditions require 7 989a 953 90 13 042D a 7 VS 0 ip 0 Ii Sn 5V2 0 This linear system for and Ii has a nontrivial solution if and only if its deterrninant is zero This requires that SU V2 pa2r0 U 7 VS 0 990 This is called the characteristic equation and it can be interpreted as the condition that gives 039 0047 P07 S Recall that V 17 S In particular7 for given a given set of parameters To and S and for a given wavenurnber 04 this equation determines whether or not that particular mode for those particular parameter values grows Rea gt 0 and is linearly unstable or decays Rea lt 0 and is linearly stable It is informative to examine this result in some special cases First consider the case S 0 so V 1 Here 021 a 07 991 August 21 2008 Industrial Math Notes DM Anderson 92 so that 039 7Gamma0a0 This corresponds to complete stability 0 is always negative when 04 31 0 Second7 the case S 1 so V 0 has 039 0 a2a2F0 039 0 992 It can be argued here as well that 039 lt 0 always when 04 31 0 Third7 a more general examination of equation 990 shows that for P0 0 039 increases monotonically with 042 and asymptotes at a value of 025 For T 31 07 039 is positive for 04 lt ac and 039 is negative for 04 gt do where do is some critical wavenumber This suggests that the presence of P0 stabilizes shorter wavelengths larger values of 04 only longer waves can be unstable Note that long waves have small values of 04 while short waves have large values of a Fourth7 one can examine the long wave limit 042 ltlt 1 by supposing the existence of an expansion of the form 039 00 04202 After some effort7 it can be shown that S 2 Ult 7P0gta7 993 for 042 ltlt 1 This shows that the system is stable to long wavelength perturbations if To gt S1 7S In other words7 the system is stable if the surface energy effect is suf ciently strong Another useful scenario to examine is called neutral stability The objective here is to obtain the boundary in parameterwavenumber space that separates stable ReU lt 0 from unstable ReU gt 0 behavior If we consider only steady modes with 039 real and examine the characteristic equation 990 with 039 0 we obtain 5V2 g V V2 4amp2 0421 0 7 VS 0 994 Noting that 042D S VS7 one can obtain the result that 2 17 SS P017 S 7 1 7 995 a P0 S A plot of 042 versus S shows that an unstable region can be identi ed on the interval PO1 P0 lt S lt 1 Recall that our analysis is not valid for S gt 1 since V lt 0 there as shown in by the base state solution 96 Homogenization 0f the Stefan Problem The one dimensional Stefan problem can be written as 61 7 2 61 9 6t 62 62 7 39 August 21 2008 Industrial Math Notes DM Anderson 93 which applies in both solid and liquid phases Here Is kpcp where H is the thermal diffusivity7 k is the thermal conductivity7 p is the density and Up is the speci c heat The phase boundary7 denoted by 2 ht is assumed to be planar Boundary conditions at the solidiliquid interface include continuity of temperature and a balance of heat TS TLT1I7 LV iimi 99710 where TM is the interface melting temperature7 L is the latent heat7 V is the normal interface velocity and kg and kL are the thermal conductivities of the solid and liquid phases There are additional boundary conditions such as T5 T0 at 2 0 and T a To0 as 2 a 00 We would like to examine how the system can be modeled if the thermal conductivities and the latent heat are assumed to vary rapidly in space so that k5 k5267 kL kL26 and L L26 The heat equation itself has been examined under these conditions by Keller and also HinchM using the multi scale method4 We follow this method here Of particular interest in this case is the treatment of the boundary conditions and how an averaged boundary condition might be derived We rst consider the one dimensional heat equation 667 gltazgtg gt 998 We employ the multiscale method4 where 2 and Z 26 are treated as independent variables so that 6 6 62 6 6 1 6 7 a 7 if 7 if 62 62 62 6Z 62 E 6Z Further we expand the temperature T27 Z7t T027 Z7t ET1Z7Z7 t and collect terms in powers of 6 At 0162 we have 999 i 6T0 i T0 E Hwy 7 0 9100 which gives To T027 t is independent of Z At 016 we have i 6 6T0 6 6T0 T1 7 a K E KZE 7 9101 We can integrate this equation to obtain 92 1 A02 t 9102 August 21 2008 Industrial Math Notes DM Anderson 94 Requiring that T1 be a periodic function we can write upon integrating this expression over a period 1 6T0 A t i 17 9103 ow ltKgt 62 lt gt where denotes the mean over a period In general we can write down a solvability condition for TZ gl en 6T0 1 b M2 T062 7 71 TogldZ 9104 If we assume that the rapid variation is periodic the boundary terms vanish Furthermore7 since To is independent of Z the solvability condition is the the right hand sides gi integrate to zero have zero mean At 01 we have 7 6 3T1 6 3T1 62TO 3T0 7 5T2 5 4 El E l 5 522 t E 927 9105 Note that we can rewrite 92 as 6 3T0 6 6T1 62TO 3T0 77 7 Z 7 A t 7 7 7 7 7 92 62 4 gt62 027 ll 62 gt62 622 6t 6 6 3T1 3T0 Ele lel Kala E7 1 7162TO 6 3T1 3T0 7 7 7 7 7 7 9106 n 622 62 l 62 6t Upon integrating this we obtain the solvability condition 6T0 1 462T 7 7 9107 at ltHgt 322 so that the averaged heat equation involves the inverse of the harmonic mean of the conduc tivity Now suppose we have the system in which in the solid and liquid bulk phases satisfy 6T5 i 6 6T5 W 7 a ltl 5ltZgtEgt 9108a dTL i 6 dTL W i a KLZEgt 9108b subject to boundary conditions at z ht TS TL TM397 6T5 dTL LZV k Z77k Z7 9109b ltgt swag mag lt gt August 21 2008 Industrial Math Notes DM Anderson 95 The leading order averaged bulk equations are QTSO 1 71 627 190 at 622 7 9110 dTLO i 1 7 16sz 3t 7 ltI Lgt 322 9111 where T50 T50z7 t and TLO TL0z7t If we formally expand the latent heat boundary condition we nd that at 016 that aTso aTLO i k Z 62 L l 62 which is identically satis ed since the leading order temperatures are independent of Z At 01 we nd that 6T 6T 6T 6T LZVbz7Z7t ksZlt 50 aglgtimzgtlt 60Tgt 9113 9112 0 ksZ However7 from the bulk equations we know that 6T0 6T1 i A0ltZ7t 32 dz 7 KZ 9114 so the boundary condition becornes LZ27Z7t A 2717 igA mL 9115 i 7167 190 ltigt716TL0 9116 LZVb27z7t K5Zlt Sgt 62 2 L 6 27 If we now divide through by LZ7 note that H kpcp and average the equation we nd that PSCps 1 15T50 PLCpL 1 15TL0 V 7 i i i 9117 ltogt lt L gtltKSgt 62 lt L gtltKLgt 62 or if pgcpg chpL pop is uniform then the boundary condition is grim ltkisgt1eltjgt1 9118 August 21 2008 Industrial Math Notes DM Anderson 96 10 Fluid Dynamics There are a number of references appropriate for this section Batchelor An Introduction to Fluid Mechanics 7 Tritton Physical Fluid Dynamics 7 Van Dyke An Album of Fluid Motion 7 Schaums Outline Series on Fluid Dynamics 7 and Fowler Chapter 6 to name just a few Many of the notes that follow are based closely on notes from graduate courses given at Northwestern University by MG Worster7 SH Davis and JB Grotberg during 198971993 We will describe a continuum approach to modeling a uid In this setting7 the uid is characterized by mean properties of a large number of individual molecules That is7 a uid particle77 refers to an in nitesimally small region mathematically but actually encompasses many individual molecules physically Each uid particle has associated properties such as temperature7 density and viscosity for example and has a well de ned position and velocity 101 Lagrangian and Eulerian Descriptions In the continuum approach7 there are two basic ways in which we can describe the motion of a collection of uid particles 7 the Lagrangian description and the Eulerian description 7 which differ fundamentally in the point of view from which the uid motion is described LAGRANGIAN In the Lagrangian description we keep track of the position7 velocity7 density7 etc of individual uid particles as they move Here we write 13 12207 t to indicate the uid velocity associated with the particle that was initially at location 20 In this sense7 the Lagrangian description is particle based Similarly7 we write p pi 07t to indicate the density of the uid particle that was initially at location 20 We de ne pathline of a uid particle as the path that that particular uid particle takes through space Each uid particle therefore has a different pathline Since two uid particles may occupy the same location in space at two different times we observe that pathlines can cross Pathlines may also change with time so that if we release a continuous stream of markers in a ow eg a stream of dye from a point source we observe a dye streak showing7 at a snapshot in time7 the location of the dye that has been released Since we 7follow7 speci c uid particles here7 the acceleration of a uid particle is given by d ail3t This is the usual acceleration formula that one would typically apply in a classical Newtonian physics problem eg What is the acceleration of an apple falling from a tree An incompressible uid is one in which the density of each uid particle remains constant in time In the Lagrangian description this is expressed as dpdt 0 EULERIAN In the Eulerian description7 we view the uid velocity7 density7 etc as functions of spatial coordinates xed in a laboratory frame of reference Here we write aalt2tgt p plt2tgt 101 August 21 2008 Industrial Math Notes DM Anderson 97 where f is a spatial coordinate in a xed frame of reference eg f z7y7z in Cartesian coordinates In contrast to the Lagrangian description7 we do not 7track7 speci c particles but rather observe particles passing through a xed location In this setting7 at a xed point in time7 we have a velocity vector 17 de ned at all spatial locations Streamlines are everywhere parallel to the velocity vector Unlike pathlines7 streamlines cannot cross except at sources or sinks The Eulerian description is usually a more convenient mathematical description as it allows for a eld variable eg the velocity vector 17 de ned everywhere in space and time The dif culty with this point of view is that since we do not follow a speci c uid particle we need to carefully modify the familiar7 Newtonian mechanics which characterize the motion of a particle as we follow that particle in order to write down things like acceleration and ultimately F ma With this in mind7 let us derive a formula for acceleration that is valid for the Eulerian description Consider a uid particle at location i at time t whose uid velocity is given by 1227 t After a small amount of time St the same uid particle is at a new location 2 62 Therefore7 the acceleration of this uid particle is given by w 6m it 7 m t 1x7 t align 6t 102 Using a generalized Taylor expansion we nd that ad up 6923 t 6t 1M t 692 V i 7 t Stag 06t2 6932 103 Noting that in the limit St a 07 the quadratic terms vanish and that Sf5t a 1257 t7 we nd that a f6f7t6ti f7t i a a 612 am 0 i 6133 f i u W a 104 Therefore7 the acceleration is given by D1 612 a E 7 n x 105 a Dt 6tltu Wu lt gt where the operator D 6 7 E 7 i 106 Dt 6t 1 V7 l is called the material derivative or substantial derivative While the Eulerian description is nicer mathematically in many ways than the Lagrangian description7 the nonlinear term 12 V1Z in the acceleration generally makes things dif cult but also very interestingl We shall focus our attention on the Eulerian description and in the next few sections begin to develop mathematical equations describing uid ow August 21 2008 Industrial Math Notes DM Anderson 98 EXAMPLE 101 A two dimensional discrete analog Consider a collection of ping pong balls or marbles rolling around on a table top A Lagrangian description of the motion would rst identify the initial position of each ping pong ball and then track each balls position velocity and mass suppose that each ping pong ball could somehow change its mass as it rolls as a function of time Here each ping pong ball has its own pathline An Eulerian description of the motion would involve describing the velocity and density of the ping pong balls as they appear at a speci c location 7 eg suppose you de ne a su ciently ne rectangular grid and keep track of the particle mass and velocity at each bop of the grid EXAMPLE 102 Describe in words a one dimensional discrete analog of motion using I a Lagrangian description and an Eulerian Description 102 Conservation of Mass Consider an arbitrary control volume V with boundary 3V xed in the laboratory frame of reference and through which uid ows Let n be the outward unit normal to the volume A mass conservation statement can be stated as The rate of change of mass in V is equal to the mass ucc through its surface 3V Mathematically7 we write this statement as d a A aVpdl iav pundS7 107 where p is the density of the uid and J the uid velocity The control volume is assumed to be xed in space and time Consequently7 the time derivative on the left hand side of the above equation can be moved inside the integral Further7 using the divergence theorem and rearranging the expression leads to the result A Finally7 since the control volume was chosen arbitrarily7 the integrand must be zero at every v my dV 0 108 point in the domain Therefore7 we arrive at the so called continuity equation for mass conservation 6 i v pa 07 109 6t which applies at each point in space and time in the uid Expanding the divergence term in this expression V 1 Vp pV 13 eg see appendix on index notation and using the de nition of the material derivative 106 allows us to write the continuity equation in the form DP 7 0 1010 Dt pV u August 21 2008 Industrial Math Notes DM Anderson 99 The condition of incompressibility is expressed as DP 7 0 1011 Dt lt gt that is7 the density of a uid particle is constant Note that in this Eulerian description this does not imply necessarily that the density at a particular location in space is constant Using equation 1010 implies that V0 0 1012 A ow that has V 11 0 is sometimes referred to as 7solenoidal7 Generally7 for slow ows ie where the uid velocity is much smaller than the speed of sound incompressibility is a reasonable assumption 103 The Stress Tensor It can be shown Batchelor7 Aris7 Tritton that there is a linear relation between the stress force per unit area exterted by the uid on a surface and the unit normal 73 of the surface pointing into the uid so that F 0 7 1013 where F is the stress and 0 is the stress tensor The stress tensor includes both pressure like forces and friction like forces 1031 Newtonian Fluid A Newtonian uid is one in which the uid is isotropic and in which there exists a linear relation between the stress and the rate of strain This leads to a symmetric stress tensor 0 ipl d7 1014 where the rst term represents the isotropic part and the second term represents the deVia toric part The deViatoric part has the property that trd 0 where tr is the trace Using index notation we write UM dij7 where 61 is the Kronecker delta de ned to be 1 ifz39 j and 0 ifz39 31 j The conditions that 1 d is symmetric and 2 that trd 0 in index notation dii 0 along with 3 the assumption that dij is a linear function of V1 ie uihj lead to 7 3n 3117 2 31 div 7 116 6 7 167166 1016 August 21 2008 Industrial Math Notes DM Anderson 100 01 2 dultV V T igmvom 1017 where u is the dynamic viscosity of the uid In the simple two dimensional ow con guration in which the uid velocity is purely in the horizontal direction so that 13 uy70 we can identify7 for a Newtonian uid7 a scalar stress 739 5792 Here we see that the stress 739 is linearly related to the rate of strain ouoy the proportionality constant is the dynamic viscosity M This constant represents the slope of the stress vs rate of strain curve as shown in Figure 5 Those familiar with elastic solids will note that in that case the stress is linearly related to the strain rather than the rate of strain 103 2 NonNewtonian Fluids For non Newtonian uids7 there can be a nonlinear relation between the stress and the rate of strain Cases of shear thinning7 and 7shear thickening7 uids are shown in Figure 5 We can interpret the local slope of the stress versus rate of strain curve as an apparent viscosity7 Therefore7 in the shear thinning case the viscosity decreases with increasing stress or increasing rate of strain An example of a shear thinning uid is nondrip thixotropic paint 7 it does not tend to drip but moves easily when you brush it ie apply a stress In the shear thickening case the apparent viscosity increases with increasing stress or increasing rate of strain An example of a shear thickening uid is a corn starchewater mixture try mixing 34 cup of water with approximately 20 tbsp of cornstarch The mixture acts much like a liquid when you stir it slowly but acts much like a solid when you try to stir it fast or quickly pull a spoon out of it A Bingham plastic see Figure 5 has a yield stress the stress vs rate of strain curve does not pass through the origin In this case a nite stress must be applied before the uid moves An example of such a uid is toothpaste 104 Momentum Equation Again consider an arbitrary control volume V xed in space and time A momentum balance statement for the uid in the control volume is given by The rate of chartge of total momerttum m V is equal to the momerttum flucc through its surface plus the body forces acting or the fluid plus the forces acting or the surface Mathematically7 this is expressed as ipodv 7 p dS fdv a dS7 1018 dt V 3V V 3V August 21 2008 Industrial Math Notes DM Anderson 101 Newtonian Shear Thinning w w w w 9 9 17 17 rate of strain rate of strain Bingham plastic SnearTnickening stress stress rate of strain rate of strain Figure 5 Stress versus rate of strain curves for Newtonian and Non Newtonian uids where frepresents a body force eg gravity Basically7 this is a force balance dPdt F where P is momentum and F the sum of all forces Noting that d a 7 0W amp vadV 7 V at dV 1019 and 3Vpuu ndS V V pu udV V u Vpu puV dV7 1020 the momentum balance becomes 6W 1 ATdV7A7uVpupuVufVadV 1021 Now7 arguing as before since the control volume is arbitrary7 the integrand must be zero everywhere in space and time so 0W at aVpapavafva 1022 Rewriting this by expanding the rst two terms on the right hand side gives 61 49p a a u u u u a piuiuuVppuVupuVufVa 1023 3t 3t August 21 2008 Industrial Math Notes DM Anderson 102 A further rearrangement of terms gives 6 6 a plt6izuvugtultaigmvppvugtfva7 1024 x so that the momentum equation becomes 12 a 7 1025 P Dt V a f 1f the uid is incompressible so that V 13 0 then 039 mud 421 ltV v12 1026 and V 0 Vp uvzm 1027 see appendix on index notation for derivation Putting this into the momentum equa tion 1025 gives D J 2d a p in uV u f7 1028a V ii 0 1028b These are the Navieristokes equations for the incompressible case In three dimensions they represent four equations for four unknowns ii U707 w and the pressure p The uid density p and dynamic viscosity M are material properties It is often convenient to de ne a kinematic viscosity V Mp which has units of length squared per time For a body force such as gravity7 7 p where is the acceleration due to gravity typically if gravity oriented vertically downwards7 g where I is a unit vector in the zidirection Gravity is a conservative body force and we can write 7 7V0 where p ip liz pgz In such a case it is possible to de ne a modi ed pressure p p p and write the momentum balance as Di p5 Special care must however be taken the density p is not constant ie the variation of density with temperature or uid composition is often an important effect that can cause Vp uvzu 1029 uid motion 105 Vorticity The vorticity is the curl of the uid velocity 13 V gtlt ii 1030 Using index notation we have M eijkdukdxj An equation for the vorticity can be found by examining the curl of the momentum equation eg see Batchelor7 p 267 August 21 2008 Industrial Math Notes DM Anderson 103 106 Boundary Conditions Typically7 uid ows are con ned by boundaries The ow of blood inside an artery is con ned by artery waII7 the ow of air in the atmosphere is con ned on much larger length scales by the earth7s surface The uid motion occurring inside a drop of rain sliding down a window pane is con ned by the interfaces between the water and air and between the water and glass Often in uid mechanics the locations of these interfaces are not known in advance and must be determined as part of the solution to a mathematical model of the ow In this section we introduce boundary conditions associated with mass conservation7 the so caIIed no inp condition and stress jump conditions that arise in uid ows MASS CONSERVATION A mass conservation boundary condition at an interface between uid 1 and uid 2 can be derived via a piIIbox7 argument Consider the continuity equation 3 lvpa 07 1031 at and integrate over a xed control volume that contains a moving uid uid interface The result of this is the interfacial jump condition P1Wl I7 73 PZWZ I fh 1032 where 7 is a unit normal vector to the interface and I7 is the local velocity of the interface This boundary condition states that there is no jump in mass ux across the interface In general it allows for mass ux across the interface as would occur for example at an evaporating interface If it is known that there is zero mass ux across the interface7 then this condition reduces to 121 71 132 71 In this case the interface is referred to as a material interface7 and the normal uid velocity on one side is matched by that on the other side of the interface In the case of a uid solid boundary with no mass transfer across the interface the mass balance boundary condition is a A A un 77 1033 Ill where I75 is the velocity of the solid boundary At a non moving solid boundary this reduces to the condition that the normal uid velocity is zero there 13 0 NO SLIP CONDITION This condition says that for viscous uids there is no jump in the tangential velocities across a uidi uid interface 0 1034 where f is a unit tangent vector to the interface This boundary condition can be derived from the momentum balance when viscous terms are present Essentially this condition ensures that there is not in nite stress derivative of velocity August 21 2008 Industrial Math Notes DM Anderson 104 At a non moving solid boundary this no slip condition is 1216 0 1035 in the uid at the solid boundary The no slip condition does not apply for inviscid uids modeled with u 0 lnvis cid uids are not physical uids7 although the inviscid limit is often realized away from boundaries in the ow and is an important and appropriate description of many uid ow situations STRESS BALANCE The jump in the stress at a uide uid interface in its most simple form is written as 039 7 yV 731m7 1036 where y is the surface tension between the two uids and V 71 is the interface curvature This is vector balance of which we are generally interested in the normal and tangential components given by 7V 7 1037a 0 103710 These conditions can be modi ed by other effects not discussed here For example7 gradients in surface tension along the interface would produce Marangoni stresses and correspondingly a term proportional to V7 in the tangential component of the stress balance condition A simple situation in which one can interpret the normal component of the stress balance is for a two dimensional interface parameterized by z h with the phase in z gt h and the 7 phase in the region 2 lt h in equilibrium Here hm 1 W 1038 and 039 71 710737 1039 If hm ltlt 1 then V 71 ihmm 1040 and the normal component of the stress balance is apt 7H 19 719 vhm 1041 Here we see if 10 gt p then we expect the interface to be convex towards the phase7 as is consistent with hm gt 0 August 21 2008 Industrial Math Notes DM Anderson 105 107 Recap Navier Stokes Equations The Navierestokes equations for incompressible ow are given by 5 1 673mm izvavzm 10425 W 12 07 104210 plus necessary boundary conditions and initial conditions In three dimensional Cartesian coordinates with uid velocity 13 U707 w and pressure p these equations determine four unknowns u7 117 w and p and are written 3U du du du 1 3p Euava 105 izaVV2U7 1043a uvwg 7VV2117 1043b uvw Z szw7 10430 2 2 27 0 1043d 108 Nondirnensionalization Let us now make equation 1042a dimensionless by introducing the following dimensionless quantities UW7 zyzgtiLiltzzycz gt wanUm paw 1044 where U is a velocity scale7 L is a length scale7 P is a pressure scale and the primed quantities are used to represent the dimensionless variables ie 13 and U have dimensions of length per time while 12 is dimensionless Many ows have what can be identi ed as a 7typical7 scales L could be the radius of a pipe or the width of an airfoil7 U could be the speed of a stirring stick moved through a container7 or P could be a externally imposed pressure difference Other choices for these scales to balance terms in an asymptotic sense 7 ie to balance the pressure and viscous terms in the momentum equation or to balance the pressure and time dependent terms The details of each problem under consideration dictate how the nondimensionalization is done Combinations of these scales and other parameters in the problem such as viscosity or density lead to dimensionless parameters Possibly the most important of these in uid dynamics is the Reynolds number Re Below we show how the Reynolds number arising in the nondimensionalization ofthe Navier Stokes equations We substitute the expressions 10 44 into equation 1042a gives U aw i H P H WU m LlU ltWu 39Vugt imvp WV U7 1045 August 21 2008 Industrial Math Notes DM Anderson 106 where V indicates differentiation with respect to dimensionless quantities Rearranging terms leads to W Ll Now suppose that we choose the pressure scale P pU2 this is one possible choice for D P u 711W 27 W pmz v 1 1046 P which effectively 7balances7 the pressure gradient term with the terms on the left hand side of the equation This leads to the dimensionless momentum equation DU 1 2a Dy 7 7V1 EV u7 1047 where the Reynolds number is de ned by Re UM 1048 V Physically7 the Reynolds number represents the ratio of inertial to viscous forces lf Re gtgt 1 corresponding to a relatively large velocity7 relatively large length scales7 andor relatively small kinematic viscosity 7 ie where inertial forces dominate viscous forces this limit leads to the Euler equation D17 Dt 7V p 1049 In certain settings this limit may need to be re assessed near boundaries where viscous forces may not be negligible The limit Re ltlt 1 corresponding to a relatively small velocity7 relatively small length scales7 andor relatively large kinematic viscosity 7 ie where viscous forces dominate inertial forces we cannot write V Z 0 7 as a general rule the pressure term must remain in the balance Instead we need to choose a different pressure scale than the choice made above In other words7 we rescale the pressure to balance the viscous term so 15 Rp This gives D17 1 Dt E 7V p V Z 1050 Here we can examine the limit Re ltlt 1 to obtain 0 7V p WW 1051 which is called Stokes equations In this case7 the time derivative is no longer in the equation and time dependence of the ow may enter the problem through boundary conditions 109 Some Simple Flows Here we examine some simple solutions of equations 1042 with boundary conditions as speci ed for each case August 21 2008 Industrial Math Notes DM Anderson 107 1091 2D Parallel Flows Here we examine ows for which there is no vertical motion and for which the the horizontal velocity is a function of the vertical coordinate and time only That is the uid velocity vector has the form U uy t 0 First we examine the continuity equation 1042b aiu 671 0 1052 3m 3y Observe that this is satis ed for any ow of the form speci ed above We say this equation is identically satis ed Next examine the momentum equation 104221 The x component of the momentum equation is 3U 3U 3U 1 3p 321i 321i 7 7 7 777 7 7 1053 at u6zvdy p6ylt6z2 3y Due to the above assumed form of the velocity eld we nd 3U 1 3p 3 7 777 7 1054 at p h ydyz The y component of the momentum equation is 31 31 31 1 3p 321 321 7 7 7 777 7 7 1055 at u6v6y p6yylt6x26y2 Due to the above assumed form of the velocity eld we nd that 6p 0 77 1056 6 lt gt which tells us that the pressure can be a function of z and It only pxt Further taking the x derivative of equation 1054 and noting that dudz 0 we nd that 3210de 0 which tells us that the pressure gradient is constant Therefore we write 619 7 E 7 1057 65 where G is assumed to be given or otherwise obtained from boundary conditions Therefore for parallel ows in two dimensions the equation governing the horizontal uid velocity uy t is a 7 G 62 10 58 p at M 611239 39 Note that an adverse pressure gradient drives a ow That is uid moves from regions of high pressure towards regions of low pressure August 21 2008 Industrial Math Notes DM Anderson 108 STEADY PARALLEL FLOWS 1 Couette Flow Couette ow is ow between two horizontal parallel plates in which the lower plate at y 0 is stationary and the upper plate at y d moves horizontally with speed U0 In this setting there is no applied pressure gradient G 0 Here we must solve 3211 0 7 1059a u 0 on y 07 1059b u U0 on y d 1059c The two boundary conditions are no slip conditions on the tangential component of the uid velocity at the two plates This problem has the unique solution 110 77 1060 which is a simple shear ow2 Recall that the stress exerted by the uid on a surface is F 039 ip MV V T For parallel ows an an 31 i i 37 t i 2a 31 O 31quot 1 A H 0 CT H V A i 0 A7 1 n7 1J7 ti07 1062 so that 0 a 7 3y 739 p1u0 1063 Then7 the tangential stress is r t a 10 64 739 7 May Therefore7 in Couette ow the tangential stress is A U Ft 1065 d This result is used to measure the viscosity of a uid in a shear viscometer 2Planar Couette ow is stable to linear disturbances but has possible subcritical instabilities 7 see Drazin and Reid7 pl 453 August 21 2008 Industrial Math Notes DM Anderson 109 1 Steady Poiseuille Flow Steady Poiseuille ow is ow in a pipe with an imposed pressure gradient The downstream end of the pipe is held at pressure p2 and the upstream end of the pipe is held at pressure p1 The pipe has radius a and length L It is assumed here that the horizontal uid velocity is a function of the radial coordinate measured from the center of the pipe In this geometry7 we identify the pressure gradient as G 7p2 7p1L The momentum equation in cylindrical coordinates is 1 d du 0 G u 7 1066a u 07 at r a7 1066b along with the condition that u is bounded at r 0 Integrating equation 1066a gives 00 G 2 u 77 7 7r c 1067 2T2 17 where 00 and 01 are integration constants Applying the boundary conditions implies that co 0 and 01 GrizZip Therefore7 steady Poiseuille ow is given by G 2 2 7 7 1068 W 4 a r lt gt which is a parabolic flow pro le3 The volumetric ow rate associated with this ow ie with this given pressure gradient 27r a Q urrdrdt97 0 0 2 1 7T 012 7 7 27quotdr7 0 4n 7 7Ta2G i 7Ta4AP 10 69 7 8n 7 8uL 39 39 Note that Q has units of length cubed per time 1010 Approximate Solutions to NS Equations Despite the exact solutions to the Navier Stokes equations shown above7 most uid ows of interest cannot be solved exactly Many ows require direct numerical simulation Often7 however7 approximate solutions can be obtained based on asymptotic analysis of the original equations We describe some such situations in the sections below 3For stability of such ows see Batchelor7 pl 186 and Drazin and Reid7 pl 190 August 21 2008 Industrial Math Notes DM Anderson 110 10101 General Lubrication Theory For ows in thin liquid lms eg a tear lm on the eye ows through relatively long and narrow gaps eg the uid region between a bearing and a solid boundary or more generally for ows in which a typical7 length scale in the vertical direction is much smaller than a 7typical7 length scale in the horizontal direction including even large scale con gurations such as lava ows lubrication theory can be used to simplify the governing equations In some cases this allows analytical solution to the ow In other cases the governing equations plus boundary conditions can be reduced to a simpler numerical problem ie reducing the coupled set of nonlinear PDEs to a single PDE To begin we consider the Navier Stokes equations in two dimensions 3U 3U 3U 16p 321i 321i 7 7 7 77 7 7 1070 at u6z 12 pdz V dz 3y 7 a a a a 4 an n 1070b 6tu6zydy T p6yy 6z26y2 g7 39 du 31 7 0 10700 where we have included gravitational forces Let p p pgy SCALING ANALYSIS Now suppose that in the problem of interest there is a vertical length scale H and a horizontal length scale L such that HL ltlt 1 Further suppose that the horizontal velocity component u has typical magnitude U the vertical velocity component 1 has typical magnitude V and the pressure 19 has typical magnitude P Let T represent a typical time scale of interest Through a scaling analysis we will identify the simpli ed lubrication equations and the conditions under which they follow from the original equations Scaling 7 continuity equation Consider the equation 1070c Based on the above as sumed scales a scaling analysis tells us that these terms have typical magnitude U V f 1071 Note that each term represented should be dimensionally correct The colon is used here simply to distinguish different terms in an equation In order to balance both terms together we need V N U ltlt U 1072 if HL ltlt 1 This tells us that the ow is nearly parallel 7 that is the horizontal velocity is much larger in magnitude than the vertical velocity This balance gives the full equation 10 70c August 21 2008 Industrial Math Notes DM Anderson 111 Scaling 7 x component of momentum equation Consider the equation 1070a A scal ing analysis gives U U2 UV P VU VU 7 7 7 7 7 7 1073 T L H pL L2 H2 Multiplying each term by HZVU and using 1072 H2 UHZ UHZ PHZ H2 7 7 1 1074 VT VL VL pVLU L2 Note that in this form7 each expression is dimensionless and can be compared to the term of unit magnitude For example7 based on our initial assumption HL ltlt 1 we know that HZLZ ltlt 1 The relative magnitude of the remaining terms depend on choices made for P7 U and T For lubrication theory we make the following choices 0 In order to retain the pressure term in the leading order balance for the horizontal ow in a similar way to how it balances in the parallel ow solutions we choose 7 uLU P H2 1075 0 We expect a time scale based on the horizontal velocity and the horizontal length scale and so choose T 3 1076 0 The above choice means that the rst two terms look like UHZ VL39 1077 For lubrication theory7 we want these terms to be small relative to 1 so we require 77 ltlt 1 1078 Under these conditions7 the leading order x component of the momentum equation is 310 3211 Scaling 7 y component of momentum equation Consider the equation 1070b A scal ing analysis gives UV VZ PAVV VV LHpH 1080 August 21 2008 Industrial Math Notes DM Anderson 112 Using the above choices for T7 V and P we nd that UZH UZH UZH VUL VUH VU 1081 L2 39 L2 39 L2 39 H3 39 L3 39 HL39 39 Next7 dividing each term by VULHS gives UH4 UH4 UH4 1 H4 H2 1082 VL3 39 VL3 39 VL3 39 39 L4 39 L2 39 39 NUHVH3L3ltlt1 Since HL ltlt 1 and UHVHL ltlt 1 by assumption7 the dominant term in this equa tion is the y derivative of the pressure Therefore7 the y component of the momentum equation reduces to 310 0 7 1083 6y In other words7 p p zt so p p zt 7 pgy Therefore7 under the two conditions 0 the geometry is 7thin7 so that H 3 ltlt 1 1084 0 the Reynolds number UHV is not too large UHH ltlt 1 10 85 V L 7 the lubrication equations are du 31 i 0 1086 dd 3y 7 a 3 3p i 1086b ayz 5957 p MW 7 My 1086c These equations apply in the bulk uid 7 boundary and initial conditions must still be applied We address below speci c problems in which the lubrication equations apply EXERCISE 101 Derive the lubrication equations for an accisymmetric flow ie in cylindrical coordinates where the radial coordinate is in the horizontal direction August 21 2008 Industrial Math Notes DM Anderson 113 1011 Axisymmetric Gravity Current Here we develop a model for the spreading a gravity current on a horizontal boundary We shall later use this model to understand and predict the spreading of corn syrup on glass This problem was rst analyzed in a series of papers by Huppert l3gtl 1gtlL5 Consider an axisymmetric gravity current occupying the region between the liquid air interface at y hrt and the liquid solid boundary y 0 The interface y hrt is an unknown and must be determined as part of the problem 7 in fact it will be the quantity of primary interest as it will tell us where the moving uid region is located We shall assume that hrt is a single valued function of the radial coordinate 7 although in reality this may not be true very near the leading edge contact line of the moving uid region As we shall see below our model of the axisymmetric gravity current will eventually neglect the effects of surface tension which are important for a detailed description of the contact line region as well as other ow regimes In the uid region the lubrication equations are 16ru 31 i a i 0 1087a 3211 310 672 i 57 1087b p p 73tpgy 10870 The full boundary conditions for this system are u 0 on y 0 no slip 1088a w 0 on y 0 no vertical ow 1088b 039 4amp4 yV on y hrt stress balance 1088c 13 ft on y hrt kinematic condition 1088d where 13 111 unit normal and tangent vectors are de ned by A 71771 A 17h t 1089 1 haw27 1 haw l and the interface velocity vector is de ned by 17 htk 1090 We shall now however simplify the boundary conditions consistent with lubrication theory Recall that in two dimensions the Newtonian stress tensor takes the form 231 311 V V T 373 941235 1091 August 21 2008 Industrial Math Notes DM Anderson 114 Therefore7 using the lubrication theory scaling we nd that W WV i 3 l 0HL 1092 32 Also7 a 071 0HL t 170 0HL 1093 Therefore7 the normal component of the stress balance condition is 71011 0HL 1094 so 10 r 0 on y hr7 t and the tangential component of the stress condition is du i O H L 1095 9 6y 1 lt gt lt gt so making the simplifying assumption that the air is inviscid we have the condition on the uid that dudy 0 on y hr7t If the ambient air pressure is pA constant atmospheric pressure then on y hr7t we must have pA 100 7 t 7 pgh Therefore7 pT7t pA pghr7t39 That is7 the form of the pressure is effectively determined by the interface shape hr7t Consequently7 the radial derivative of the pressure is given by 310 31 7 7 1097 6T M 6T Realizing that the pressure gradient can be written in terms of I17 we revisit the horizontal component of the momentum equation 1087b to nd that 3211 31 672 P937 1098 subject to boundary conditions u 0 on y 0 and dudy 0 on y hr7t This system can be solved for u in terms of h so 3le 777921 7 9 1099 M where hr7 t is still to be determined Next7 now that u is known in terms of h we can use equation 1087a to solve for i We nd that 116 Tu ur7y7t 70 Eardy7 10100 August 21 2008 Industrial Math Notes DM Anderson 115 which satis es the boundary condition 1 0 on y 0 The above results give us p7 u and i in terms of the interface position h The remaining condition that determines h is the kinematic boundary condition Using 13 Iii1251727 10101 17 7 1 127 10102 we nd that the kinematic condition can be written ht i 7 uh7 on y hrt 10103 Note that all three ofthese terms scale in the same way N UHL Therefore7 we apply 10 103 as the equation that determines h lnserting equation 10100 into this result leads to h ht 7 lawman71111 0 1 r 37quot 3 h 777 CI 10104 T 6T 0 m y lt gt where the second equality follows from application of Leibnitz rule for differentiation of integrals Now7 h all h M d if 7 217 d 0 Tu y afar 0 y y y igraih hi1 fh37 2M 37 3 1 1 351 7 g rh gr 10105 So the equation 3h p91 3 361 7 777 h 7 10106 3t 3ur6r T 37 7 determines the evolution of the interface position hr7t Note that this equation is now completely decoupled mathematically from the ow The interface shape hr7 t is subject to the condition that the total volume of uid remains constant so that TNlttgt v5 27r rhr7tdr7 10107 0 and the contact condition hrNt7 t 07 where rNt is a quantity to be determined Below we discuss appropriate initial conditions Scaling Analysis of Reduced Model It turns out that there is a similarity solution for this problem and we shall solve it in that way here In order to see this and7 in fact7 predict August 21 2008 Industrial Math Notes DM Anderson 116 the correct similarity scaling we rst do a scaling analysis on the reduced problem given by equation 10106 and condition 10107 Let H represent a typical vertical length scale7 R represent a typical horizontal length scale and T represent a suitable time scale Here we are seeking a balance of terms in the equation we7ve already reduced the problem and so examining the equations should suggest relationships among T7 R and H From the evolution equation 10106 we have H pg H4 7 77 10108 T M R2 balancing these two terms suggests that VRZ T 7 10109 9H3 lt gt Next7 from the volume constraint 10107 v5 RZH 10110 or interpreting this condition as giving a vertical length scale we have Vb H 10111 Solving equations 10109 and 10111 for H and R in terms of T reveals that 14 3 18 H 7 R 9 T 10112 gT V These two results actually provide extremely valuable information about the dynamics of the uid 7 they indicate how R the typical drop radius and H the typical drop height scale with time in this lubrication limit That is7 based on this analysis alone we can predict that R N T lS7 or 1 lnR glnT 10113 Careful experiments measuring R and T should be able to test the prediction of slope 18 between lnR and ln T We can learn more about this scaling if we revisit the conditions HL ltlt 1 and UHVHL ltlt 1 under which lubrication theory is valid From our solution for u we identify U gHSVR Requiring HL ltlt 1 we nd from equations 10109 and 10111 that V Tgtgt 7 10114 9161 August 21 2008 Industrial Math Notes DM Anderson 117 Requiring UHVHL ltlt 1 we similarly nd that 13 T gtgt E 10115 gV We will see below that for our parameters of interest the second condition is stronger Other than these inequalities we have no further balances that tell us how to determine the time scale This again suggests the possibility of a similarity solution for this problem Similarity Solution for Reduced Model The above scaling analysis suggests the possibility of a similarity solution to equation 10 106 and associated condition 10107 and condition hrNt7t 0 based on the following similarity variable and rescaling of the interface height h V 18 T VV 14 7 7 h t J 10116 n M M 72 gt M lt gt where f77 is a function of 77 to be determined We reformulate equation 10106 through a change of variables Note that unlike the similiary solution presented earlier for the heat equation solution7 the independent variable h here depends explicitly on both 77 and t First7 an 1 u 18 r n a g W 7 10117a 677 18 1 a 1011710 Transforming terms in 10106 we nd 9 6t 71754 EM mml 1011821 amp m 37 T gt 77 3 VV 14 V 18 1 M ny 10118b August 21 2008 Industrial Math Notes DM Anderson 118 3h VV TIRE nfgf n7 101180 3 136771 i i Vivi f3 L 18 10118d 37 T 37 i dn gt 77 9163 t187 39 Therefore7 equation 10106 becomes 1Vl3141 1 791d Vlb 3 V 181 3607 g lel g F W f 93 W7 1 1 1 d 77 1f77 gUfnll Ufol7 10119 or 1 1 2 11 s 4nf77877 f773dn nf f i 07 1 1 g 7721 nfsf 0 10120 The volume constraint 10107 can also be rewritten in terms of 77 If we de ne V 18 TN W E W7 10121 and note that 1 8 dn 51 10122 the volume constraint becomes 1 27r0W nfndn 10123 Summarizing these results gives us the problem ltn2fnf3f 07 10124a 27r0W77f77d77 17 1012410 mm o7 lame plus a boundedness condition as described below To solve this problem we rst integrate equation 10124a to obtain 1 1 n2f nf3f A 10125 August 21 2008 Industrial Math Notes DM Anderson 119 where A is an integration constant However if we evaluate this expression at 77 0 and assume that f 0 is not unbounded we must have A 0 This gives 77f m 0 10126 In order to satisfy this equation the quantity in parentheses must vanish This condition can be written as 1 2 1 3 7 7 0 10127 1677 l 9 f l 7 which can be integrated once to give 1 2 1 3 7 7 B 10128 1671 9f 7 where B is an integration constant Using the condition 10124c we nd B nJZV16 so that 13 10 g 10129 From this result we can con rm that f 0 is indeed bounded However f is in nite at 771V indicating that we have left out some physics at the contact point This tells us that lubrication theory breaks down near the contact point where the interface shape and ow quantities are no longer slowing varying in the radial direction Despite this we can see that f3f 7712V 7 77213 a 0 as 77 a 771V and that f a 0 there Finally the volume constraint 10124c can be used to determine 771V Here 9 13 711V 1 27f i 7707120 7 77913610 16 0 9 13 3 2 2 43 27TltEgt lt gt lt0N U 0 7 3 9 13 ltEgt 771 10130 711V This gives nN 437r3816918 m 0779 So V3t 18 W My gt 10131 14 13 VV 9 13 hrt 70 T6 7727772 10132 Note that the result for rNt was already predicted up to a multiplicative constant by the scaling analysis result in equation 10112 August 21 2008 Industrial Math Notes DM Anderson 120 These results have been validated carefully experimentally by Huppert Recall that for our predictions to be valid we require V 13 V T gtgt 3 and T gtgt m 10133 9 glb For an experiment with corn syrup spreading on a glass surface7 using g N 1030771527 Vb N 10007713 and V N 100771257 we nd that the rst condition requires T gtgt 02sec and the second condition requires T gtgt 10 837 indicating that the rst is a stronger condition EXERCISE 102 Using the emperimental data generated in class try to verify the sealing r N t l8 for the azisyrnrnetrie gravity current 1012 Unsteady Parallel Flow A good reference for the topics to follow is Batchelor7 pp 1867195 10121 Oscillatory Flat Plate Consider a viscous uid occupying an in nite region y gt 0 above a solid boundary at y 0 that is oscillating with velocity U0 cos wt in the horizontal direction In this situation7 the ow is parallel to the horizontal boundary so that 17 uy7t70 There is no imposed pressure gradient From the incompressible Navier Stokes equations we nd that the con tinuity condition is identically satis ed7 the vertical component of the momentum equation is identically satis ed7 and all that remains is the horizontal component of the momentum equation Therefore7 the problem under consideration satis es 3n 3 7 7 10134 6t Veg27 l where V up is the kinematic viscosity with corresponding boundary conditions u a 07 as y a 007 10135a u U0 cosuit7 on y 0 10135b This is a linear PDE for the horizontal velocity u and the forcing7 or inhomogeneity7 comes from the boundary condition at y 0 Of interest will be to identify a steady state7 solution in which initial transients from any initial conditions associated with starting the plate motion have decayed away We can identify a solution to this problem by rst anticipating that a solution of the form u ReUOeiwthy 10136 August 21 2008 Industrial Math Notes DM Anderson 121 where Re7 indicates the real part and noting that hy must then satisfy z wh Vh 7 10137a h 7 07 as y 7 007 10137b h 17 on y 0 10137c This ODE for hy has exponential solutions expmy of the form h exp 402 10138 V That is7 a linear combination of these two solutions makes up the general solution Now7 since i1 we have 10 0 h ex if i 10139 plt V41 In order for h 7 0 as y 7 00 we must choose the exponential solution with the minus sign Further7 to satisfy the condition at y 0 we have 1 0 0 h ex 7 7 10140 0 p lt gt Returning to our expression for u shows that 1 0 0 R U wt 7 i u e 06 exp lt Vy l7 U0 exp cos wt 7 10141 We can interpret the above solution as a vertically travelling wave with wavelength 2 2 wavelength 7 27139 17 10142 0 V w and amplitude that decays with distance from the boundary The decaying amplitude in dicates a loss of energy due to the viscosity of the uid The rate of dissipation per unit volume of uid can be de ned as D 20 51751de 10143 V eg see Batchelor7 p 153 Since the only nonzero contribution to 617 is 612 621 37 we nd that 617617 ig72 ig72 2 so that the rate of dissipation per unit area 32 on the plate is 00 du 2 DA 00 a dy 10144 August 21 2008 E Industrial Math Notes DM Anderson Oscillating Flat Plate Velocity Profiles y 032v12 47 3 2 1 G l l 1 08 06 04 02 0 02 04 06 08 1 uU 0 Figure 6 Sketch of velocity pro les for an oscillating at plate for one period in time 122 August 217 2008 Industrial Math Notes DM Anderson 123 Of most interest is a time averaged rate of dissipation over one period of oscillation and so we examine 1 T DA 7 DAdt7 10145 T 0 where T 27Tw These calculations can most easily be handled with the solution in exponential form and so we write 011 2 2 web 2 g U0 Ree 10146 Now7 in general7 if a Re aoem 7 b Re beam 7 10147 then ab aRcoswtiaIsinwt bRcoswtibIsinwt 7 0 0 0 0 a b cos2 wt 7 15196 a b cos wtsin wt 1 a b sin2 wt7 1 1 i055 10507 1 ERe 10193 10148 Therefore7 1 00 61 611 D 7U2 7 d 10149 ltAgt 21100777767721 lt gt Now7 3h w w a 71z4 exp 71z4ygt7 10150 so aha1 i 2ltwgtex lt72 W 3y 3y 7 2V p 2Vy 7 w 20 iexp lt y 10151 V V Finally7 20 20 ltDAgt g Ur exp 7 0 din V 1 V 7 U2 7 2p Ow 2w7 1 V0 7 U207 10152 2p 0 2 7 August 21 2008 Industrial Math Notes DM Anderson 124 which looks like a kinetic energy pUg times a velocity factor Mylo2 7 that is7 a kinetic energy ux One can associate the dissipation with an irreversible addition of heat to the uid EXERCISE 103 Show that the above result for ltDAgt is equivalent to the time averaged rate of work W by the plate per uriit area oh the fluid where W is aforee per uriit area times a the plate velocity so that W 75y 0U0 cos wt where 75y 0 7ug7y 0 10122 Impulsively Started Flat Plate Consider a viscous uid occupying an in nite region y gt 0 above a solid boundary at y 0 that at time t 0 begins moving with constant speed U0 As in the previous problem7 the ow here is parallel to the horizontal boundary and we have 17 uy7 t7 0 This problem is formulated mathematically as 71 viiiZ 10153a u 7 07 as y 7 007 10153b u 07 everywhere for t lt 07 10153c u U07 at y 0 for t gt 0 10153d Note that this is mathematically equivalent to an impulsively heated boundary where u represents temperature and V represents the thermal diffusivity There are different ways to reach a solution of the above problem One could seek a solution via a Laplace transform in time ie an eigenfunction expansion with respect to the exponential e pt One could also seek a solution via a Fourier transform in space ie an eigenfunction expansion with respect to the exponential eiky Alternatively7 one could perform a scaling analysis that would reveal the possibility of a similarity solution to this problem Here we describe a Fourier transform approach to equations 10153 ln particular7 we write the Fourier transform and inverse Fourier transform as uyt ukteikydk 1015421 mm QLm uyte ikydy 1015410 7139 700 Note that here we have written the y integration on the interval 7007 00 while the domain of physical interest is 000 If we identify the initial condition uy0 2U01 7 where is the Heaviside function with H 1 for y 2 0 and H 0 for y lt 0 and note that the Fourier transform gives the mid value at a discontinuity this approach will provide the solution to the original problem 7 we go through the details below August 21 2008 Industrial Math Notes DM Anderson 125 Inserting the Fourier transform into equation 10153 we nd that the reduced problem is the ODE for 12k7t given by 61 6 sz a 10155 subject to the initial condition 120570 Luy 70e ikyldy 7 10156 where we now take uy70 2U01 7 as noted above Therefore7 we can write the solution 12167 t as mm k70equot k2t 10157 Therefore7 substituting these results into the Fourier transform 1015421 our solution uy7t can be expressed as oo 1 co uy7t LOO LOO uy70e ky dy E Vk2telkydk 10158 5050 Now7 inserting the initial condition for uy7 0 it follows that 00 1 0 I I I uy7t LOO LOO 2er lky dy equot k2telkydk7 E 00 0 eiikydy1671k2t6ikydk7 7T 700 00 E 00 57Vk2t5ikydk7 7T 700 0 E 00 00 Eikateikyeikgdk d5 7139 0 00 E 0 U eivk2tiky dk d5 7139 0 00 E U 67vtk2iiky Vtldk d5 7139 0 00 E 00 U ammomeelty gt2lt4vtgtdk d5 7139 0 00 E fewerm U e VtkiyE2vt12dk d 10159 7139 0 700 Now note that eevt1keiltysvlt2vtgt12dk eevtpzdp 1 0 2 7 6 5 d5 x Loo August 21 2008 Industrial Math Notes DM Anderson 126 8 2 esds o gerfcw 10160 l lw lw Therefore U0 7T 00 2 u t 77 NM lt4Wgtd 10161 6 gt T M 0 s lt gt Further manipulations of the above expression using the de nition of the complementary error function reveals that uyt errfcy2m 10162 The above solution for shows that the independent variables appear only in the combi nation This suggests that a similarity solution approach is possible In order to see this before actually identifying the solution above a scaling analysis could have been done We outline the details of such an analysis here The equations 10153 are linear and driven by a boundary term with velocity scale U0 Therefore anticipate that Mm Uofy7t 10163 where fyt is a dimensionless function satisfying 6f 62f i 7 10164 6 Val27 a f 0 as y a 00 10164b f 0 everywhere for t lt 0 10164c f 1 at y 0 for t gt 0 10164d Introduce a dimensionless length 77 yL where L is a length scale and a dimensionless time 739 tT where T is a time scale In general the function f must depend on dirnensionless variables so f f77739 lnserting this into the PDE gives 1 3f V 32f if 77 10165 T 3739 L2 3772 This equation suggests taking L VT The boundary conditions however do not provide any further information relating L and T In particular there is no intrinsic time scale the only time scale in the problem is the elapsed time itself Therefore choose T t which implies that f MINE 1 axW 10166 August 21 2008 Industrial Math Notes DM Anderson ImpulsiverStarted Flat Plate Velocity Profiles 01 02 03 04 05 00 uU0 Figure 7 Sketch of velocity pro les for an impulsively started at plate 127 August 217 2008 Industrial Math Notes DM Anderson 128 It is apparent here that the function f can only depend on the combination yx mt Therefore7 with u Uog77 where 77 yxVt we transform the original equations 10153 using 377 i 77 du i 1 3277 i H 1 a U09 77 ii a U09 WW7 672 U09 707 10167 leads to the ODE g ing 07 10168 subject to boundary conditions 9 a 0 as 77 a 00 and g 1 on 77 0 Solving this equation leads to g77 erfc772 Therefore7 u errfc772errfcy2 7 10169 as was discovered previously in equation 10162 using the Fourier transform method Using a similar scaling analysis7 one can anticipate the form of the tangential shear stress on the boundary y 0 ln particular7 739S nail3y quL In terms of the similarity scaling found for L we anticipate then that MUO 739S N 7 10170 M To test this7 examine the full solution 10162 which shows that 9U MUO 2 7 2 MUO 7395 7 0 if 6 114 7 10171 776W gt 2m 77 lt gt We see that the scaling analysis predicts the actual tangential stress up to multiplicative constant 7 in this case a factor of August 21 2008 Industrial Math Notes DM Anderson 129 References 1 m H no H o H H H m H H H an H 0 H H H 00 H co m o M Abramowitz and lA Stegun7 Handbook of Mathematical Functions With Formulas Graphs and Mathematical Tables7 National Bureau of Standards7 Applied Mathematics Series 55 Washington7 DC7 1970 Gl Barenblatt7 Scaling Self Similarity and Intermediate Asymptotics Cambridge7 1996 GK Batchelor7 An Introduction to Fluid Mechanics Cambridge7 1967 CM Bender and SA Orszag7 Aduanced Mathematical Methods for Scientists and En gineers7 McGraw Hill7 New York7 1978 HS Carslaw and JC Jaeger7 Conduction of Heat in Solids7 2nd Ed Oxford7 1959 SH Davis7 Theory of Solidi cation Cambridge7 2001 Drazin and Reid7 Hydrodynamic Stability Cambridge7 1981 LC Evans7 Partial Di erential Equations7 American Mathematical Society7 1998 AA Evett7 1966 Permutation symbol approach to elementary vector analysis AC Fowler7 Mathematical Models in the Applied Sciences7 Cambridge7 1997 EJ Hinch7 Perturbation Methods7 Cambridge7 1995 MH Holmes7 Introduction to Perturbation Methods7 Springer7 New York7 1991 Huppert7 HE 1982 Nature 3007 4277429 Huppert7 HE 1982 J Fluid Mech 1217 43758 Huppert7 HE 1986 The intrusion of uid mechanics into geology7 J Fluid Mech 1737 5577594 J Kevorkian and JD Cole7 Perturbation Methods in Applied Mathematics7 Springeri Verlag7 New York7 1981 Kurz and Fisher7 AH Nayfeh7 Perturbation Methods Wiley7 New York7 1973 Schaum7s Outline on Fluid Dynamics MR Spiegel7 1968 Mathematical handbook of formulas and tables7 Schaum s Outline Series McGraw Hill7 New York August 21 2008 Industrial Math Notes DM Anderson 21 22 23 24 25 26 27 81 Rubinow7 Introduction to Mathematical Biology7 Dover7 1975 DJ Tritton7 Physical Fluid Dynamics Oxford7 1988 A Urnantsev7 Crystallography 307 153 1985 A Umantsev and SH Davis7 Phys Rev A 457 7195 1992 Van Dyke7 Ari Album of Fluid Motiori DJ Wollkind and RN Maurer7 J Cryst Growth 427 24 1977 E Zauderer7 Partial Differential Equations 130 August 21 2008 Industrial Math Notes DM Anderson 131 Appendix A See also Schaum7s outline 77 Mathematical Handbook o Trig identities sinA i B sinAcosB i cosAsinB A1a cosA i B cosAcosB sinAsinB A1b sinAsin B cosA 7 B 7 cosA 37 A2a cosAcos B cosA 7 B cosA 37 A2b sinAcos B sinA 7 B sinA B A2c sin2 A g 7 cos 2A7 A3a cos2 A g cos 2A7 A3b cos3 A 2 cos A cos 3A7 A3c sins A 3 sin A 7 i sin 3A A3d Expansions for z lt 1 1 V1z 1 z A4 The error function and complementary error functions are de ned by 2 m 752 i 2 00 752 erf 0 e d57 erfcz7 m 6 d5 A5 These have the properties erf0 07 erfoo 1 and erfc 1 7 erfz Leibnitz Rule d ism 1504 6F d EO Fm am 0 an Fm wing A6 Appendix B lndex notation is a useful tool for manipulating vector and tensor quantities and various differential operations on them One reference on this is a paper by Evett 1966 August 21 2008 Industrial Math Notes DM Anderson 132 In this context we shall de ne 11 as the 2 component of a vector 6 and A as the i j component of the tensor A We de ne the Kronecker delta 61 as 1 ifz39 j and 0 ill 31 j We also de ne the quantity 1 if ijk12372317312 if ijk32172137132 B1 if any 2397 j7 or k are equal H Eijk It follows from these de nitions that emnerpq imprint 7 6mq6np B2 The idea is to use index notation to represent vector quantities ie Vp and use the representation to derive further identities ie V Vp V210 First7 the 2 component of the gradient of a scalar Vp is written in index notation as 6 W 7 7 226 B3 Quantities such as the dot inner product of two vectors 6 and lican be written as 3 39 g aibi albl agbg agbg Note here that in index notation7 summation is implied over indices that are repeated In a related fashion7 the divergence operator on a vector is written as 91M 6111 duz 6113 V a 7 7 7 B5 u lt 3x1 3 gig The 2 component of the cross product of two vectors 6 and l is given by 3 X eijkajbk In a similar way7 the 2 component of the curl of a vector 6 is given by a 60k V X 1 Eijkai j For example7 the rst component of the curl 239 1 is given by Bag 602 603 602 7 7 7 7 7 B8 61de 6132 gig gig 63 Note that repeated indices j and k in equation B7 indicate that there are also terms involving 61117 61127 61137 61217 61227 6131 and 6133 which all vanish by the de nition of 6 You