Notes Shallow water � Simplified linear analysis before had dispersion relation g c = k tanh kH • For shallow water, kH is small (that is, wave lengths are comparable to depth) • Approximate tanh(x)=x for small x: c = gH � Now wave speed is independent of wave number, but dependent on depth • Waves slow down as they approach the beach cs533d-term1-2005 1 cs533d-term1-2005 2 What does this mean? PDE’s � We see the effect of the bottom � Saving grace: wave speed independent of k means we can solve as a 2D PDE • Submerged objects (H decreased) show up � We � ll derive these “shallow water equations” as places where surface waves pile up on • When we linearize, we � ll get same wave speed each other � Going to PDE � s also let � s us handle non-square domains, • Waves pile up on each other (eventually changing boundaries • The beach, puddles, … should break) at the beach • Objects sticking out of the water (piers, walls, …) with the right • Waves refract to be parallel to the beach reflections, diffraction, … • Dropping objects in the water � We can � t use Fourier analysis cs533d-term1-2005 3 cs533d-term1-2005 4
Kinematic assumptions Conservation of mass � We � ll assume as before water surface is a height field y=h(x,z,t) � Integrate over a column of water with cross- � Water bottom is y=-H(x,z,t) section dA and height h+H � Assume water is shallow (H is smaller than wave lengths) and calm • Total mass is � (h+H)dA (h is much smaller than H) • Mass flux around cross-section is • For graphics, can be fairly forgiving about violating this… � On top of this, assume velocity field doesn � t vary much in the y � (h+H)(u,w) direction � Write down the conservation law • u=u(x,z,t), w=w(x,z,t) • Good approximation since there isn � t room for velocity to vary much in � In differential form (assuming constant density): � y(otherwise would see disturbances in small length-scale features on ( ) + � � ( h + H ) u ( ) = 0 � t h + H surface) � Also assume pressure gradient is essentially vertical • Good approximation since p=0 on surface, domain is very thin • Note: switched to 2D so u=(u,w) and � =( � / � x, � / � z) cs533d-term1-2005 5 cs533d-term1-2005 6 Pressure Conservation of momentum � Look at y-component of momentum equation: � Total momentum in a column: � � = � � � p � ( ) v t + u � � v + 1 h u h + H � y = � g + � � 2 v u � � H � Assume small velocity variation - so dominant � Momentum flux is due to two things: terms are pressure gradient and gravity: • Transport of material at velocity u with its own � p 1 � y = � g � � ) � � momentum: ( � h u u � H � Boundary condition at water surface is p=0 • And applied force due to pressure. Integrate again, so can solve for p: pressure from bottom to top: ( ) p = � g h � y = � g � � h h ( ) ( ) = � g h � y h + H 2 p � H � H 2 cs533d-term1-2005 7 cs533d-term1-2005 8
Pressure on bottom Shallow Water Equations � Not quite done… If the bottom isn � t flat, � Then conservation of momentum is: there � s pressure exerted partly in the � � � ( h + H ) + � g � t � � ) + � � � � � ( ( h + H ) 2 ( h + H ) 2 � � � g ( h + H ) � H = 0 � u u u horizontal plane � � • Note p=0 at free surface, so no net force there ( ) n = � H x , � 1, � H z � Together with conservation of mass � Normal at bottom is: � � Integrate x and z components of pn over ( ) + � � ( h + H ) u ( ) = 0 � t h + H bottom • (normalization of n and cosine rule for area we have the Shallow Water Equations projection cancel each other out) ( ) � H dA � � g h + H cs533d-term1-2005 9 cs533d-term1-2005 10 Note on conservation form Simplifying Conservation of Mass � Expand the derivatives: � At least if H=constant, this is a system of � ( h + H ) conservation laws + u � � ( h + H ) + ( h + H ) � � u = 0 � t � Without viscosity, “shocks” may develop D ( h + H ) = � ( h + H ) � � u • Discontinuities in solution (need to go to weak integral Dt form of equations) � Label the depth h+H with � : D � • Corresponds to breaking waves - getting steeper and Dt = � � � � u steeper until heightfield assumption breaks down � So water depth gets advected around by velocity, but also changes to take into account divergence cs533d-term1-2005 11 cs533d-term1-2005 12
Simplifying Momentum Interpreting equations � Expand the derivatives: � So velocity is advected around, but also � � ) t + � � � uu � + � g accelerated by gravity pulling down on ( �� u 2 � 2 � � � g � � H = 0 � � � higher water ( ) + �� u � � u + � g � � � � � g � � H = 0 �� u t + � u � t + � u � � � u � For both height and velocity, we have two operations: � Subtract off conservation of mass times velocity: �� u t + �� u � � u + � g � � � � � g � � H = 0 • Advect quantity around (just move it) • Change it according to some spatial � Divide by density and depth: u t + u � � u + g � � � g � H = 0 derivatives � Our numerical scheme will treat these � Note depth minus H is just h: separately: “splitting” u t + u � � u + g � h = 0 Du Dt = � g � h cs533d-term1-2005 13 cs533d-term1-2005 14 Wave equation Deja vu � Only really care about heightfield for � This is the linear wave equation, with wave speed c 2 =gH rendering � Which is exactly what we derived from the � Differentiate height equation in time and dispersion relation before (after linearizing the plug in u equation equations in a different way) � Neglect nonlinear (quadratically small) � But now we have it in a PDE that we have some terms to get confidence in • Can handle varying H, irregular domains… h tt = gH � 2 h cs533d-term1-2005 15 cs533d-term1-2005 16
Shallow water equations Advection � To recap, using eta for depth=h+H: � Let � s discretize just the material derivative (advection equation): D � Dt = � � � � u Dq q t + u � � q = 0 Dt = 0 or Du Dt = � g � h � For a Lagrangian scheme this is trivial: just move the particle that stores q, don � t � We � ll discretize this using “splitting” change the value of q at all • Handle the material derivative first, then the ( ) = q x 0 ,0 ( ) q x ( t ), t right-hand side terms next • Intermediate depth and velocity from just the advection part � For Eulerian schemes it � s not so simple cs533d-term1-2005 17 cs533d-term1-2005 18 Scalar advection in 1D First try: central differences � Let � s simplify even more, to just one � Centred-differences give better accuracy dimension: q t +uq x =0 • More terms cancel in Taylor series � Further assume u=constant � Example: � � � q i � t = � u q i + 1 � q i � 1 � � � And let � s ignore boundary conditions for � � 2 � x now • 2nd order accurate in space • E.g. use a periodic boundary � Eigenvalues are pure imaginary - rules out � True solution just translates q around at Forward Euler and RK2 in time speed u - shouldn � t change shape � But what does the solution look like? cs533d-term1-2005 19 cs533d-term1-2005 20
Testing a pulse A pulse (initial conditions) � We know things have to work out nicely in the limit (second order accurate) • I.e. when the grid is fine enough • What does that mean? -- when the sampled function looks smooth on the grid � In graphics, it � s just redundant to use a grid that fine • we can fill in smooth variations with interpolation later � So we � re always concerned about coarse grids == not very smooth data � Discontinuous pulse is a nice test case cs533d-term1-2005 21 cs533d-term1-2005 22 Centered finite differences Centred finite differences � A few time steps (RK4, small � t) later � A little bit later… • u=1, so pulse should just move right without changing shape cs533d-term1-2005 23 cs533d-term1-2005 24
Centred finite differences What went wrong? � A fair bit later � Lots of ways to interpret this error � Example - phase analysis • Take a look at what happens to a sinusoid wave numerically • The amplitude stays constant (good), but the wave speed depends on wave number (bad) - we get dispersion • So the sinusoids that initially sum up to be a square pulse move at different speeds and separate out � We see the low frequency ones moving faster… • But this analysis won � t help so much in multi- dimensions, variable u… cs533d-term1-2005 25 cs533d-term1-2005 26 Modified PDE’s Interpretation q t + uq x = � u � x 6 � Another way to interpret error - try to account for q xxx it in the physics 6 � Look at truncation error more carefully: � Extra term is “dispersion” q i + 1 = q i + � x � q � x + � x 2 � 2 q � x 2 + � x 3 � 3 q ( ) � x 3 + O � x 4 • Does exactly what phase analysis tells us 2 6 • Behaves a bit like surface tension… q i � 1 = q i � � x � q � x + � x 2 � 2 q � x 2 � � x 3 � 3 q ( ) � x 3 + O � x 4 � We want a numerical method with a different sort of 2 6 truncation error q i + 1 � q i � 1 = � q � x + � x 2 � 3 q ( ) • Any centred scheme ends up giving an odd truncation error --- � x 3 + O � x 3 2 � x dispersion 6 � Let � s look at one-sided schemes � Up to high order error, we numerically solve q t + uq x = � u � x 2 q xxx 6 cs533d-term1-2005 27 cs533d-term1-2005 28
Recommend
More recommend
