Changes between Version 84 and Version 85 of u/erica/RoeSolver
- Timestamp:
- 05/23/13 16:01:49 (12 years ago)
Legend:
- Unmodified
- Added
- Removed
- Modified
-
u/erica/RoeSolver
v84 v85 10 10 The Roe solver is an approximation means for the numerical flux of the Godunov method, which is derived through linearizing a hyperbolic system of equations. For instance, the Euler equations in conservative form are written 11 11 12 {{{#!Latex 13 \vec{U}_t + \vec{F}(\vec{U})_x = \vec{0} 14 }}} 12 [[latex($\vec{U}_t + \vec{F}(\vec{U})_x = \vec{0}$)]] 15 13 16 14 which using the chain rule is identical to 17 15 18 {{{#!Latex 19 \vec{U}_t + \frac{\partial \vec{F}} {\partial \vec{U}} ~\vec{U}_x = \vec{0} 20 }}} 16 [[latex($\vec{U}_t + \frac{\partial \vec{F}} {\partial \vec{U}} ~\vec{U}_x = \vec{0}$)]] 21 17 22 18 If we assume that 23 19 24 {{{#!Latex 25 \hat{A} = \frac{\partial \vec{F}}{\partial \vec{U}} 26 }}} 20 [[latex($\hat{A} = \frac{\partial \vec{F}}{\partial \vec{U}} $)]] 27 21 where, A-hat is a Jacobian matrix of averaged/constant values, we can derive an expression for the numerical flux in terms of 1) wave strengths (alpha), and the 2) eigenvalues (lambda) and 3) right eigenvectors (K) of the 'averaged' Jacobian: 28 22 29 {{{#!Latex 30 \vec{F}_{i+1/2}= \frac{1}{2}(\vec{F}_L + \vec{F}_R) - \frac{1}{2}\Sigma_{i=1,m} \tilde{\alpha}_i \tilde{\lambda}_i \tilde{\vec{K^{(i)}}} 31 }}} 32 23 [[latex($\vec{F}_{i+1/2}= \frac{1}{2}(\vec{F}_L + \vec{F}_R) - \frac{1}{2}\Sigma_{i=1,m} \tilde{\alpha}_i \tilde{\lambda}_i \tilde{\vec{K^{(i)}}}$)]] 33 24 where 34 25 35 {{{#!Latex 36 \vec{F}_k = \vec{F}(\vec{U}_k) 37 }}} 38 26 [[latex($\vec{F}_k = \vec{F}(\vec{U}_k)$)]] 39 27 So the goal is to compute the wave speeds and associated eigenvalues and eigenvectors of the Jacobian matrix. There are 2 methods by which we can do this: 1) The 'Roe' approach, which (non-trivially) constructs an averaged Jacobian directly that must satisfy conditions such as hyperbolicity and conservation, and 2), the newer 'Roe-Pike' approach, which avoids solving for the Jacobian and instead develops algebraic expressions for the sought quantities based on averages of the initial data. It is the 2nd, more widely used approach, that we will explore here. 40 28 … … 45 33 From the equation above of the Euler equations in Jacobian form, we can derive the equations for wave speeds: 46 34 47 {{{#!Latex 48 \triangle \vec{U} = \vec{U}_R - \vec{U}_L = \Sigma _{i=1,m} \hat{\alpha}_i(\vec{W}) \hat{\vec{K}}^{(i)}(\vec{W}) 49 }}} 35 [[latex($\triangle \vec{U} = \vec{U}_R - \vec{U}_L = \Sigma _{i=1,m} \hat{\alpha}_i(\vec{W}) \hat{\vec{K}}^{(i)}(\vec{W})$)]] 50 36 51 37 The eigenvalues and vectors are derived straightforwardly from the Jacobian, using typical methods. … … 57 43 For the x-split, 3-dimensional Euler equations, we have the following eigenvalues and vectors: 58 44 59 {{{#!Latex 60 \lambda_1 = u-a, ~\lambda_2 = \lambda_3 = \lambda_4 = u, ~\lambda_5 = u+a 61 }}} 62 45 [[latex($\lambda_1 = u-a, ~\lambda_2 = \lambda_3 = \lambda_4 = u, ~\lambda_5 = u+a $)]] 63 46 and 64 47 65 48 66 {{{#!Latex 67 \vec{K}^1 = <1, u-a, v, w, H-ua>, 68 }}} 69 {{{#!Latex 70 \vec{K}^2 = <1,u,v,w,V^2/2>, 71 }}} 72 {{{#!Latex 73 \vec{K}^3 = <0,0,1,0,v>, 74 }}} 75 {{{#!Latex 76 \vec{K}^4 = <0,0,0,1,w>, 77 }}} 78 {{{#!Latex 79 \vec{K}^5 = <1, u+a, v, w, H+ua> 80 }}} 49 [[latex($\vec{K}^1 = <1, u-a, v, w, H-ua>,$)]] 50 [[latex($\vec{K}^2 = <1,u,v,w,V^2/2>,$)]] 51 [[latex($\vec{K}^3 = <0,0,1,0,v>, $)]] 52 {[[latex($\vec{K}^4 = <0,0,0,1,w>, $)]] 53 {[[latex($\vec{K}^5 = <1, u+a, v, w, H+ua> $)]] 81 54 82 55 where u, v, w, are velocities in x, y, z directions, respectively, a is the local sound speed, given by 83 56 84 {{{#!Latex 85 a^2 = \frac{\gamma P}{\rho} 86 }}} 57 [[latex($a^2 = \frac{\gamma P}{\rho}$)]] 87 58 88 59 H is the enthalpy given by 89 60 90 {{{#!Latex 91 H = \frac{E+p}{\rho} 92 }}} 93 61 [[latex($H = \frac{E+p}{\rho}$)]] 94 62 E is the total energy per unit volume, 95 63 96 {{{#!Latex 97 E = \frac{1}{2}\rho \vec{V}^2 + \rho e 98 }}} 99 64 [[latex($E = \frac{1}{2}\rho \vec{V}^2 + \rho e$)]] 100 65 with 101 66 102 {{{#!Latex 103 \vec{V}^2 = u^2 + v^2 + w^2 104 }}} 67 [[latex($\vec{V}^2 = u^2 + v^2 + w^2$)]] 105 68 The Euler equations are a set of non-linear PDE's, which comprise an eigenvalue problem. The eigenvalues for the equations are functions of the solution to the equations themselves. This means that the waves, which propagate with speeds = to their eigenvalues, distort the solution, and the solution distorts them over space and time. Thus the solution to the Riemann 106 69 107 70 and e being the specific internal energy, which for ideal gases is, 108 71 109 {{{#!Latex 110 e = \frac{p}{(\gamma - 1)\rho} 111 }}} 112 72 [[latex($e = \frac{p}{(\gamma - 1)\rho}$)]] 113 73 Recall that the numerical flux is given by, 114 74 115 {{{#!Latex 116 \vec{F}_{i+1/2}= \frac{1}{2}(\vec{F}_L + \vec{F}_R) - \frac{1}{2}\Sigma_{i=1,m} \tilde{\alpha}_i \tilde{\lambda}_i \tilde{\vec{K^{(i)}}} 117 }}} 118 75 [[latex($\vec{F}_{i+1/2}= \frac{1}{2}(\vec{F}_L + \vec{F}_R) - \frac{1}{2}\Sigma_{i=1,m} \tilde{\alpha}_i \tilde{\lambda}_i \tilde{\vec{K^{(i)}}}$)]] 119 76 where the tilde's denote derived expressions evaluated at the average primitive variables: 120 77 121 {{{#!Latex 122 \tilde{\alpha}_i = \hat{\alpha}_i(\tilde{\vec{W}}), ~\tilde{\lambda}_i = \lambda_i(\tilde{\vec{W}}), ~\tilde{\vec{K}^{(i)}} = \vec{K}^{(i)}(\tilde{\vec{W}}) 123 }}} 78 [[latex($\tilde{\alpha}_i = \hat{\alpha}_i(\tilde{\vec{W}}), ~\tilde{\lambda}_i = \lambda_i(\tilde{\vec{W}}), ~\tilde{\vec{K}^{(i)}} = \vec{K}^{(i)}(\tilde{\vec{W}})$)]] 124 79 125 80 Note that the vector equation for the numerical flux makes up 3 equations for the 3D Euler equations, one for each component of the following vectors: 126 81 127 82 the vector of numerical fluxes in conserved form, 128 {{{#!Latex 129 \vec{F}_{i+1/2} = <\rho_{i+1/2}, ~p_{i+1/2}, ~E_{i+1/2}> 130 }}} 131 83 [[latex($\vec{F}_{i+1/2} = <\rho_{i+1/2}, ~p_{i+1/2}, ~E_{i+1/2}>$)]] 132 84 vector of left-fluxes in conserved form, 133 {{{#!Latex 134 \vec{F}_L = <\rho_L, ~\rho_L u_L, ~0.5\rho_L u_L^2 + \frac{p_L}{\gamma - 1}> 135 }}} 85 [[latex($\vec{F}_L = <\rho_L, ~\rho_L u_L, ~0.5\rho_L u_L^2 + \frac{p_L}{\gamma - 1}>$)]] 136 86 137 87 vector of right-fluxes in conserved form, 138 {{{#!Latex 139 \vec{F}_R = <\rho_R, ~\rho_R u_R, ~0.5\rho_R u_R^2 + \frac{p_R}{\gamma - 1}> 140 }}} 141 88 [[latex($\vec{F}_R = <\rho_R, ~\rho_R u_R, ~0.5\rho_R u_R^2 + \frac{p_R}{\gamma - 1}>$)]] 142 89 and eigenvectors, 143 90 144 {{{#!Latex 145 \vec{K}^{(1)} = <1, ~u-a, ~H-ua> 146 }}} 91 [[latex($\vec{K}^{(1)} = <1, ~u-a, ~H-ua>$)]] 147 92 148 {{{#!Latex 149 \vec{K}^{(2)} = <1, ~u, ~0.5 u^2> 150 }}} 93 [[latex($\vec{K}^{(2)} = <1, ~u, ~0.5 u^2>$)]] 151 94 152 {{{#!Latex 153 \vec{K}^{(3)} = <1, ~u+a, ~H+ua> 154 }}} 95 [[latex($\vec{K}^{(3)} = <1, ~u+a, ~H+ua>$)]] 155 96 156 97 Finally, the derived expressions for the wave-speeds are as follows, 157 98 158 {{{#!Latex 159 \hat{\alpha}_1 = \frac{\triangle p - \hat{\rho} \hat{a} \triangle u}{2 \hat{a} ^2} 160 }}} 99 [[latex($\hat{\alpha}_1 = \frac{\triangle p - \hat{\rho} \hat{a} \triangle u}{2 \hat{a} ^2}$)]] 161 100 162 {{{#!Latex 163 \hat{\alpha}_2 = \triangle \rho - \triangle p / \hat{a}^2 164 }}} 101 [[latex($\hat{\alpha}_2 = \triangle \rho - \triangle p / \hat{a}^2$)]] 165 102 166 {{{#!Latex 167 \hat{\alpha}_3 = \hat{\rho} \triangle v 168 }}} 103 [[latex($\hat{\alpha}_3 = \hat{\rho} \triangle v$)]] 169 104 170 {{{#!Latex 171 \hat{\alpha}_4 = \hat{\rho} \triangle w 172 }}} 105 [[latex($\hat{\alpha}_4 = \hat{\rho} \triangle w$)]] 173 106 174 {{{#!Latex 175 \hat{\alpha}_5 = \frac{\triangle p + \hat{\rho} \hat{a} \triangle u}{2 \hat{a} ^2} 176 }}} 107 [[latex($\hat{\alpha}_5 = \frac{\triangle p + \hat{\rho} \hat{a} \triangle u}{2 \hat{a} ^2}$)]] 177 108 178 109 and the corresponding expressions for rho, u, v, w, H, and a that we would plug into these wave speeds and the above expressions for the eigenvalues and eigenvectors are listed in Toro, eqn. 11.118. These in turn would be used in the numerical flux for the Godunov method. … … 180 111 Note in the above formalism, the following definition is used, 181 112 182 {{{#!Latex 183 \triangle r = r_R - r_L 184 }}} 113 [[latex($\triangle r = r_R - r_L$)]] 185 114 186 115