CS205b/CME306Lecture 131 Systems of Conservation LawsSupplementary Reading: Osher and Fedkiw, §14.5.1, §14.5.2Consider a simple hyperbolic system of N equationsφt+ [f(φ)]x= 0. (1)1.1 DiscretizationAssume that the system (1) has N equations. Then the Jacobian, J =∂f∂φ, will be and N × Nmatrix. Futhermore, we know that J is diagonalizable, since the system is hyperbolic. Let usdenote the eigenvalues, r ight eigenvectors, and left eigenvectors of J as λp, Rp, Lprespectively forp = 1, . . . , N. Recall f rom the previous lecture that we can choose the lef t and right eigenvectorsso that forR =R1R2··· RN, L =L1L2...LNthe relation RL = LR = I holds. That is, R and L are chosen to be inverses.In the previous lecture we looked at a linear, constant coefficient system. In that case, theJacobian was a constant matrix. In general, the Jacobian, and hence its eigensystem, will bespatially varying.As in the scalar case, our discretization is of the form(φi)t+fi+1/2− fi−1/2∆t= 0.For grid point i0, we need to compute the numerical flux functions at xi0+1/2and xi0−1/2. Let uslook in detail at computing Fi0+1/2.The first step is to evaluate the eigensystem at the point xi0+1/2. Since we only have φ at the gridpoints, we obtain φ at the cell walls using the standard average, φi0+1/2= (φi0+ φi0+1)/2. Thenfor p = 1, . . . , N, we find the component of the numerical flux function in the p-th characteristicfield. In the p-th characteristic field we have an eigenvalue λp(φi0+1/2), left eigenvector Lp(φi0+1/2),and right eigenvector Rp(φi0+1/2).We put φ values and f(φ) values into the p-th characteristic fieldby taking the dot product with the left eigenvector,ui= Lp(φi0+1/2) · φi1fi= Lp(φi0+1/2) · f(φi)where uiand fiare scalars. Once in the characteristic field we perform a scalar version of theconservative ENO scheme obtaining a scalar numerical flux function Fpi0+1/2in the scalar field. Wetake this flux out of the ch aracteristic field by multiplying with the right eigenvector,Fpi0+1/2= Fpi0+1/2Rp(φi0+1/2)where Fpi0+1/2is the portion of the numerical flux function Fi0+1/2from the p-th field. Once wehave evaluated the contribution to the numerical flux function from each fi eld, we get the totalnumerical flux by summing the contributions from each field,Fi0+1/2=XpFpi0+1/2completing 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. (2)where h is the h eight of the water, and u is the velocity. The first equation is the equation forconservation of mass, and the second is the equation for conservation of momentum. In order todiscretize the system using the procedur e we described above, we must first find the Jacobian andits eigensystem analytically.In computing the Jacobian, it is very important to remember that we take the conserved vari-ables (in this case h and hu) to be th e independent variables. To make this fact more apparent,we can rewrite the equations asu1u2t+u2u22u−11+12gu21x= 0. (3)and then define h = u1, u = u2u−11. Below we compute th e Jacobian matrix.J =∂f∂φ= ∂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2Note: if you find the treatment of h and hu as independent variables in the above computationconfusing, you may prefer rewrite the system as in (3), compute the Jacobian in terms of u1andu2, and then substitute f or 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+ 4g h2= 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HenceL = R−1=1−2√ghu −√gh −1−u −√gh 1orL1=−u2√gh+1212√gh, L2=u2√gh+12−12√gh.3 Compressible FlowThe inviscid Euler equations for one phase compressible flow in the absen ce of chemical reactionsin one spatial dimension areφt+ f(φ)x= 0which can be written in detail asρρuEt+ρuρu2+ p(E + p)ux= 03where ρ is the density, u are the velocities, E is the total energy per unit volume, and p is thepressure. Th e total energy is the sum of the internal energy and the kinetic energy,E = ρe + ρ(u2)/2where e is the internal energy per unit mass.3.1 Ideal Gas Equation of StateFor an ideal gasp = ρRTwhere R = Ru/M is the specific gas constant with Ru≈ 8.31451J/(molK) the universal gasconstant and M the molecular weight of the gas. Also valid for an ideal gas iscp− cv= Rwhere cpis the specific heat at constant pressure and cvis the specific h eat at constant volume.Gamma is the ratio of specific heats,γ = cp/cv.For an ideal gas, one can writede = cvdTand assuming that cvdoes not depend on temperature (calorically perfect gas), integration yieldse = eo+ cvTwhere eois not uniquely determined, and one could choose any value for e at 0K. We take e0= 0arbitrarily for simplicity.Note thatp = ρRT =Rcvρe =cp− cvcvρe = ρ(γ − 1)e = (γ − 1)ρe,so our equation of state isp = (γ − 1)ρe,orp = (γ − 1)E