Lecture 9Wednesday, April 27, 2005Supplementary Reading: Osher and Fedkiw, §14.5.21 DiscretizationWe are interested in constructing a discretization for the hyperbolic system ofconservation laws~Ut+~F~Ux= 0. (1)Assume that the system (1) has N equations. Then the Jacobian, J =∂~F∂~U,will be and N ×N matrix. Futhermore, we know that J is diagonalizable, sincethe system is hyperbolic. Let us denote the eigenvalues, right eigenvectors, andleft eigenvectors of J as λp, Rp, Lprespectively for p = 1, . . . , N. Recall fromthe previous lecture that we can choose the left and right eigenvectors so thatforR =R1R2. . . RN, L =L1L2...LNthe relation RL = LR = I holds. That is, R and L are chosen to be inve rses.In the previous lecture we looked at a linear, constant coefficient system. Inthat case, the Jacobian was a constant matrix. In general, the Jacobian, andhence its eigensystem, will be spatially varying.As in the scalar case, our discretization is of the form~Uit+~Fi+1/2−~Fi−1/24t= 0.For grid point i0, we need to compute the numerical flux functions at xi0+1/2and xi0−1/2. Let us look in detail at computing Fi0+1/2.The first step is to evaluate the eigensystem at the point xi0+1/2. Sincewe only have~U at the grid points, we obtain~U at the cell walls using thestandard average,~Ui0+1/2= (~Ui0+~Ui0+1)/2. Then for p = 1, . . . , N, wefind the com ponent of the numerical flux function in the p-th characteristicfield. In the p-th characteristic field we have an eigenvalue λp(~Ui0+1/2), left1eigenvector~Lp(~Ui0+1/2), and right eigenvector~Rp(~Ui0+1/2).We put~U valuesand~F (~U) values into the p-th characteristic field by taking the dot productwith the left eigenvector,ui0=~Lp(~Ui0+1/2) ·~Ui0(2)fi0=~Lp(~Ui0+1/2) ·~F (~Ui0) (3)where ui0and fi0are scalars. Once in the characteristic field we perform ascalar version of the conservative ENO scheme obtaining a scalar numerical fluxfunction Fi0+1/2in the scalar field. We take this flux out of the characteristicfield by multiplying with the right eigenvector,~Fpi0+1/2= Fi0+1/2~Rp(~Ui0+1/2) (4)where~Fpi0+1/2is the portion of the numerical flux function~Fi0+1/2from the p-thfield. Once we have evaluated the contribution to the numerical flux functionfrom each field, we get the total numerical flux by summing the contributionsfrom each field,~Fi0+1/2=Xp~Fpi0+1/2(5)completing the evaluation of our numerical flux function at the point xi0+1/2.2 Shallow Water EquationsThe shallow water equations are given byhhut+huhu2+12gh2x= 0. (6)where h is the height of the water, and u is the velocity. The first equationis the equation for conservation of mass, and the second is the equation forconservation of momentum. In order to discretize the system using the proce-dure we described above, we must first find the Jacobian and its eigensystemanalytically.In computing the Jacobian, it is very important to reme mber that we takethe conserved variables (in this case h and hu) to be the independent variables.To make this fact more apparent, we can rewrite the equations asu1u2t+u2u22u−11+12gu21x= 0. (7)2and then define h = u1, u = u2u−11. Below we compute the Jacobian matrix.J =∂~F∂~U= ∂~F1∂u1∂~F1∂u2∂~F2∂u1∂~F2∂u2!= ∂(hu)∂h∂(hu)∂(hu)∂∂h(hu)2h−1+12gh2∂∂(hu)(hu)2h−1+12gh2!=0 1−(hu)2h−2+ gh 2 (hu) h−1=0 1−u2+ gh 2uNote: if you find the treatment of h and hu as independent variables in the abovecomputation confusing, you may prefer rewrite the system as in (7), computethe Jacobian in terms of u1and u2, and then substitute for h and u at the end.Next we find the eigensystem for the Jacobian. We havedet(λI − J) =λ −1−u2+ gh λ − 2u= λ2− 2uλ + u2− gh⇒ λ =2u ±p4u2− 4u2+ 4gh2= u ±pghNext we find the right eigenvectors.Jab=u ±pghab⇒b−u2+ gha + 2ub=u ±pghabTherefore, we haveR1=1u +√gh, R2=1u −√ghThenR =R1R2=1 1u +√gh u −√gh3HenceL = R−1=1−2√ghu −√gh −1−u −√gh 1orL1=−u2√gh+12,12√gh, L2=u2√gh+12, −12√gh.3 Compressible FlowThe inviscid Euler equations for one phase compress ible flow in the absence ofchemical reactions in one spatial dimension are~Ut+~F (~U)x= 0 (8)which can be written in detail asρρuEt+ρuρu2+ p(E + p)ux(9)where ρ is the density, u are the velocities, E is the total energy per unit volume,and p is the pressure. The total energy is the sum of the internal energy andthe kinetic energy,E = ρe + ρ(u2)/2 (10)where e is the internal energy per unit mass.3.1 Ideal Gas E quation of StateFor an ideal gasp = ρRTwhere R = Ru/M is the specific gas constant with Ru≈ 8.31451J/(molK) theuniversal gas constant and M the molecular weight of the gas. Also valid foran ideal gas iscp− cv= Rwhere cpis the specific heat at constant pressure and cvis the specific heat atconstant volume. Gamma is the ratio of specific heats,γ = cp/cv.For an ideal gas, one can writede = cvdT (11)4and as suming that cvdoes not depend on temperature (calorically perfect gas),integration yieldse = eo+ cvT (12)where eois not uniquely determined, and one could choose any value for e at0K. We take e0= 0 arbitrarily for simplicity.Note thatp = ρRT =Rcvρe =cp− cvcvρe = ρ(γ − 1)e = (γ − 1)ρeSo, our equation of state isp = (γ − 1)ρe,orp = (γ − 1)E