Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 76
Текст из файла (страница 76)
Using the following approximation:Fig. 13.9 Stencil of theCrank Nicholson SchemeCCC50213 Temporal Discretization: The Transient Term/ / þ /2ð13:40Þthe algebraic equation [Eq. (13.38)] becomes0aC /C þ 0:5@aC /C þ10 aF /F A ¼ bC 0:5@ aC þ 2aC /C þXF NBðCÞ1XAaF /FF NBðC Þð13:41ÞThus, the stability condition becomesaC þ 2aC 0:ð13:42ÞFor the one dimensional transient advection problem displayed in Fig. 13.5,Eq. (13.42) results inDt 2q2q2DxCC VCC DxC DyC¼ ;qC uC DyCm_ eveð13:43Þwhere it has been assumed that the advection term is discretized using the upwindscheme. Using the CFL number for convection defined above, Eq. (13.43) isexpressed asCFLconv 2ð13:44ÞThe larger CFL limitation is pleasing, but the improved accuracy is just moreimportant as it allows for accurate solutions to be achieved without the need toresort to very small time steps, especially that the second order derivative is noweliminated from the error.
More details on accuracy analysis will be presented inlater sections.13.2.5Implementation DetailsThe CN scheme can also be derived by summing the Forward and Backwardtransient Euler schemes [4], as shown next.Forward Euler ! ðqC /C Þt ðqC /C ÞtDtVC ¼ L /tCDtð13:45Þ13.2The Finite Difference ApproachBackward Euler !503 ðqC /C ÞtþDt ðqC /C ÞtVC ¼ L /tCDtð13:46ÞForward Euler + Backward Euler: ðqC /C Þt ðqC /C ÞtDtðq / ÞtþDt ðqC /C ÞtVC þ C CVC ¼ L /tC L /tCDtDt ð13:47Þðq / ÞtþDt ðqC /C ÞtDt! C CVC þ L /tC ¼ 02Dt! Crank Nicolson!This formulation points to a simple implementation of the CN scheme within animplicit scheme framework as a two-step procedure.
In the first step a BackwardEuler formulation is used to implicitly find ðq/Þt fromðqC /C Þt þDt t L /C ¼ ðqC /C ÞtDtVCð13:48Þwhile in the second step the CN value at time step t þ Dt is found explicitly as ðq / Þt ðqC /C ÞtDtðqC /C ÞtþDt ðqC /C ÞtVC ¼ L /tC ¼ C CVCDtDttþDtttDt) ðqC /C Þ ¼ 2ðqC /C Þ ðqC /C Þð13:49ÞIn this derivation it was assumed that the transient time step Dt is divided intotwo equal local time steps ðDtlocal Þ, with Dtlocal equals half the set time step Dt.It is important to note that while the CN scheme is second order accurate, it isstill an explicit scheme, which is constrained by a CFL like condition, as explainedabove.13.2.6Adams-Moulton SchemeThe development of the second order Adams-Moulton scheme [15, 16] requiresexpanding the values of T at t Dt and t 2Dt using Taylor series expansionsaround t, yielding@T ðtÞ@ 2 T ðtÞ 4Dt2þ 2Dt þ@t@t2 2!@T ðtÞ@ 2 T ðtÞ Dt2þ Dt þT ðt DtÞ ¼ T ðtÞ @t@t2 2!T ðt 2DtÞ ¼ T ðtÞ ð13:50Þ50413 Temporal Discretization: The Transient TermThe first derivative is obtained by combining the two equations in such a way thatthe second order derivative is eliminated, resulting in the following equation:@T ðtÞ 3T ðtÞ 4T ðt DtÞ þ T ðt 2DtÞ¼@t2Dtð13:51Þwhich, upon substituting in Eq.
(13.3), yields 3ðqC /C Þt 4ðqC /C ÞtDt þ ðqC /C Þt2DtVC þ L /tC ¼ 02Dtð13:52ÞExpanding the spatial term, the final form of the algebraic equation is obtained asX aF /F ¼ bC aC /C að13:53ÞaC þ aC /C þC /CF NBðC Þwith the coefficients given by3qC VC2Dt2q VCaC ¼ CDtqVC CaC ¼2DtaC ¼ð13:54ÞIt is clear that the aC coefficient has a positive sign implying that an increase in /Cwould lead to a decrease in /C .
This is mitigated by the large aC coefficient, whichhas the right influence. Thus while the scheme is stable, it is not bounded withunphysical oscillations expected in certain circumstances.Example 1The thermal conductivity of a solid sphere of volume 1 m3 is so high that itsresistance to conduction is very small as compared to its resistance to convection heat transfer with the surroundings.
Thus temperature gradientswithin the sphere are negligible and the temperature of the sphere is spatiallyuniform at any instant. The initial temperature of the sphere is Th and that ofthe surroundings is T1 . The density, specific heat, sphere surface area, andconvection heat transfer coefficient with the surroundings are ρ, c, As, andh1 , respectively. Neglecting heat transfer by radiation, the energy equationfor the sphere is given byqcVdT¼ h1 AS ðT T1 ÞdtDefining a dimensionless temperature as13.2The Finite Difference Approach505/¼T T1Th T1the energy equation and initial condition becomed/h1 AS/ and /ð0Þ ¼ 1¼qcVdtFor a value of h1 AS =qcV ¼ 1, compare dimensionless temperature valuesobtained analytically with numerical ones generated using the first orderexplicit scheme, the first order implicit scheme, and the second order CNscheme applied via the two-step procedure at times 0.1, 0.2, and 0.3 using atime step with size of 0.1.SolutionThe governing equation reduces tod/¼ /dtsubject to/ð0Þ ¼ 1By separation of variables and application of the initial condition the analytical solution is found as/ðtÞ ¼ etThus the analytical solution at times 0.1, 0.2, and 0.3 are given by/exact ð0:1Þ ¼ e0:1 ¼ 0:9048/exact ð0:2Þ ¼ e0:2 ¼ 0:8187/exact ð0:3Þ ¼ e0:3 ¼ 0:7408The numerical solution is obtained with V ¼ 1, Lð/n Þ ¼ /n , andL /nþ1 ¼ /nþ1 :The error in the numerical solution is found usingerror ¼ j/numerical /exact jNumerical solution using the first order explicit scheme/tþDt ¼ ð1 DtÞ/t50613 Temporal Discretization: The Transient Term/explicit ð0:1Þ ¼ ð1 0:1Þ/ð0Þ ¼ 0:9 1 ¼ 0:9/explicit ð0:2Þ ¼ ð1 0:1Þ/ð0:1Þ ¼ 0:9 0:9 ¼ 0:819>=>;/explicit ð0:3Þ ¼ ð1 0:1Þ/ð0:2Þ ¼ 0:9 0:81 ¼ 0:72983>< errorexplicit ð0:1Þ ¼ 4:8 10) errorexplicit ð0:2Þ ¼ 8:7 103>:errorexplicit ð0:3Þ ¼ 1:18 102Numerical solution using the first order implicit scheme1/t1 þ Dt91>/implicit ð0:1Þ ¼ð1Þ ¼ 0:9091 >8>>1 þ 0:1>errorimplicit ð0:1Þ ¼ 4:3 103>>>>=<1errorimplicit ð0:2Þ ¼ 7:7 103/implicit ð0:2Þ ¼ð0:9091Þ ¼ 0:8264 )>>1 þ 0:1>>>:>>errorimplicit ð0:3Þ ¼ 1:05 102>1>/implicit ð0:3Þ ¼ð0:8264Þ ¼ 0:7513 ;1 þ 0:1/tþDt ¼Numerical solution using the second order CN schemeIn this case the solution is obtained using Eqs.
(13.48) and (13.49) that arereduced to1/ ðt Þ1 þ Dt=2/CN ðt þ DtÞ ¼ 2/ ðt þ Dt=2Þ /ðtÞ/ ðt þ Dt=2Þ ¼where the total time step Dt has been divided into two equal time steps ofvalue Dt=2. Applying the above equations, the solutions are found as9>>>>>>>>>>/CN ð0:1Þ ¼ 2/ ð0:05Þ /ð0Þ ¼ 0:90476 >>>8>>>errorCN ð0:1Þ ¼ 7:551 1051>>><=/ð0:1Þ ¼ 0:861678 >/ ð0:15Þ ¼1 þ 0:05) errorCN ð0:2Þ ¼ 1:366 104>>>>:>/CN ð0:2Þ ¼ 2/ ð0:15Þ /ð0:1Þ ¼ 0:81859 >>errorCN ð0:3Þ ¼ 1:854 104>>>>1>>/ð0:2Þ ¼ 0:779615 >/ ð0:25Þ ¼>>1 þ 0:05>>>>;/CN ð0:3Þ ¼ 2/ ð0:25Þ /ð0:2Þ ¼ 0:7406/ ð0:05Þ ¼1/ð0Þ ¼ 0:952381 þ 0:0513.313.3The Finite Volume Approach507The Finite Volume ApproachThe Finite Volume approach for the discretization of the transient term is verysimilar to the discretization of the convective term [4], except that the integration iscarried over temporal rather than spatial element (Fig.
13.10).Integration of Eq. (13.3) over the time interval ½t Dt=2; t þ Dt=2 yieldstþDt=2ZtþDt=2Z@ ðqC /C ÞVC dt þ@ttDt=2Lð/C Þdt ¼ 0ð13:55ÞtDt=2|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl}Term ITerm IIWith VC treated as a constant, Term I turned into a difference of face fluxes, andTerm II evaluated as a volume integral using the mid point rule, Eq. (13.55)becomes VC ðqC /C ÞtþDt=2 VC ðqC /C ÞtDt=2 þL /tC Dt ¼ 0ð13:56ÞEquation (13.56) is the semi-discretized transient equation, which can be written inthe more standard form by dividing all terms by the temporal element volume, Dt,leading toElement(temporal domain)(xC,t + t ) =x C ,t +t=2( x ,t ) =C(xt+t2CttL(tCx C ,tt=2,tt) =Ctt+ tCtt2Ctbt tCtFig. 13.10 Element in the transient domaint)50813 Temporal Discretization: The Transient Term ðqC /C ÞtþDt=2 ðqC /C ÞtDt=2VC þ L /tC ¼ 0Dtð13:57ÞTo derive the full discretized equation, an interpolation profile expressing the facevalues at ðt Dt=2Þ and ðt þ Dt=2Þ in terms of the element values at ðtÞ, ðt DtÞ,etc., is needed.
The selection of this profile can heavily rely on the understandinggained from the discretization of the convection term. The choice will obviouslyaffect the accuracy and robustness of the method. In that regard, it is worth mentioning that the integration of the spatial operator is second order in time, but theaccuracy of the operator itself is determined by the options used during itsdiscretization.Independent of the profile used, the flux will be linearized based on old and newvalues asFluxT ¼ FluxC/C þ FluxC /C þ FluxVð13:58Þwhere again superscript ° refers to old values.
With the linearization completed, thecoefficients of the algebraic equation can then be assembled intoaCbCaC þ FluxCbC FluxC /C FluxVð13:59ÞIn what follows the discretization for a number of interpolation profiles is presented.13.3.1First Order Transient SchemesThe first order implicit and explicit Euler schemes will be constructed next byadopting an upwind [14, 17] and a downwind [4, 18] transient interpolation profile,respectively.13.3.2First Order Implicit Euler SchemeThe transient first order implicit Euler scheme is obtained by using a first-order“upwind” interpolation profile [14, 17].
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
