New User Special Price Expires in

Let's log you in.

Sign in with Facebook


Don't have a StudySoup account? Create one here!


Create a StudySoup account

Be part of our community, it's free to join!

Sign up with Facebook


Create your account
By creating an account you agree to StudySoup's terms and conditions and privacy policy

Already have a StudySoup account? Login here


by: Ms. Helen Sipes


Ms. Helen Sipes
GPA 3.98

Hailiang Liu

Almost Ready


These notes were just uploaded, and will be ready to view shortly.

Purchase these notes here, or revisit this page.

Either way, we'll remind you when they're ready :)

Preview These Notes for FREE

Get a free preview of these Notes, just enter your email below.

Unlock Preview
Unlock Preview

Preview these materials now for free

Why put in your email? Get access to more of this material and other relevant free materials for your school

View Preview

About this Document

Hailiang Liu
Class Notes
25 ?




Popular in Course

Popular in Mathematics (M)

This 20 page Class Notes was uploaded by Ms. Helen Sipes on Sunday September 27, 2015. The Class Notes belongs to MATH 690N at Iowa State University taught by Hailiang Liu in Fall. Since its upload, it has received 49 views. For similar materials see /class/214488/math-690n-iowa-state-university in Mathematics (M) at Iowa State University.




Report this Material


What is Karma?


Karma is the currency of StudySoup.

You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!

Date Created: 09/27/15
Literature Review Level Set Method for Computing Solutions to Incompressible TwoPhase Flow1 Wen Wang May 13 2008 A level set method for computing the motion of incompressible two phase ow was described in this paper The motion of air bubbles in water and falling water drops in air were tested with the developed method 1 Equation of motion Vp 1 1 U U VU 7 EV 2MD 7 Eaa6dn F V U 0 where F is body force After non dimensitionalization7 momentum equation is Vp 1 1 U u U if 7 2 D 7 6 d t V p pReV u pB0rltang where the dimensionless parameters 32 2 ffE Ml 2R xgp 4p 9R p M Re 1 f 307f Pf Hf Hf 7 Re is Reynolds number and B0 is Bond numberthe ratio of body force to surface tension force R is the initial radius of the bubble and g represents a unit garvitational force D is the rate of deformation tensor D VU VUT Rewrite the equation in operator form ULU7E p to Projection method Let V LU If the IV for momentum equation is well posed7 there exists a unique decom postion Hodge decompostion V Vd VLpp with Vd is divergence free and Vd 1 VLpp De ne a projection operator PV E Vd the momentum equation is just Ut PLU To calculate projection7 we can take the curl of both sides of the equation pV de th to obtain V X W V X de Since Vd is divergence free7 there exists a stream function 11 such that Vd V x II Then the momentum equation becomes V X W VpV 1 For 2d problem 11 07 0711 above equation can be easily solved for stream function 1 3 Level set description The level set function 1 is positive in liquid and negative in the bubble and the bubble interface is the zero level set of j The level set equation is rUV O with initial value 10 to be the signed normal distance function Therefore7 the nal total set of equation is U PltLUgt gtt 7U W L gt0 p mM lt0 Pb pf2pf7 if gt 0 L gt0 M MMW7 if gtlt0 WbMHM 7 0 4 Smoothing Since the Poisson equation will yield instabilities at the interface if a unsmoothed p is applied7 a smoothed equation is used instead PbP PbP p 2 f A f Pf 2W 1 if gt a p gt Pbva if gt lt 704 6 Apsin otherwise where 04 is the nite thickness of the interface 04 is a constant for all time while it decreases at a rate of 001 The surface tension force 1 gram0n can be cast in the level set formulation and smoothed using 1 1 EWM EM gt5lt gtV with the curvature V 1 V 7 1W1 The delta function 6 gt is approximated with a discrete delta function 1 15 W zerlltm llv 1W1 07 otherwise 04 is denoted as the thickness of the surface and the usual value is 3Az2 In step 27 the contribution of surface tension to V x pV is then 1 1 E ltHlt gt6lt gt ygtm 7 ltnlt gtgt6lt gt mgtyi 7 where the Heaviside equation is H gtmHy gtyHml g7 gta Haw 7 ltiia sin7r gtozgt otherwise 5 Keeping j a distance function After some period of time 1 can become irregular when a steep gradient in j arises Without reinitialization the mass of the bubble will not be conserved after some time Conventional routines for reinitializing a distance function is to explicitly nd the contour j 0 and reset 1 at all points close to the front This takes 0n3 operations and can distort the front One approach is evolve the equation 1 1 7 lV gtl until 1 reaches a steady state at each time step7 given a region 9 with j 2 0 on 9 and j 0 on 89 If not using a distance function7 the level set function 1 can be solved to steady state by gtt S gt01 Mi 2 We 0 gtox where S is the sign function and it needs to be smoothed as 0 S 7 6 0 62 The above equation has the property that 1 remains unchanged at the interface therefore the zero level set of 10 and j are the same Away from the interface 1 will converge to lV gtl 1 The total algorithm is summarized as Step 1 lnitizlize ltjgtx7 0 to be the signed normal distance to the front Step 2 Solve UtPUi tUV gt0 for one time step Denote the updated 1 as WM2 and updated U by Un12 Step 3 Construct a new distance function by solving 4 S l21ilV gtl gtltx 0 gt 12ltxgt Denote the steady state solution by gtn1 Step 4 The zero level set of gtn1 gives the new interface position and it s a distance function Step 5 Repeat step 2 and 3 a Discretization a Discretize in time 2nd order Adams Bashforth method m U k32PLU 712PLU 1 b Convection terms since V U 07 we have U V gt WM may U VU fm gy with 2 LL U U f7ltwgtygiltvzgt The computation of quZJ similarly for umHZ gti12j7 is using a 2nd order ENO scheme c Viscous and curvature term Viscous stress tensor D is approximated using central differencing and similaryly for 11 and W Similar fashion is applied for the curvature term d Discretization of projection To de ne the discrete approximation of pro jection7 the discrete divergence and gradient operators must be de ned rst The matrix system for 11 is solved by using a preconditioned conjugate gradient algorithm using an imcomplete Cholesky factorization as a preconditioner e Reinitialization of j The reinitializing equation for j is written as t w Wgt S gto W S gtoWgtlWgtl Computation of At obeys the CFL contions due to the convection terms7 restriction due to stiff source terms and viscous terms 2 A255 WW2 my mm 3 p 87139 14 7 u h M 1 Ate man U At m2nAts At Ate 7 Main results Air bubble ow with R5 100 1000 was studied under a range of Bond numbers from 200 to 25 The simulation results with Eulerian grid are very close to the Lagrangian scheme results A falling water drop problem was also studied and the collision of two water drops was simulated There is no extra code needed for handling merging7 breakup7 dilation or contraction of the interface References 1 Sussman7 M7 Smereka7 P and Osher7 S7 1994 A level set approach for computing solutions to incompressible two phase ows77 J Comp PhysVol 1147 146 159 Math 690N Final Report Yuanhong Li May 05 2008 Accurate tracking of a discontinuous thin and evolving turbulent ame front has been a challenging subject in modelling a premixed turbulent combustion Level set method can be applied to improve the remixed turbulent combustion simulation by providing a better capturing of this highly trasient ame front This report presents the application of the level set met od derived from the Geequation 2 to simulate a ame propagation in premixed internal combustion engine combustions 1 Problem description ln a premixed engine combustion it is Well accepted that a corrugated ame front divides the ow eld into a burnt region and an unburnt region ln the level set method this ame front corresponds to an isoscalar surface or Zeroelevel set surface as CzJ GD 1 As shown in Figure 1 G gt CD is the region of burnt gas and G lt 03 is that of the unburnt mixture The choice of CD is arbitary but xed for a particular combustion event The level set model will be derived from the above equation I e a Inlmrm Figurf 1 A schematic representation ofthe ame front as an isorscalar surfac Gm Gquot 2 Algorithm description The well known level set model G equation was initially introduced by Markstein 2 in 1964 In 1985 Williams 3 introduced the method to describe ame front evolution Later Peters 1 made signi cant contributions for modelling turbulent combustion using this approach In this section a transport equation for an instantaneous ame location G is rst derived from equation 1 followed by the derivation of the transport equation for the mean 21 Equation of Motion Di erentiating equation 1 with respect to 75 results in a transport equation for the scalar G 8G 15 at VG dt 70 2 where is the propagation velocity and equals the sum of the ow velocity 17f at the front and the burning velocity in the normal direction as dd 4emb a and the vector normal to the front in the direction of the unburnt gas 7 is de ned by VG if 4 Substitution of equation 3 and 4 into equation 2 gives equation of motion as 8G 55fVG8LVG 22 Reformulation The burning velocity SL is usually not well de ned since it can be complicated by many factors such as ame stretch and curvature By taking into account the e ects of ame stretch and curvature the motion equation is reformulated as 8G a EvfVG3TVG7DltVG 6 where ST is the turbulent turning velocity D is Markstein di usivity and n is the ame curvature de ned by VG niVniV7W 7 23 Level set description Mean quantities are usually solved in CFD applications ie mean pressure mean velocities and mean temperature are solved in the averaged Navier Stokes governing equations To be consistent with CFD convention a mean G is sought as well According to Peters 1 the equation for the mean location of turbulent ame front can be written as ac E where 7 is the velocity of the moving vertex u is the unburnt gas density and is the gas density at the mean location of the turbulent ame which is de ned by G5 t 0 and F is mean ame front curvature aiayv wqumWQ 24 Needed level set techniques The level set equation 8 will not be solved until the fully developed turbulent ame is es tablished as dictated by an ignition kernel length At this moment an initial value GKE 0 is required which is computed by using the signed distance from each vertex to the kernel front as follows When the cell vertex 24 is outside the ignition kernel 57 t 1W4 MW 1114 M14 214 Zpi427 9 and when the cell vertex 24 is inside the ignition kernel Ci 75 WW4 m4 1114 ypi42 214 Zpi427 10 where i4yi4and 22 4 are the coordinates of vertex 2394 and pi4ypi4and 2p2394 are the positions of the ignition marker particles In addition reinitialization of the level set function is needed in order to reduce the numerical errors caused by steepening and attening e ects Also the so called narrow band concept 6 is employed to reduce computation time where only the vertices near the ame surface is updated and reinitialized 3 Discretization Solving the level set equation involves the discretizations of a temporal term and several spatial terms With physics in mind the convection term and di usion terms must be discretized using di erent schemes 31 Discretization schemes Here a second order EN0 7 is used to discretize the spatial gradients involved in the terms of VGY and lV l Central di erences are used to discretize each term of Dtle l A rst order forward Euler scheme is used for the temporal term Thus equation 8 can be discretized as follows Gm 7 C3 3 TL Mmaxsgn0v magi 0V 11 At n n z D 2 z Dtb wc myk D204 D2736 ifma39thz k 7 ugerteaw minwzm 7 ugerteam TL TL TL TL maxv2 jk vvertez7 0C m1nvijk vvmezv 0D maXw Lk 7 103mm 0E minw fjk 7 103mm 0Fl7 where A7 B7 C7 D7 E7 F are the terms resulting from the 2nd order ENO scheme as A1l A ngfk min modDfg 7 szizh 12 A 7 B D2716 1 min modDZJff 7 szizh 13 A C Dg k 3 min modD1quot 7 Digit 14 A D D3116 and min modDZJfHJ7 Digit 15 A E ngffk 2117 min modD fz7 BLTEa 16 F i Dz A2h i zz ZZ A i j k 2 m1andDm k 7Di jk 7 where A117A1T7Ay1f7Ay1d7A21b7A21t are the distances between vertices 239 and 239 7 17239 and 1 17 j and j 7 17 j and j 17 k and k 7 17 k and k 17 respectively7 the switch function min mod is given by mmmodmvy 21mm min rvt lyl if any gt 0 otherwise and the di erencing operators in m coordinate direction are de ned by Giggle G2 71jk Di 18 2k Asp 7 62 1 39k i 62 39k Dr 7 7 19 2k Ax 7 D i G 1jkA11 G jkA1l A1Tgt Gi71jkA1T 20 i j k 7 1A2 Am A A 2 7 5 1 17 11 1 Diz Gi71jkA11 A5321 Gi72jkA1l G jkA2l 7 21 7 19 A ZA1z Amy 7 Am Ax22Axll G 1jkA1T AMT Gi2jkA 39E1T GijkA E2T Dzz i lAiTA1T AMT AMT A2T2A1T 23316 7 22 711 11 72 2 y7y 711711 yy 272 7272 zz Other terms such as DiijmWDi j mDi j wDi j k 7Di j k 7Di j k 7Di j k Di j k 7Di j k can be de ned similarly The gradient terms VJF7 V are di erenced as follows V minA7 02 maXB702 minC397 02 maXD7 02 23 minE7 02 maxF7 027 V minB7 02 maXA02 minD7 02 maXC397 02 24 minF02 maXE02 32 Re initialization The re initialization of the level set function is required to reduce the numerical errors caused by steepening and attening gradient e ects The condition lV l 1 is applied to the level sets for 615 75 a 0 The algorithm proposed by Sussman 8 is used to ensure that this condition is satis ed by the level sets for 615 75 a 0 In the algorithm a new level set functin C31 575 is constructed in such a way that its zero level set is the same as the level set function 05t at the ame front while it is set to be the signed normal distance to the ame front away from the zero level set Speci cally i t is de ned as 6315775 3 o1lVC 1i3tl7 25 where SG0 is the sign function designed as slt ogt lt26 1 gft 52 where 5 is a small number And at t 0 1f0 60f0 4 Test Problems and expected results A simple sphere growth77 problem 4 is considered to test the numerical schemes used for the level set method In this experiment the sphere with an initial radius of 10mm grows radially at a prescribed speed of 39mm 10m3 And no ow convection is considered So the level set equation reduces to BC 5 SmoleGl 27 To measure the shape change of the sphere during its grow an e ective radius error parameterRTm is de ned as 7 2 ers ZltT gthewy 7 where n is the total number of piercing points of the spheres surface with cell edges 7 is the calculated radius and rthemy is the theoretical radius of the growing sphere 41 Computational setup The reduced equation is discretized by the schemes described above To test the e ects of grid size on the scheme prediction accuracy three di erent meshes with average cell sizes of 62 31 and 23mm are used The initial sphere radius of 10mm speci es both the initial and boundary conditions for the calculations 42 Expected results Since only sphere growth is considered here the calculated surface would be a nearly perfect sphere surface at any later time if the schemes were accurate Figure 2 shows the comparison between the calculated radius and the theoretical radius Only slight di erences are seen in the radius As the sphere grows the error increases But the relative error BTW rthemy is still relatively small as can be seen in Figure 3 The predicted radius for the di erent grids are shown in Figure 4 For the test case the numerical schemes are not sensitive to the grid size Radlus ant me ms Figure 2 Comparison ofthe govnh ofsphae radius usmg the level see mezhod and th Ihzoretlcal ndms for nprcscnbad speed of SW 10 ms Man 00257 oom some R ml m m DODS on as H vs 2 25 0 35 Txme IRS Figure 3 Relanve error ofpredmted sphere gowth Themeuca rad us Average cell sue 2 3mm Radius grawlh N u 5 Tune m5 Fxgure 4 Effects ofdnffzrem mesh sues on numeneal predmnon aeemey References N3 03 4 U G 1 00 Peters N Turbulent Combustion Cambridge university press 2000 Markstein G Nonsteady ame propagation Pergamon PressOXford 1964 Williams FA The Mathematics of Combustion SIAM Philadelphia PA pp97 1311985 Tan Z and Reitz RD An ignition and combustion model based on the level set method for spark ignition engine multidimensional modeling Combustion and Flame 145 pp1 152006 Osher S and FedkiW R Level Set Methods and Dynamic lmplicit Surfaces Springer Verlag New York Inc 2003 Chopp DL Computing minimal surfaces via level set curvature ow Journal of Computa tional Physics 106 pp77 911993 Harten A Engquist B Osher S and Chakravarthy S Uniformly high order accurate 9 anqu 39rm39nn millatory nheme 2 Journal of Computational Physics 71 pp231 3031987 Sussman M Smereka and S Osher A level set approach for computing solutions to incom pressible 2 phase ow Journal of Computational Physics 114 pp146 1591994 Von Neumann Analysis for Higher Order Methods 1 Introduction Von Neumann analysis is a widely used method to study how an initial wave is propagated with certain numerical schemes for a linear wave equation or heat equation It is the most straightforward analysis method although only limited to linear model equations With the new development of compact higher order method such as discontinuous Galerkin method spectral volume method spectral difference method etc To analyze the properties of this type of schemes for 1D case or 2D quadrilateral or triangular meshes it is necessary to extend this traditional analysis method to these higher order schemes The current presentation focuses on linear wave equation starting with the simplest case the 1D finite difference upwind scheme then extending to 1D higher order scheme analysis 2D quadrilateral mesh and finally 2D triangular mesh All the analysis here still deals with uniform mesh 2 Basic information about Von Neumann analysis 1D upwind scheme 1D wave equation u aux 0 The exact solution isutx fx at with u0 x If the initial condition is a harmonic wave 6 then the exact solution is u t x elklxiatl For wave propagation properties we need to look at the derivative of the solution and the derivative of t the exact solution isa z39kau 1 Now replace the spatial derivative with upwind scheme Q a 311 a Assuming agt0 and substitute in the initial wave at 6x Ax 1 6711 a 3 u aTul ZESIIICACZ coskAx 1ul Then define llnondimensional wave numberquot K kAx t 3M i sinKicosK lul Define the expression inside the square parenthesis as llFourier footprint r of the upwind scheme ie KW K sinK icosK 1 And compare with the exact derivative at ith grid point 6 l l39kaul 411911 at Ax Easy to see rem K K So what kind of information about the scheme can we get from Fourier footprint A The order of accuracy If for a certain scheme the Fourier footprint is rK and rK K 0rp the scheme is pth order accurate Eg for the upwind scheme rum K K 0K2 then upwind scheme is firstorder accurate What if the expression of r K is complicated and hard to do the Taylor expansion 15 K ChooseK K which is small enough evaluate rK K and r3 3 then W L B Dispersion property r 22F 310g2 9le Ixle WI The real part of the Fourier footprint represents the phase the wave propagated with the exact propagation speed K In this upwind case there is a phase lag between the wave propagation and the exact value C Dissipation Property The imaginary part of the Fourier footprint denotes the dissipation property of the wave propagated if the negative part is positive it will cause the amplitude of the wave to grow exponentially and finally make the solution blow up and therefore for the stability of the scheme the imaginary part of the Fourier footprint should normally be negative For the upwind scheme it is easy to see it is stable Dispersion Exact Finite diff 1st order Re0mega 0 05 1 15 2 25 3 K Diffusion Exact Finite diff 1st order m0mega D CFL limit Plot the number ir K in a complex coordinate we will have a closed loop and in the upwind case it is just a circle with the center at 1 O and the radius 1 First for the scheme to be stable the loop should normally be completely on the left hand side of the complex plane Although weak instability may be consolidated by certain time integration scheme the CFL limit for this to happen is often very stringent and therefore it is desirable to have the whole loop on the left hand side of the plane ie we always have negative dissipation Secondly since irKCFL must lie inside the stability region of the time integration scheme used The stability region can be derived from absolute stability theory of ODEs It is desirable to have a smallersized loop which can provide higher CFL limit for the same time integration scheme imwr 73 r2 5 r2 4 5 4 yo 5 o Rewr 3 Extension to 1D compact higher order schemes For compact higher order methods such DG SV and SD a common feature is that each cell carries more than one piece of information and these pieces of information carried by each cell is called degrees of freedom DOFs For DG the DOFs can be chosen in different ways resulting a family of methods for SV the DOFs are chosen as average solutions at sub cells and for SD the DOFs are chosen as the solutions at certain points inside the cell To analyze these schemes all the DOFs of a certain cell form a vector Then we still plug in the spatial discretization scheme to substitute the spatial derivative and after some simplification it will be in the form 6 A 417 497 lie The initial condition is still given as a harmonic wave e and then 41671 147 4147445 ArleilK 147 147161K 139 139 With a 71 l a a a a ltAe m44emgtw WW And then the eigenvalues of matrix GK is the Fourier footprint of the scheme in question Clearly there are multiple eigenvalues but only one of them denotes the principal mode the other modes are substantially damped which means in analysis they can be ignored Moreover it is obvious that GK and GK 2mr has the same eigenvalue so actually all the modes can be obtained by translating the principal mode by 21171 This means in terms of stability if the principal mode is stable there is no need to worry about the stability of other modes In practice only principal mode is studied To pick out the principal mode it is only needed to look for the eigenvalue that approximates the exact value when K is small ie rp K K gt 0 asK gt 0 Once the principal eigenvalue is chosen one can get the same kind of properties The order of accuracy dispersion dissipation properties and CFL limit the same way as is done for the upwind case Below are examples for several 4 h order schemes a Extenxinn m 2n quadrihleml mesh 2D wave equation u acosa sin aVu 0 and the initial condition is given as a plane wave u 0 x y elklxcosgwsmgl Generally the initial wave direction and the wave propagating direction is not the same but for the simplicity of analysis we always lett9 0 Then again we have the property utymt aiKu ij 1 i 1j L7 i1j ij 1 For 2D quad meshes the stencil include five cells so after plug in the spatial scheme it should be in the form 6M 6 a a a a s s E AHJMHJ Awiluwil Amuw Ant Hm A71uz1 Plug in the initial condition it can also be written in the form 6M 6 iiGK w The next step is just to obtain the eigenvalue of GK choose the principal mode and analyze the all the properties of the scheme just as the 1D case 5 2D Triangular mesh analysis The most difficult obstacle in analyzing triangular mesh is that it is hard to plug in the initial solution and relate all the DOFs in the whole stencil to the DOFs of the particular cell Therefore to do this the basic building block contains 2 triangles such that this block can be repeated and form the whole computational domain 51 51 1quot 22 1 I a T Ml M171 Ltle Ml um MM Then plug in the spatial scheme a a a a a a 6 EAzrluzrl 147171uz71 Awuw 14744un 1471uz1 And the expression can be manipulated to be like after plug in initial wave After this everything is the same as the first order case The following are examples of resulted dispersion and dissipation plot of 3rd order spectral volume schemes with different partition TE 03 ca 07 M 06 06 05 05 04 D4 03 03 02 02 M 3 M o 02 04 06 05 1 o 02 0A 05 us 03 US 07 07 05 5 as as 04 04 a 0393 oz 02 M SV3L a1 SV4P o 02 04 06 08 1 o 02 04 06 as 1 Acknowledgement The material of this presentation mainly comes from the lecture by Kris Van den Abeele during his Visit to ISU in Jan 7 Mar 2007 more details and published papers can be found at his website httpmechvubacbethermodynamicsmembersKris20Van20den20Abeelehtm


Buy Material

Are you sure you want to buy this material for

25 Karma

Buy Material

BOOM! Enjoy Your Free Notes!

We've added these Notes to your profile, click here to view them now.


You're already Subscribed!

Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'

Why people love StudySoup

Steve Martinelli UC Los Angeles

"There's no way I would have passed my Organic Chemistry class this semester without the notes and study guides I got from StudySoup."

Anthony Lee UC Santa Barbara

"I bought an awesome study guide, which helped me get an A in my Math 34B class this quarter!"

Jim McGreen Ohio University

"Knowing I can count on the Elite Notetaker in my class allows me to focus on what the professor is saying instead of just scribbling notes the whole time and falling behind."


"Their 'Elite Notetakers' are making over $1,200/month in sales by creating high quality content that helps their classmates in a time of need."

Become an Elite Notetaker and start selling your notes online!

Refund Policy


All subscriptions to StudySoup are paid in full at the time of subscribing. To change your credit card information or to cancel your subscription, go to "Edit Settings". All credit card information will be available there. If you should decide to cancel your subscription, it will continue to be valid until the next payment period, as all payments for the current period were made in advance. For special circumstances, please email


StudySoup has more than 1 million course-specific study resources to help students study smarter. If you’re having trouble finding what you’re looking for, our customer support team can help you find what you need! Feel free to contact them here:

Recurring Subscriptions: If you have canceled your recurring subscription on the day of renewal and have not downloaded any documents, you may request a refund by submitting an email to

Satisfaction Guarantee: If you’re not satisfied with your subscription, you can contact us for further help. Contact must be made within 3 business days of your subscription purchase and your refund request will be subject for review.

Please Note: Refunds can never be provided more than 30 days after the initial purchase date regardless of your activity on the site.