Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 77
Текст из файла (страница 77)
As shown in Fig. 13.11, the value of ρϕ atthe temporal element face is set equal to the value at the centroid of the upwindelement to giveðqC /C ÞtþDt=2 ¼ ðqC /C Þt and ðqC /C ÞtDt=2 ¼ ðqC /C ÞtDtð13:60Þ13.3The Finite Volume Approach509t+ tt+ t/2ttL( )tt/2ttFig. 13.11 First order transient upwind interpolation profile resulting in the implicit first ordertransient Euler schemeUsing Eq. (13.60), Eq. (13.57) becomes ðqC /C Þt ðqC /C ÞtDtVC þ L /tC ¼ 0Dtð13:61Þwhich is the first order implicit Euler scheme.
The scheme is linearized as follows:qC VCDtq VCFluxC ¼ CDtFluxV ¼ 0FluxC ¼13.3.2.1ð13:62ÞNumerical DiffusionAs this is a first order scheme, it is expected, based on the knowledge gained fromconvection schemes, to produce numerical diffusion. Its value can be determined bytrying to recover the original governing equation using a Taylor series expansionaround time t. The value of ðq/ÞtDt can be expressed asðq/ÞtDt @ ðq/Þ@ 2 ðq/Þ Dt2þ O Dt3¼ ðq/Þ Dt þ@t t@t2 t 2tð13:63Þ51013 Temporal Discretization: The Transient Termwhich can be rearranged into 2 ðq/Þt ðq/ÞtDt @ ðq/ÞDt @ ðq/Þ¼ O Dt22Dt@t t2@t|fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl}tNumericaldiffusiontermð13:64ÞSubstituting Eq. (13.64) into the discretized equation gives 2 @ ðq/Þ1 tDt @ ðq/Þþ L /C ¼þ O Dt22@t t VC2@t|fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl}tNumericaldiffusiontermð13:65ÞIn effect a numerical diffusion term has been added to the equation that scales withthe time step in a similar fashion to the upwind scheme for the advection term.
Sowhile the scheme is unconditionally stable, the solution it yields is really a stationary solution for large time steps.13.3.3First Order Explicit Euler SchemeThe transient first order explicit Euler scheme is obtained by using a first-order“downwind” interpolation profile [4, 18]. As shown in Fig.
13.12, the value of ρϕ att+ tt+ t/2ttL( )tt/2ttFig. 13.12 First order transient downwind interpolation profile resulting in the explicit first ordertransient Euler scheme13.3The Finite Volume Approach511the temporal element face is set equal to the value at the downwind elementcentroid, yieldingðqC /C ÞtþDt=2 ¼ ðqC /C ÞtþDtand ðqC /C ÞtDt=2 ¼ ðqC /C Þtð13:66ÞUsing Eq. (13.66), Eq.
(13.57) becomes ðqC /C ÞtþDt ðqC /C ÞtVC þ L /tC ¼ 0Dtð13:67Þwhich is the first order explicit Euler scheme. The scheme is linearized as follows:qC VCDtq VCFluxC ¼ CDtFluxV ¼ 0FluxC ¼ð13:68ÞNote that now the new time is at t þ Dt and that the spatial operator of Eq. (13.67)has to be evaluated at time t. Thus, it is possible to evaluate the right hand sidecompletely and find the value of ρϕ at time t þ Dt without the need to solve a set oflinear algebraic equations. This is the explicit scheme and corresponds to theassumption that ρϕ prevails over the entire time step.13.3.3.1Numerical Anti-DiffusionAgain performing a simple Taylor expansion around time t yields @ ðq/Þ@ 2 ðq/Þ Dt2þ O Dt3Dtþ2@t t@tt 2ð13:69Þ 2 ðq/ÞtþDt ðq/Þt @ ðq/ÞDt @ ðq/Þ¼þþ O Dt2Dt@t t2@t2 tð13:70Þðq/ÞtþDt ¼ ðq/Þt þwhich can be rewritten asSubstitution into Eq.
(13.67) gives 2 @ ðq/Þ1 tDt @ ðq/Þþ L /C ¼ þ O Dt2@t t VC2@t2 t|fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl}Numericalantidiffusiontermð13:71Þ51213 Temporal Discretization: The Transient TermWhere now the second order differential term has a negative sign, akin to a negativediffusion or anti-diffusion, with compression effects on profiles, very similar to theDownwind scheme in advection. Again the anti-diffusion term scales with the timestep. When used in combination with an upwind convection scheme and a Courantnumber of 1, it can be shown that the numerical diffusion of the advection scheme andthe numerical anti-diffusion of the explicit Euler scheme for a CFLconv equals to 1 areof equal magnitudes and of opposite signs.
Thus they cancel each other producingnearly an exact solution. Nonetheless this is not practical as ensuring a CFLconv of 1on anything but simple one dimensional grids is not an option for real problems.A related issue to the anti-diffusion behavior is numerical instabilities, whichincreases with increasing Dt placing a very strong restriction on the time step.
Thiscan be evaluated by applying the negative neighboring coefficient rule.13.3.4Second Order Transient Euler SchemesSimilar to advection schemes, second order transient schemes can be constructedwith a linear interpolation profile. The choice could be a symmetric profile (centraldifference) yielding the Crank-Nicolson (CN) scheme [2], or an upwind one (secondorder upwind scheme) [4, 19, 20] resulting in the Adams-Moulton scheme [15, 16],an implicit scheme also known as the Second Order Upwind Euler (SOUE).13.3.5Crank-Nicholson (Central Difference Profile)With the ρϕ computed using linear interpolation between the “Upwind” and“Downwind” nodes, the Crank-Nicholson scheme shown in Fig.
13.13 is obtained.For a uniform time step, this is expressed mathematically as11ðq / ÞtþDt þ ðqC /C Þt2 C C21tDt=2t 1ðqC /C Þ¼ ðqC /C Þ þ ðqC /C ÞtDt22ðqC /C ÞtþDt=2 ¼ð13:72ÞSubstituting in Eq. (13.57), the discretized equation becomes ðqC /C ÞtþDt ðqC /C ÞtDtVC þ L /tC ¼ 02Dtð13:73ÞThe linearization coefficients for the CN scheme can be written asqC VC2DtFluxC ¼ 0q VC /FluxV ¼ C2Dt CFluxC ¼ð13:74Þ13.3The Finite Volume Approach513t+ tt+ t/2Ltt( )tt/2ttFig. 13.13 Second order transient central difference interpolation profile resulting in the transientCN schemeThe stencil shown in Fig. 13.9 indicates that the scheme is explicit with the value atlevel t þ Dt computed explicitly from the values at times t and t Dt.
Thus itsstability is constrained by a CFL limit.Again in a similar fashion to the finite difference formulation, it can be reformulated in a two-step procedure using Eqs. (13.48) and (13.49), i.e., a first orderimplicit Euler step followed by a modified explicit Euler step in the form ofextrapolation.13.3.5.1Numerical AccuracyExpanding ðq/Þ at t þ Dt and t Dt via Taylor expansions around time t yieldsðq/ÞtþDt @ ðq/Þ@ 2 ðq/Þ Dt2 @ 3 ðq/Þ Dt3þþ O Dt4 ð13:75Þ¼ ðq/Þ þDt þ@t @t2 2@t3 6tttt @ ðq/Þ@ 2 ðq/Þ Dt2 @ 3 ðq/Þ Dt3tDttþ O Dt4 ð13:76ÞDt þðq/Þ ¼ ðq/Þ 23@t t@t@tt 2t 6Subtracting Eq.
(13.75) from Eq. (13.76), the following equation is obtained: ðq/ÞtþDt ðq/ÞtDt @ ðq/Þ @ 3 ðq/Þ Dt2¼ O Dt3þ32Dt@t t@tt 6ð13:77Þ51413 Temporal Discretization: The Transient TermSubstitution into Eq. (13.73) gives @ ðq/Þ1 t@ 3 ðq/Þ Dt2þ O Dt3þ L /C ¼ @t t VC@t3 t 6ð13:78Þconfirming that the scheme is second order accurate. The third order derivative is adispersive term that results in instability.13.3.6Second Order Upwind Euler (SOUE) SchemeUsing the second-order “upwind” interpolation profile depicted in Fig. 13.14, theinterface ρϕ values are approximated as31ðqC /C Þt ðqC /C ÞtDt223tDt=2tDt 1ðqC /C Þ¼ ðqC /C Þ ðqC /C Þt2Dt22ðqC /C ÞtþDt=2 ¼ð13:79ÞSubstituting in Eq.
(13.57), the discretized ρϕ field equation is obtained as 3ðqC /C Þt 4ðqC /C ÞtDt þ ðqC /C Þt2DtVC þ L /tC ¼ 02Dtð13:80Þwhich is the implicit second order upwind Euler (SOUE) scheme. In this schemethe values of ρϕ have to be stored for two of the older time steps, with its linearization coefficients given byt+ t/2Lttt/2ttt 2 tFig. 13.14 Second order upwind Euler scheme( )t13.3The Finite Volume Approach5153qC VC2Dt2q VCFluxC ¼ CDtqVC /CFluxV ¼ C2DtFluxC ¼13.3.6.1ð13:81ÞNumerical AccuracyThe scheme is second order as can be shown from a Taylor series evaluation.Expanding ðq/ÞtDt and ðq/Þt2Dt around time t giveðq/ÞtDt ¼ ðq/Þt @ ðq/Þ@ 2 ðq/Þ Dt2 @ 3 ðq/Þ Dt3þ O Dt4 ð13:82ÞDtþ23@t t@t@tt 2t 6ðq/Þt2Dt ¼ ðq/Þt @ ðq/Þ@ 2 ðq/Þ@ 3 ðq/Þ 8Dt32þ O Dt42Dtþ2Dt23@t t@t@ttt 6ð13:83ÞMultiplying Eq.
(13.82) by 4 and subtracting the resulting equation fromEq. (13.83), an expression for the SOUE is obtained as 3ðq/Þt 4ðq/ÞtDt þðq/Þt2Dt @ ðq/Þ @ 3 ðq/Þ Dt2¼ O Dt32Dt@t t@t3 t 3ð13:84ÞCombining Eq. (13.84) and Eq. (13.80), the recovered equation for ρϕ becomes @ ðq/Þ1 t @ 3 ðq/Þ Dt2þ O Dt3þ L /C ¼@t t VC@t3 t 3ð13:85Þwhich has a third order numerical dispersion term but no numerical diffusion.13.3.7Initial Condition for the FV ApproachThe implementation of the finite volume formulation is straight forward except forthe initial time step. As shown in Fig. 13.15, the first temporal element is aboundary element in time, as such it does not have an upwind neighbor.
Rather thevalue at the lower element face is used directly at the face resulting in a gradientthat is half the correct numerical value. This comes about because it is computed astþDt=2and /tCinitial , which are located half athe difference between the values at /Cinitial51613 Temporal Discretization: The Transient Termtt/2tinitialFig. 13.15 Boundary temporal elementtime step ðDt=2Þ apart, while dividing their difference by a full time step ðDtÞleading to a non-negligible initial error.This is easily demonstrated by considering the first temporal element in thediscretized equation of the first order implicit Euler scheme.
Using Eq. (13.57) thediscretized ρϕ field equation is obtained asðqC /C Þtinitial þDt=2 ðqC /C ÞtinitialtþDt=2VC þ L /Cinitial¼0Dtð13:86ÞFor the first temporal element, the upwind interpolation yields a gradient computedas the difference between the ρϕ values at tinitial þ Dt=2 and tinitial divided by Dt.However for the case of a regular element (Fig. 13.16), the gradient is actuallybetween the ρϕ values at tinitial þ 3Dt=2 and tinitial þ Dt=2, divided by Dt. Thettt /2tinitialFig.
13.16 Treatment of initial condition and the virtual element of centroid tinitial13.3The Finite Volume Approach517difference between the two gradients is substantial, and any scheme that starts withthe gradient of Eq. (13.86) will result in a large initial error that will affect thesolution at the following steps. This error can be avoided if a grid similar toFig. 13.16 is adopted. In this case the solution of the finite difference and finitevolume methods will be basically similar, as for a regular grid.Adopting this approach, the upwind values at the faces of the first temporalelement spanning the time interval ½tinitial þ Dt=2; tinitial þ 3Dt=2 are obtained asðqC /C Þtinitial þ3Dt=2 ¼ ðqC /C Þtinitial þDtðqC /C Þtinitial þDt=2 ¼ ðqC /C Þtinitialð13:87ÞSubstituting in Eq. (13.57), the discretized ρϕ field equation becomesðqC /C Þtinitial þDt ðqC /C ÞtinitialVC þ L /tCinitial þDt ¼ 0Dtð13:88Þwhich is similar to the equation obtained for any internal element.Example 2Repeat example 1 using the CN scheme applied via Eq.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
