The development of High-Temperature Gas-Cooled Reactors (HTGRs) and Fluoride-Cooled High-Temperature Reactors (FHRs) that utilize pebble fuel has drastically increased the demand for improving the capabilities to simulate the packed beds found in these reactors. The complex flow fields found in a pebble bed make computational fluid dynamics (CFD) simulations time consuming and costly. Intermediate fidelity porous media models, however, are capable of approximating these flow fields in a much more computationally efficient manner. These models require the use of closures to capture the effects of complex flow phenomena without modeling them explicitly. This research employs data obtained from high-fidelity CFD simulations of a pebble bed to improve the drag closures used in porous media models in the near-wall region of the bed. Specifically, NekRS, a GPU-enabled spectral element CFD code, was used to simulate a bed of 1,568 pebbles at Reynolds numbers of 2,500, 5,000, and 10,000. The case was divided into five concentric subdomains to extract radial profiles of the average porosity, velocity, and wall shear in each subdomain. A model consistent with the high-fidelity model was created in Idaho National Laboratory's Pronghorn porous media code and the KTA correlation was chosen as the drag closure of comparison. Averages for the velocity, friction factor, and form factor were extracted in each region of the Pronghorn results to provide the comparative data. It was found that the KTA correlation significantly overestimates the velocity in the near-wall region. An investigation of the drag coefficients between the two codes revealed that the KTA correlation underestimated the form factor in the outermost region while overestimating it in the inner four regions. Manually increasing this form factor in the Pronghorn case predicted a velocity distribution in better agreement with the NekRS result.

An investigation of the pressure drop was also performed, in which the NekRS results were compared to the KTA correlation, the manually adjusted Pronghorn values, and the Wentz/Thodos correlation. It was found that the KTA correlation overpredicted the pressure drop at all three Reynolds numbers. The Wentz/Thodos correlation performed similarly to the KTA correlation, although it underpredicted the pressure drop in Re = 5,000 and Re = 10,000 cases. Although the manually increased drag coefficients saw better velocity agreement with the NekRS result, the increase in the form factor worsened the overprediction of the pressure drop.

This analysis in this work has revealed the underlying inaccuracy in the near-wall region of the KTA correlation and has set up the process for using high-fidelity simulation to predict more accurate drag coefficients a priori, rather than with a manual velocity-matching approach. This process will allow for the devlopment of an improved drag closure for use in porous media models.

§ INTRODUCTION

Packed beds are an integral part of emerging generation-IV nuclear reactors. Some of these reactors plan to utilize pebble fuel in which the fissile material is embedded into spheres of graphite. These graphite spheres are loaded into the reactor vessel, creating a randomly packed bed of spheres. The coolant in these designs, typically either an inert gas or a molten salt, flows through the packed bed, creating a complex flow field. Understanding the intricacies of this complex flow is critical to advancing the deployment of pebble bed reactors.

This work focuses in particular on an investigation of the wall-channel effect in packed beds. This effect is found in the near-wall region of a packed bed where the packing of the pebbles is disrupted by the presence of the wall. This disruption causes them to pack less efficiently, creating more space between pebbles and thus increasing the porosity near the wall. The increased porosity has a number of implications on the flow velocity, pressure drop, and heat transfer in this region and thus must be understood in order to accurately model Pebble Bed Reactors (PBRs).

Accurately modeling the wall-channel effect in a packed bed requires an advanced understanding of the pressure drop in each region of the system. Many correlations have been previously developed to predict the pressure drop. These correlations typically represent the pressure loss as the sum of viscous (Darcy) and inertial (Forccheimer) energy losses. One of the first major correlations to predict the pressure drop in packed beds in this manner was developed by Ergun [1]. Mehta and Hawley [2] later found that Ergun's correlation saw significant discrepancy from experimental results for $d_{bed}/d_{peb} < 50$ as a result of the larger role that wall effects play in low bed-to-pebble diameter ratio cases. They suggested improvements to Ergun's correlation by redefining the equivalent diameter to include the size of the bed to account for wall effects. This correlation was later improved by Riechelt [3] through the use of experimental studies to account for a larger range of Reynolds numbers, porosities, and bed-to-pebble diameter ratios. Hicks [4] also provided improvements to the Ergun equation by suggesting that the constants in front of the viscous and inertial terms are not constant, but rather depend on the Reynolds number. More recently, Choi [5] developed a model based on the Ergun equation for small bed-to-pebble diameter rations. This correlation included a wall correction factor for the inertial loss term. Wu et. al. [6] performed a study to investigate the effect of bed height on the pressure drop and developed a correlation to reflect this effect. The German Nuclear Safety Standards Comission (KTA) developed a correlation for specific application to High Temperature Gas-Cooled Reactors (HTGRs) [7]. This correlation is the main comparative correlation for this study for this reason. Wentz and Thodos have also developed a correlation for a wide range of Reynolds numbers and porosities [8]. This wide range of validity makes this correlation valuable for comparison in this study as well.

Other studies similar to this work have also been conducted to compare CFD data to available correlations. Yildiz et. al. [9] have performed DNS simulations to compare the results to multiple of the previously mentioned correlations. Atmakidis and Kenig [10] have also performed a study analyzing the wall effects in regular and irregular packed beds to compare to available correlations. Jun et. al. [11]have simulated bypass flow rates in the HTR-PM HTGR.

Previous work in this line of study focused on comparing the velocity distributions between high-fidelity LES simulation and intermediate-fidelity porous media simulation [12]. In this work, the KTA correlation was used as the drag coefficient in Idaho National Lab's Pronghorn porous media code. Two cases of 1,568 pebbles and 45,000 pebbles were created and simulated in NekRS at a Reynolds number of 20,000. A visualization of the 1,568 pebbles case with the velocity field can be found in Figure <ref>.

1,568 pebbles case from the previous work in this line of study with instantaneous velocity field shown on the left side.

The cases were separated into five concentric subdomains in postprocessing to extract porosity and velocity averages as a function of the radial position. The four outermost rings had a thickness of 0.25$d_{peb}$ each and the inner ring served contained the rest of the bed. The use of five rings was chosen to provide sufficient resolution of the near-wall region while also allowing for the combining of rings in the future if a coarser description of the wall region is desired. Two cases consistent with the NekRS cases were created in Pronghorn. It was found that the KTA correlation was causing the velocity in the near-wall region to be vastly overestimated and the distance needed to redistribute amongst the five subdomains was significantly longer than what was seen in the NekRS case. Manual changes to the drag coefficients in the near-wall region demonstrated significantly better agreement with the NekRS results, and a reduction in the non-flow direction drag coefficients improved the flow redistribution distance as well. This work demonstrated that capturing the wall-channeling effect in the porosity alone of a porous media model is not sufficient to reproduce the effect on the velocity field. It was determined that an improved correlation would be necessary to improve the ability to model the drag in the outer region. In this previous work, pebble-to-pebble contact points were handled by slightly shrinking the diameter of each pebble such that the pebbles do not touch in the resulting mesh. This approach, however, creates pressure gradients that are significantly lower than what is expected from a packed spherical bed. The major change to the mesh for this work is that the pebble contact is now included in the mesh, and a small chamfer is placed at the point of contact to prevent singular mesh points from appearing.

This work aims to continue the investigation of porous media codes' ability to model the wall-channel effect by examining the pressure drop in high-fidelity and intermediate-fidelity simulations. An investigation of the pressure drop provides additional insight into the ratio of friction and form losses in the different regions of a packed bed and will identify the areas of discrepancy that are seen with available correlations. Large-scale simulation with the Nek5000/NekRS spectral element code provide the benchmark data, and the KTA and Wentz/Thodos correlations are used are the pressure drop correlations of choice.

§ METHODS

§.§ NekRS

The high-fidelity code chosen for this work is Argonne National Laboratory's NekRS spectral element code. This code is a relatively new GPU-oriented version of the established open-source spectral element code Nek5000 [13]. NekRS is written in C++ and it is able to link to Nek5000 to leverage the existing extensive pre- and post-processing utilities. It was created as part of the DOE's Exascale Computing Project. NekRS is able to realize high throughput on advanced GPU nodes and demonstrates excellent scalability, and further details of NekRS performance for nuclear application have been provided in a recent publication [14].

Vadlidation of the Nek5000/NekRS code has previously been performed in a number of existing works. Specifically for pebble beds, work by Yildiz et. al [9] has demonstrated good agreement of first-order statistics between Nek5000 and experimental results. This work also demonstrated that the meshing method of applying chamfers at pebble contact areas does not significantly influence the porosity of the resulting bed. Additional verification and validation studies for Nek5000/NekRS can be found in works by Lai et. al [15] and Obabko et. al [16].

Calculations performed as part of this work are wall-resolved Large Eddy Simulations (LES). A high pass filter approach is used for to mimic the effect of dissipation by the subgrid scales [17].

The mesh for this study was created using a novel voronoi cell approach developed as part of the Cardinal multi-physics project [18]. First, pebble centers were generated by performing a Discrete Element Method (DEM) calculation with a target bed porosity of 40%. This method then generates an all-hexahedral mesh for the pebble void regions based on a tessellation of the voronoi cells that are based on the previously generated pebble centers. This method assumes that all pebbles are initially perfectly spherical, then a small chamfer is placed in locations where pebbles contact each other or the wall. The resulting mesh contains 1,568 pebbles at an average porosity of roughly 44%. A cross section of the mesh can be found in Figure <ref>. Further information on the meshing strategy can be found in [19].

The model was simulated using NekRS on the Summit supercomputer across 228 GPU's. The case was initally simulated at Reynolds numbers of 2,500, 5,000, and 10,000 based on the pebble diameter. Cases were run until the ring-average velocity results converged to a constant value.

§.§ Postprocessing of CFD data

A key aspect of this work has been the development of a post-processing framework to reduce the CFD data to a form usable in comparisons with porous media models.
Postprocessing was  performed on the time-averaged fields to extract porosity, velocity, pressure, and wall-shear data from five concentric regions. These quantities were extracted from an axial slice of height $2.5 d_{peb}$ in the center of the bed to avoid the influence of entrance or exit effects.

In order to determine the friction and form coefficients from this data, the following process was used. First, the wall shear stress was used to determine the friction factor using the following formula:

\begin{equation}
f = \frac{{8}\tau_{wall}}{\rho{v_s}^2}
\end{equation}

Then, the pressure gradient from frictional losses was determined with the Darcy-Weisbach equation.

\begin{equation}
\frac{\Delta{P}}{L}_{friction}= f\left(\frac{\rho}{2}\right)\left(\frac{{v_s}^2}{D_h}\right) \qquad \textnormal{where} \quad D_h = \frac{\epsilon{d_{peb}}}{1-\epsilon}
\end{equation}

The pressure gradient from form losses was then calculated by subtracting the frictional gradient from the total pressure gradient

\begin{equation}
\frac{\Delta{P}}{L}_{form}=\frac{\Delta{P}}{L}_{total}-\frac{\Delta{P}}{L}_{friction}
\end{equation}

Finally, the form coefficient was determined from the form pressure gradient.

\begin{equation}
C_{form} = \frac{\frac{\Delta{P}}{L}_{form}}{0.5\rho{v_s}^2}
\end{equation}

§.§ Pronghorn

Idaho National Laboratory's Pronghorn porous-media code was used as the intermediate-fidelity code for this work. Pronghorn is a finite element thermal-hydraulics code built on the Multiphysics Object Oriented Simulation Environment (MOOSE) framework. It is intended to provide steady-state and transient simulation results in short code execution times to assist in design scoping studies or to provide boundary conditions for system-level analysis. Pronghorn is capable of multiple solver models, but the compressible porous media Navier-Stokes model was chosen for this work. The strong form of this model is given in [20].

\begin{equation} \label{eqn1}
\epsilon\frac{\partial{\rho}}{\partial{t}}+\nabla\cdot(\epsilon\rho_{f}\vec{V})=0,
\end{equation}

\begin{equation}
\begin{aligned}
\epsilon\frac{\partial{(\rho_{f}\vec{V})}}{\partial{t}}+\nabla\cdot(\epsilon\rho_{f}\vec{V}\vec{V})+\epsilon\nabla{P} -\epsilon\rho_{f}\vec{g}+W\rho_{f}\vec{V} -\nabla\cdot(\mu\nabla\vec{V})=0,
\end{aligned}
\end{equation}

\begin{equation}
\begin{aligned}
\epsilon\frac{\partial{(\rho_{f}{E}_{f})}}{\partial{t}}+\nabla\cdot(\epsilon{H}_{f}\rho_{f}\vec{V})-\nabla\cdot(\kappa_{f}\nabla T_{f})-\epsilon\rho_{f}\vec{g}\cdot\vec{V} +\alpha(T_{f}-T_{s})-\dot{q_{f}} = 0,
\end{aligned}
\end{equation}

\begin{equation}
(1-\epsilon)\rho_{s}{C}_{p,s}\frac{\partial{T_{s}}}{\partial{t}}-\nabla\cdot(\kappa_{s}\nabla{T}_{s})+\alpha(T_{s}-T_{f})-\dot{q_{s}} = 0.
\end{equation}

A model consistent with the NekRS model was created in Pronghorn to provide comparative data. A nondimensionally identical RZ mesh was used to create these models. This mesh included the five subdomains where the porosity values extracted from the NekRS model were explicitly assigned. The Pronghorn mesh is pictured in Figure <ref> overlayed over a centerplane slice of the NekRS model. Helium was selected as the coolant fluid, and the inlet velocity was scaled to achieve each of the three average Reynolds numbers of the NekRS cases. Averaging was performed in the center of the model in the same region as the NekRS model. The average velocity, pressure, friction coefficient, and form coefficient were extracted for each of the five subdomains.

Centerplane slice of the NekRS mesh with the Pronghorn mesh overlayed on the right side. The averaging region where values were extracted from in both cases is indicated with the red dashed line.

§.§ KTA Drag Correlation

The KTA drag correlation developed by the German Nuclear Safety Standards Commission was the primary correlation that was investigated for this study. This correlation was developed specifically for application to HTGR's and is valid only over a narrow range of porosities from $0.36 < \epsilon < 0.42$. It should be noted that the bed used for this study had an average porosity slightly outside of this range.  The KTA correlation also does not account for any wall-channeling effects. The range of valid Reynolds numbers, however, is quite large at $10 < {Re}_{m} < 100,000$. The correlation can be found in Equation <ref>.

\begin{equation}
\frac{\Delta{P}}{L} = \left(\frac{320}{{Re}_{m}}+\frac{6}{{{Re}_{m}}^{0.1}}\right)\left(\frac{1-\epsilon}{\epsilon^3}\right)\left(\frac{\rho^2{v_{s}}^2}{{d}_p}\right)\left(\frac{1}{2\rho}\right)
\label{eq::KTA}
\end{equation}

§.§ Wentz and Thodos Drag Correlation

The second drag correlation used to compare in this work was the Wentz and Thodos correlation. This correlation was selected due to the large range of valid porosities compared to other correlation, with the range quoted as $0.354 < \epsilon < 0.882$, meaning that the simulation case is in the range of the correlation. This correlation as well does not explicitly account for any wall-channeling effects. The range of valid Reynolds number is smaller than the KTA correlation, but the cases used for this work are within the range of $2,540 < Re_{m} < 64,900$. It was developed using pebbles of 3.1242 cm diameter and a bed-to-pebble diameter ratio of 11.382. The correlation can be found in Equation <ref>.

\begin{equation}
\frac{\Delta{P}}{L} = \left(\frac{\mu{v_s}}{{d_{p}}^2}\right)\left(\frac{(1-\epsilon)^2}{\epsilon^3}\right)\left(\frac{0.396{Re}_{m}}{{{Re}_{m}}^{0.05}-1.20}\right)
\label{eq::WentzThodos}
\end{equation}

§ VERIFICATION AND VALIDATION ACTIVITIES

Validation of the NekRS results was performed through the use of comparative experimental data provided by Texas A&M University for a reduced experimental bed.

A section of the top of the test setup was recreated using the known locations of the centers of the 67 spheres in the reconstruction region. This case was simulated in NekRS at the three Reynolds numbers used in the experiment: 757, 1,289, and 1,820. Each case was run for sufficiently long in order to establish a converged field and then the field was averaged in time. The three average velocity fields can be seen in Figure <ref>.

Experimental data on a slice of the domain was gathered using a Particle Image Velocimetry (PIV) technique. Data was provided through several sampled lines for comparison to the NekRS results. The slice and corresponding sampling lines may be seen in Fig. <ref>. Comparisons of the velocity magnitude over these lines can be found in Figs. <ref> and <ref>. NekRS is generally capable of matching the shape of the velocity distribution, but some discrepancies in the magnitude are visible. Additionally, there are some slight differences in the positions of the pebbles seen in this comparison. This may be attributed to uncertainty in the PIV measurements, although further work is being performed to determine the source of this disagreement.

Time-averaged velocity fields at the centerplane of the TAMU cases. Re = 757 (left), Re = 1,289 (center), Re = 1,820 (right)

Turbulent Kinetic Energy fields at the centerplane of the TAMU cases. Re = 757 (left), Re = 1,289 (center), Re = 1,820 (right)

Velocity magnitude slice from the NekRS result of the Re = 1,820 case used to compare to the experimental results. The two line samples used for this comparison are shown in red and pink.

In-plane velocity magnitude comparisons between the experimental and NekRS results at line 1 for Re = 757 (left), Re = 1,289 (center), Re = 1,820 (right).

In-plane velocity magnitude comparisons between the experimental and NekRS results at line 2 for Re = 757 (left), Re = 1,289 (center), Re = 1,820 (right).

§ ANALYSIS AND RESULTS

The pebble bed was simulated with NekRS at the three specified Reynolds numbers and the averaged fields can be found in Figure <ref>. The wall-channeling effect is visible in these images, as the velocity is noticeably higher in the near-wall region. From these average flow fields, the porosity, velocity, pressure, and wall-shear data was extracted. Velocity, pressure, friction factors, and form factors were retrieved from the corresponding Pronghorn cases for comparison. Figure <ref> shows the radial velocity profiles for each of the three Reynolds numbers. Similarly to the previous work, the KTA correlation was found to drastically overestimate the velocity in the outermost ring. The velocity in ring three was also moderately overestimated and the velocity in ring 4 was consistently underestimated. Form coefficients were manually altered to match the NekRS velocity results in the outermost, producing Pronghorn results that saw much better agreement with the NekRS case. A better understanding of why these changes cause the velocity profile to better reproduce the NekRS result can be obtained by investigating the friction and form coefficients between the two codes. Details on how the friction and form losses have been extracted from NekRS are provided in Section <ref>.

Time-averaged velocity fields at the centerplane of the 1,568 pebbles cases. Re = 2,500 (left), Re = 5,000 (center), Re = 10,000 (right)

Radial velocity profiles for the NekRS, Pronghorn KTA, and Pronghorn with manual adjustments cases

The friction factor in each ring at each of the three Reynolds numbers is pictured in Figure <ref>. The KTA correlation appears to overestimate the friction factor in all but the outermost ring for the Re = 2,500 case while a general underprediction of the friction factor is found for Re = 5,000 and Re = 10,000. Additionally, it is evident that the shape of the friction factor curve is much flatter in NekRS than what is determined by the KTA correlation. These values are also presented normalized by the first ring in Figure <ref>. Although there is a significant disagreement between the NekRS and KTA results, this factor was not altered for this study due to the fact that the form losses (i.e., Forccheimer) are much more dominant in this case. Since only the form coefficients were altered in the manual Pronghorn case, there is only a slight difference in the friction factors compared to the KTA values. This change is a result of the velocity in each ring being different than when the KTA correlation is used for both friction and form.

Friction factor in each ring for the NekRS, Pronghorn KTA, and Pronghorn with manual adjustments cases for Re=2,500 (left), Re=5,000 (center) and Re=10,000 (right)

Friction factor in each ring normalized by the first ring for the NekRS, Pronghorn KTA, and Pronghorn with manual adjustments cases for Re=2,500 (left), Re=5,000 (center) and Re=10,000 (right)

Figure <ref> shows the form coefficients for each of the various cases. Much like the friction coefficient, the form coefficient was overestimated compared to the NekRS result in all but the outermost ring. This issue is worsened when looking at the manually adjusted coefficients, as the coefficients in ring five needed to be increased to match the NekRS velocity results. A better understanding of why the manual coefficients produce a velocity distribution in agreement with NekRS can be obtained by analysing the normalized form factors in Fig <ref>. These results clearly display that the form factor in the near-wall region is much lower compared to the first ring with the KTA correlation than with NekRS. The manually adjusted Pronghorn values were adjusted such that the velocity in ring five matched the NekRS result, and the normalized form factor plots reveal that making these adjustments caused the normalized form factors in this region to match the NekRS result very well. This result is entirely expected and serves to further validate the methodology for producing a new correlation that may be used in Pronghorn to determine the correct coefficients a priori.

Form factor in each ring for the NekRS, Pronghorn KTA, and Pronghorn with manual adjustments cases for Re=2,500 (left), Re=5,000 (center) and Re=10,000 (right)

Form factor in each ring normalized by the first ring for the NekRS, Pronghorn KTA, and Pronghorn with manual adjustments cases for Re=2,500 (left), Re=5,000 (center) and Re=10,000 (right)

The pressure drop for the NekRS, KTA, and Wentz/Thodos cases can be found in Table <ref>. The KTA correlation overpredicted the pressure drop at all three Reynolds numbers. Based on the previous analysis of the friction and form coefficients, this result is expected since the form coefficients were overpredicted in all but the outermost region. Likewise, the manually adjusted KTA case overpredicts the pressure drop even more since the form coefficients in ring  five needed to be increased to better match NekRS's velocity distribution. This result implies that although the form coefficient in the near-wall region needs to be increased with respect to the form coefficient in ring one, the coefficients in all rings need to be reduced if the pressure drop is to be matched between NekRS and the new correlation. The Wentz/Thodos correlation was also investigated due to its large range of porosity and Reynolds number validity. It performed, however, roughly similarly to the KTA correlation with around 7% error compared to the NekRS result. The KTA correllation was much more consistent in its error, varying from 3-7% error while the Wentz/Thodos correlation got siginificantly less accurate for higher Re, varying from around 1% for the Re = 2,500 case to greater than 10% in the Re =  10,000 case.

Source     Re = 2,500     Re = 5,000     Re = 10,000
NekRS     7.46     6.87     6.37

KTA     7.91     7.19     6.58

Wentz/Thodos     7.59     6.50     5.67

Comparison of pressure gradients predicted by NekRS, KTA correlation, and Wentz/Thodos correlation. Given pressure drops are nondimenionalized as $P_0/{d_{peb}}$ where $ P_0 = P/{\rho}{v}^2$

After the initial analysis was performed on the three Reynolds numbers, data for additional Reynolds numbers was gathered ranging from Re = 625 to Re = 10,000. A comparison of the nondimensional pressure gradients determined with NekRS and from the two correlations is found in Figure <ref>. The same trend experienced with the initial 3 Reynolds numbers is seen again, where the Wentz/Thodos correlation appears to have better agreement with the NekRS results for low Reynolds numbers, but deviates further at higher Reynolds numbers. The KTA correlation, meanwhile, is fairly consistent in its overprediction of the pressure gradient for Re $>$ 3,750.

Dimensionless pressure gradient versus Reynolds number from NekRS, KTA correlation, and Wentz/Thodos correlation.

§ CONCLUSIONS

In this work, a pebble bed mesh was created with an average porosity of 0.44 and simulated using the spectral-element code NekRS. Cases were run for sufficiently long to produce a converged average flow field, and the average fields were extracted. Postprocessing was performed to determine the axial pressure profile, radial velocity profile, radial porosity profile, and radial wall shear profile of each of the three cases. With the pressure profiles known, the pressure gradients in the center of each case were determined and used for comparison to two correlations. The KTA correlation was chosen as the first comparative correlation due to its applicability to HTGR's. However, the porosity of the pebble mesh was slightly outside of the range of valid porosities for this correlation. The Wentz and Thodos correlation was chosen as the second comparative correlation, as it is said to have a wider valid porosity range that should be valid even in the near-wall region. The three Reynolds numbers used for this study were in the valid range of both correlations.

An investigation of the pressure drop revealed that the KTA tended to overestimate the pressure drop with all three Reynolds numbers. The Wentz and Thodos correlation showed roughly equal agreement with the NekRS results at around 7% error. It overpredicted the pressure drop for a Reynolds number of 2,500, but underpredicted the pressure drop for the Reynolds numbers of 5,000 and 10,000. Separating the pressure drop into the friction and form losses revealed that the pressure drop is largely dominated by form effects, with the form losses accounting for roughly 87% of the total pressure drop in the inner rings and roughly 95% in the near-wall region.

The velocity results from the two codes displayed an overprediction of the velocity in the near-wall region when using the KTA correlation in Pronghorn. The form coefficient in ring five was then manually increased so that the velocity in the near-wall region would match the results from the NekRS case. A closer investigation into the friction and form coefficients revealed the reason for this improved agreement. The friction factor from the KTA correlation was overpredicted compared to the NekRS result. It also was more sensitive to the porosity and Reynolds number of the ring, with the friction factor profile being much less flat than the profile from NekRS. The friction factor was not altered for this study since the form effects were shown to be dominant, so no significant change was observed when using the manually adjusted correlation. The form factor, like the friction factor, also displayed overprediction when using the KTA correlation. This overprediction explains why the KTA correlation produced a larger pressure gradient than what was seen in NekRS. When the form factor results were normalized about the innermost ring, it was evident that the form coefficient in ring five was lower than what was seen from the NekRS result. The manually adjusted value for ring five improved the velocity distribution. This in turn also leads to a slight worsening of the overprediction for the pressure drop since the form coefficient in ring five was increased even further. The normalized form factor plots, however, revealed that the adjusted values saw an improved agreement with the normalized NekRS results. This result aligns with intuition, as it is the magnitude of the coefficients that determines the pressure drop and the ratio between the coefficients in the rings that determines the velocity distribution.  With the process for this analysis now understood, it will be possible to determine improved coefficients for use in Pronghorn a priori, rather than manually adjusting them until the velocity results are more closely matched.

Given that the KTA correlation is being used here slightly outside its range of validity, it is not surprising that it overpredicts the pressure drop. Coefficients would need to be decreased in magnitude across all regions to obtain better agreement with the NekRS result. Nonetheless, a crucial result of this work is that the form coefficient in the near-wall region needs to be increased with respect to the innermost bulk ring to prevent the overprediction of the velocity in this outer region. With the ability to determine form coefficients directly from the NekRS results, additional Reynolds number cases were simulated to generate more data to obtain a correlation for the form coefficient in the near-wall region. The data gathered from this study suggests that the the Forccheimer term of the KTA equation may altered in the near-wall region to the following:

\begin{equation}
C_{Form} = \frac{8.4}{Re_m^{0.1}}
\end{equation}

The resulting form factors and velocity profiles when using this correction the outermost ring are pictured in Figures <ref> and <ref>.

Forccheimer term in the near-wall region from NekRS, KTA, and the corrected KTA correlations.

Radial velocity profiles for the NekRS, Pronghorn KTA, and Pronghorn with corrected KTA correlation.

The coefficients in the other four regions will still need to be reduced if the full NekRS velocity profile and pressure drop is to be matched. Future work will investigate the influence of the pebble packing ratio on the resulting form coefficient in the near-wall region, as this may significantly alter the agreement seen with the KTA correlation. Additionally, future work will aim to investigate and improve correlations for heat transfer in the near-wall region.

[A]CFDComputational Fluid Dynamics
[A]PBRPebble Bed Reactor
[A]HTGRHigh-Temperature Gas-Cooled Reactor
[A]FHRFlouride-Cooled High-Temperature Reactor
[A]MOOSEMultiphysics Object-Oriented Simulation Environment
[A]DEMDiscrete Element Method
[A]LESLarge Eddy Simulation
[A]GPUGraphics Processing Unit
[A]DOEDepartment of Energy
[A]DNSDirect Numerical Simulation
[v]$D_h$Hydraulic Diameter
[v]$\rho$Fluid Density
[v]$Re$Reynolds Number
[v]$Re_m$Modified Reynolds Number ($\frac{Re}{1-\epsilon})$
[v]$\tau_w$Wall Shear Stress
[v]$v_i$Interstitial Velocity
[v]$v_s$Superficial Velocity
[v]$d_{peb}$Pebble Diameter
[v]$f$Darcy Friction Factor
[v]$C_{form}$Forccheimer Form Factor

[1]
S. Ergun.
Fluid flow through packed columns.
Journal of Chemical Engineering Progress, 48(2):89–94, 1952.

[2]
D. Mehta and M. Hawley.
Wall effect in packed columns.
Industrial & Engineering Chemistry Process Design and
Development, 8(2):280–282, 1969.

[3]
W. Reichelt.
Zur berechnung des druckverlustes einphasig durch- stromter kugel-
und zylinderschuttungen.
Rheological Acta, 17:676–692, 1978.

[4]
R.E. Hicks.
Pressure drop in packed beds of spheres.
Industrial & Engineering Chemistry Fundamentals,
9(3):500–502, 1970.

[5]
Y.S. Choi, S.J. Kim, and D. Kim.
A semi-empirical correlation for pressure drop in packed beds of
spherical particles.
Transport in Porous Media, 75(2):133–149, 2008.

[6]
J.S. Wu, B. Yu, and M. Yun.
A resistance model for flow through porous media.
Transport in Porous Media, 71(3):331–343, 2008.

[7]
K. Ausschusses and KTA.
Reactor core design of high-temperature gas-cooled reactors - part 3:
Loss pressure through friction pebble bed cores.

[8]
C.A. Wentz and G. Thodos.
Pressure drops in the flow of gases through packed and distended beds
of spherical particles.
AlChE, 9(1):81–84, 1963.

[9]
M. A. Yildiz, G. Botha, H. Yuan, E. Merzari, R. Kurwitz, and Y. A. Hassan.
Direct numerical simulation of the flow through a randomly packed
pebble bed.
Journal of Fluids Engineering, 142(4), 2020.

[10]
E.Y. Kenig and T. Atmakidis.
Cfd-based analysis of the wall effect on the pressure drop in packed
beds with moderate tube/particle diameter ratios in the laminar flow regime.
Chemical Engineering Journal, 155(1-2):404–410, 2009.

[11]
S. Jun, Z. Yanhua, and L. Fu.
Various bypass flow paths and bypass flow ratios in htr-pm.
Energy Procedia, 39:258–266, 2013.

[12]
D. Reger, E. Merzari, P. Balestra, S. Schunert, and Y.A. Hassan.
Comparison of pebble bed velocity profiles between high-fidelity and
intermediate-fidelity codes.
In Internaltional Conference On Nuclear Engineering, August

[13]
Paul Fischer, J Lottes, S Kerkemeier, O Marin, K Heisey, E Obabko, E Merzari,
and Y Peet.
Nek5000 user documentation.
Argonne National Laboratory, Lemont, IL, Report No.
ANL/MCS-TM-351, 2016.

[14]
Elia Merzari, Paul Fischer, Misun Min, Stefan Kerkemeier, Aleksandr Obabko,
Dillon Shaver, Haomin Yuan, Yiqi Yu, Javier Martinez, Landon Brockmeyer,
et al.
Toward exascale: overview of large eddy simulations and direct
numerical simulations of nuclear reactor flows with the spectral element
method in nek5000.
Nuclear Technology, 206(9):1308–1324, 2020.

[15]
J. Lai, E. Merzari, Y.A. Hassan, P. Fischer, and O. Marin.
Verification and validation of large eddy simulation with nek5000 for
cold leg mixing benchmark.
Nuclear Engineering and Design, 358, 2020.

[16]
A. Obabko, P. Fischer, O. Marin, E. Merzari, and D. Pointer.
Verification and validation of nek5000 for t-junction, matis,
siberia, and max experiments.
In Internation Conference on Nuclear Reactor Thermal
Hydrualics, September 2015.

[17]
Saleh Rezaeiravesh, Ricardo Vinuesa, and Philipp Schlatter.
On numerical uncertainties in scale-resolving simulations of
canonical wall turbulence.
Computers & Fluids, page 105024, 2021.

[18]
E. Merzari, H. Yuan, M. Min, D. Shaver, R. Rahaman, P. Shriwise, P. Romano,
A. Talamo, Y. Lan, D. Gaston, R. Martineau, P. Fischer, and Y. A. Hassan.
Cardinal: A lower length-scale multiphysics simulator for pebble-bed
Nuclear Technology, 2021.

[19]
Y. Lan, P. Fischer, E. Merzari, and M. Min.
All-hex meshing strategy for densely packed spheres.

[20]
A. Novak, S. Schunert, R. Carlsen, P. Balesta, D. Andrs, J. Kelly,
R. Slaybaugh, and R. Martineau.
Pronghorn Theory Manual, 2020.


# Recent experimental progress on probing Lorentz Violation in pure gravity
for $d\\!=\\!6$

Tao Jin  Jia-Rui Li  Yu-Jie Tan  and Cheng-Gang Shao  MOE Key Laboratory of
Fundamental Physical Quantities Measurements $\&$ Hubei Key
Laboratory of Gravitation and Quantum Physics,
PGMF and School of Physics, Huazhong University of Science and Technology,
Wuhan 430074, Peoples Republic of China

###### Abstract

Short-range experiment with striped design of the test and source masses
offers uniquely sensitive probes for Lorentz Violation. In our previous work,
we proposed to combine the horizontally and vertically striped experiments to
constrain the violating parameters. Here, we further point out adopting the
vertically striped structure with the experimental setup placed at two azimuth
angles, which has more experimental operability. For the experiment, we have
completed the assembly and measurement of the torsion pendulum.

##  1 Introduction

Lorentz symmetry, the idea that physical laws are unchanged under rotations
and boosts, is the cornerstone of both General Relativity and the Standard
Model. In SME frame, the pure gravity Lorentz Violation (LV) effects cause
corrections to the Newton potential, which are related to the violating
coefficients[1]. As the corrected potentials for $d\\!\geq{}\\!6$ are
inversely proportional to power of distance, the short-range gravity
experiments are sensitive to test these LV effects. By analyzing the
experimental data of testing the gravitational inverse square law, our group
collaborating with Indiana University has given the best LV constraints for
$d=6$ [2]. Moreover, to obtain a better constraint, we designed a special
experimental scheme for LV test, in which both the test and source masses are
in striped structure to enhance the LV signal[3, 4], and the experiment is
under going. In this work, we mainly introduce the experimental progress and
the related theoretical analysis.

##  2 Experimental scheme of LV test

The experimental principle has been described in the previous work[3, 4], and
the experimental apparatus is similar to the inverse-square experiment
HUST-2012, shown in Fig.1. The LV force between test and source masses varies
with earth’s rotation by the sidereal time. Meanwhile, we introduce another
modulation by changing the distance between the test and source masses at
frequency $\omega_{s}$ for reducing the torque noise.

Figure 1: Schematic diagram of the LV experimental setup (not to scale). The
pendulum twist is measured by an autocollimator, and controlled by two
differential capacitive actuators. The rotating copper cylinder is set to
calibrate the sensitivity of the pendulum.

The measured LV torque in the experiment can be expressed as:

$\tau_{\textup{measured}}^{z}(T)=\tau_{\textup{LV}}(T)\cos(\omega_{s}T+\varphi),$
(1)

with

$\tau_{\textup{LV}}(T)=C_{0}+\sum_{m=1}^{4}{C_{m}\cos(m\omega_{\oplus}T)+S_{m}\sin(m\omega_{\oplus}T)}.$
(2)

Here, $\omega_{\oplus}\approx 2\pi/{}$(23 h 56 min) is the Earth’s sidereal
frequency, and $\varphi$ is initial phase. The nine Fourier amplitudes
$C_{0},C_{m},S_{m}$ can be expressed in terms of the LV coefficients
$k_{jm}^{N(6)}$ and the transfer coefficients $\Gamma_{j}$, shown in Table.2.
$\Gamma_{j}$ is related to the geometries of the test and source masses, the
colatitude $\chi$, and the angle $\theta$ between the $x$ axis and local
south. To pursuit the good constraints, any experimental design should make
$\Gamma_{j}$ as large as possible. According to the typical design parameters,
we numerically calculated the transfer functions for the test and source
masses in horizontal stripes, as well as vertical case. The result shows
$\theta=3\pi/5$ and $\pi$/6 in the vertically striped case are the good choice
for the experiment (the corresponding values of the transfer coefficients are
listed in Table.2, which is more sensitive than that of HUST-2012 almost one
order of magnitude). Thus, we can use one experimental setup to construct two
experiments by setting the angle at these two values, respectively. Combining
the experimental data, the 14 LV parameters can be constrained independently.

Transfer coefficients $\Gamma_{j}$ for the vertically striped experiment.
Quantity Expression Transfer HUST-2012 New experimental design coefficients
$\theta=3\pi/5$ $\theta=\pi/6$ $C_{0}$ $\Gamma_{1}k_{2,0}+\Gamma_{2}k_{4,0}$
$\Gamma_{1}$ -0.9 -25.3 -10.7 $\Gamma_{2}$ 0.0 19.0 -8.2 $C_{2}$
$\Gamma_{3}\textup{Re}k_{2,2}+\Gamma_{4}\textup{Im}k_{2,2}$ $\Gamma_{3}$ -0.4
17.8 12.9 $+\Gamma_{5}\textup{Re}k_{4,2}+\Gamma_{6}\textup{Im}k_{4,2}$
$\Gamma_{4}$ 0.3 -14.5 36.4 $S_{2}$
$\Gamma_{4}\textup{Re}k_{2,2}-\Gamma_{3}\textup{Im}k_{2,2}$ $\Gamma_{5}$ -0.9
-14.1 -14.7 $+\Gamma_{6}\textup{Re}k_{4,2}-\Gamma_{5}\textup{Im}k_{4,2}$
$\Gamma_{6}$ -0.4 4.1 4.9 $C_{4}$
$\Gamma_{7}\textup{Re}k_{4,4}+\Gamma_{8}\textup{Im}k_{4,4}$ $\Gamma_{7}$ -3.7
1.5 22.9 $S_{4}$ $\Gamma_{8}\textup{Re}k_{4,4}-\Gamma_{7}\textup{Im}k_{4,4}$
$\Gamma_{8}$ -0.4 14.3 33.1 $C_{1}$
$\Gamma_{9}\textup{Re}k_{2,2}+\Gamma_{10}\textup{Im}k_{2,2}$ $\Gamma_{9}$ -1.6
20.5 13.3 $+\Gamma_{11}\textup{Re}k_{4,2}+\Gamma_{12}\textup{Im}k_{4,2}$
$\Gamma_{10}$ -0.5 -4.0 18.7 $S_{1}$
$\Gamma_{10}\textup{Re}k_{2,2}-\Gamma_{9}\textup{Im}k_{2,2}$ $\Gamma_{11}$ 1.0
18.9 8.5 $+\Gamma_{12}\textup{Re}k_{4,2}-\Gamma_{11}\textup{Im}k_{4,2}$
$\Gamma_{12}$ -0.4 -19.4 2.9 $C_{3}$
$\Gamma_{13}\textup{Re}k_{4,4}+\Gamma_{14}\textup{Im}k_{4,4}$ $\Gamma_{13}$
-0.1 -6.4 16.8 $S_{3}$
$\Gamma_{14}\textup{Re}k_{4,4}-\Gamma_{13}\textup{Im}k_{4,4}$ $\Gamma_{14}$
-0.4 16.7 -22.6

##  3 Error analysis and experimental progress

The LV coefficients $k_{jm}^{N(6)}$ can be constrained by a fit of the Fourier
amplitudes $C_{0},C_{m},S_{m}$ data. Thus, the constrained levels of LV
coefficients are determined by the uncertainties of the Fourier amplitudes.
The value of the dc component $C_{0}$ extracted from recorded data
$\tau_{\textup{measured}}^{z}(T)$ is restricted by the systematic errors
arising from uncertainties in the dimensions and locations of the test and
source masses. The sidereal-harmonic amplitudes $C_{m},S_{m}$ are extracted
from $\tau_{\textup{LV}}(T)$ demodulated from the recorded data
$\tau_{\textup{measured}}^{z}(T)$, and their uncertainties are dominated by
the statistical uncertainty which is at the same level $0.4\times 10^{-16}$ Nm
for each harmonic.

Figure 2: Assembling and measurement of pendulum.

At present, the experimental components, such as striped tungsten slices and
glass blocks, have been machined; the geometry parameters of the tungsten
slices and some core parts have been measured in a double-blind way, and the
result shows the machining accuracy meets the experimental requirements; the
torsion pendulum has been assembled and measured, shown in Fig.2. Next, we
will suspend the pendulum to test the background, assemble the source mass and
put it into the vacuum chamber, and finally record the experimental data after
everything is normal.

##  4 Conclusion

In this work, we mainly review and slightly improve the experimental scheme we
proposed to test LV in the previous work, which can improve the current
constraints on the LV coefficients by more than one order of magnitude.
Currently, the experiment is proceeding as planned, and we hope the
experimental results can be obtained within two years.

## Acknowledgments

This work was supported by the National Natural Science Foundation of China
(Grant Nos. 12175076 and 11925503).

## References

  * [1] V.A. Kostelecký and M. Mewes, Phys. Lett. B 766, 137(2017)
  * [2] C.-G. Shao, Y.-J. Tan, W.-H. Tan, S.-Q. Yang, J. Luo, M.E. Tobar, Q.G. Bailey, J.C. Long, E. Weisman, R. Xu, V.A. Kostelecký, Phys. Rev. Lett. 117, 071102(2016)
  * [3] C.-G. Shao, Y.-F Chen, Y.-J Tan, J. Luo, S.-Q Yang and M.E Tobar, Phys. Rev. D 94, 104061(2016).
  * [4] Y.-F. Chen, Y.-J. Tan and C.-G. Shao, Symmetry 9, 217(2017).


# A note on holographic dark energy with varying $c^{2}$ term

A. Sheykhi1,2,111<EMAIL_ADDRESS>S. Ghaffari2 and N. Roshanshah1
1 Physics Department and Biruni Observatory, College of Sciences, Shiraz
University, Shiraz 71454, Iran
2Research Institute for Astronomy and Astrophysics of Maragha (RIAAM), P.O.
Box 55134-441, Maragha, Iran

###### Abstract

We reconsider the holographic dark energy (HDE) model with a slowly time
varying $c^{2}(z)$ parameter in the energy density, namely
$\rho_{D}=3M_{p}^{2}c^{2}(z)/L^{2}$, where $L$ is the IR cutoff and $z$ is the
redshift parameter. As the system’s IR cutoff we choose the Hubble radius and
the Granda-Oliveros (GO) cutoffs. The latter inspired by the Ricci scalar
curvature. We derive the evolution of the cosmological parameters such as the
equation of state and the deceleration parameters as the explicit functions of
the redshift parameter $z$. Then, we plot the evolutions of these cosmological
parameters in terms of the redshift parameter during the history of the
universe. Interestingly enough, we observe that by choosing $L=H^{-1}$ as the
IR cutoff for the HDE with time varying $c^{2}(z)$ term, the present
acceleration of the universe expansion can be achieved, even in the absence of
interaction between dark energy and dark matter. This is in contrast to the
usual HDE model with constant $c^{2}$ term, which leads to a wrong equation of
state, namely that for dust $w_{D}=0$, when the IR cutoff is chosen the Hubble
radius.

##  I Introduction

According to nowadays observations, present acceleration of the universe
expansion has been well established Riess . Within the framework of general
relativity, the responsible component of energy for this accelerated expansion
is known as dark energy (DE) with negative pressure. However, the nature of DE
is still unknown, and some candidates have been proposed to explain it. The
earliest and simplest candidate is the cosmological constant with the time
independent equation of state $\omega_{\Lambda}=-1$ which has some problems
like fine-tuning and coincidence problems. Therefore, other theories have been
suggested for the dynamical DE scenario to describe the accelerating universe.

An interesting attempt for probing the nature of DE within the framework of
quantum gravity, is the so-called HDE proposal. This model which has arisen a
lot of enthusiasm recently Coh ; Li ; Huang ; Hsu ; HDE ; Setare1 ; wang0 ;
SheyHDE , is motivated from the holographic hypothesis Suss1  and has been
tested and constrained by various astronomical observations Xin . In
holographic principle a short distance cutoff could be related to a long
distance cutoff (infrared cutoff) due to the limit set by formation of a black
hole. Based on the holographic principle, it was shown by Cohen et al. Coh
that the quantum zero-point energy of a system with size $L$ should not exceed
the mass of a black hole with the same size, i.e.,

$L^{3}\rho_{D}\leq LM_{p}^{2},$ (1)

where $M^{2}_{p}=8\pi G$ is the reduced Planck mass and $L$ is the IR cutoff.
The largest $L$ allowed is the one saturating this inequality so that we get
the

$\rho_{D}=3c^{2}M_{p}^{2}/L^{2},$ (2)

where $3c^{2}$ is a dimensionless constant. There are many models of HDE,
depending on the IR cutoff, that have been studied in the literatures Horava ;
Fischler ; Nojiri ; Gao . The simple choice for IR cutoff is the Hubble
radius, i.e., $L=H^{-1}$ which leads to a wrong equation of state (EoS) and
the accelerated expansion of the universe cannot be achieved. However, as soon
as an interaction between HDE and dark matter is taken into account, the
identification of IR cutoff with Hubble radius $H^{-1}$, in flat universe, can
simultaneously drive accelerated expansion and solve the coincidence problem.
Pavon . Then, Li Li  showed that taking the particle horizon radius as IR
cutoff it is impossible to obtain an accelerated expansion. He also
demonstrated that the identification of $L$ with the radius of the future
event horizon gives the desired result, namely a sufficiently negative
equation of state to obtain an accelerated universe.

It is worth noting that, for the sake of simplicity, very often the $c^{2}$
parameter in the HDE model is assumed constant. However, there are no strong
evidences to demonstrate that $c^{2}$ should be a constant and one should bear
in mind that it is more general to consider it a slowly varying function of
time. It has been shown that the parameter $c^{2}$ can play an essential role
in characterizing the model. For example, it was argued that the HDE model in
the far future can be like a phantom or quintessence DE model depending
whether the parameter $c^{2}$ is larger or smaller than $1$, respectively
Radicella . By slowly vary function with time, we mean that
$(\dot{c}^{2})/(c^{2})$ is upper bounded by the Hubble expansion rate, i.e.,
Radicella

$\frac{({c}^{2})^{\dot{}}}{c^{2}}\leq H,$ (3)

where dot indicates derivative with respect to the cosmic time. In this case
the time scale of the evolution of $c^{2}$ is shorter than $H^{-1}$ and one
can be satisfied to consider the time dependency of $c^{2}$ Radicella .
Considering the future event horizon as IR cutoff, the HDE model with time
varying parameter $c^{2}$, has been studied in Malekjani . It was argued that
depending on the parameter $c^{2}$, the phantom regime can be achieved earlier
or later compared to the usual HDE with constant $c^{2}$ term Malekjani . In
this paper, we reconsider the HDE model with the slowly varying parameter
$c^{2}(z)$ by taking into account the Hubble horizon $L=H^{-1}$ and GO cutoff,
$L=(\alpha H^{2}+\beta\dot{H})^{{-1}/{2}}$, as the system’s IR cutoffs. We
shall study four parameterizations of $c(z)$ as follows

$\displaystyle\text{GHDE1}:~{}~{}~{}~{}c(z)=c_{0}+c_{1}\frac{z}{(1+z)},$ (4)
$\displaystyle\text{GHDE2}:~{}~{}~{}~{}c(z)=c_{0}+c_{1}\frac{z}{(1+z)^{2}},$
(5)
$\displaystyle\text{GHDE3}:~{}~{}~{}~{}c(z)=\frac{c_{0}}{1+c_{1}\ln(1+z)},$
(6)
$\displaystyle\text{GHDE4}:~{}~{}~{}~{}c(z)=c_{0}+c_{1}\left(\frac{\ln(2+z)}{1+z}-\ln
2\right),$ (7)

where GHDE stands for the Generalized HDE model. The above choices for $c(z)$
are, respectively, inspired by the parameterizations known as Chevallier-
Polarski-Linder parametrization (CPL) Chevallier , Jassal-Bagla-Padmanabhan
(JBP) parametrization Jassal , Wetterich parametrization Wetterich , and Ma-
Zhang parametrization Ma . Setting $c_{1}=0$ in all these four
parameterizations, the original HDE with constant $c$ parameter is recovered.

This paper is organized as follows. In section II, we drive the basic
equations for the HDE with time varying $c^{2}$ parameter. In this section, we
also consider the Hubble radius as IR cutoff and derive the evolution of EoS
and deceleration parameters by choosing $c(z)$. In section III, we repeat the
study for the GO cutoff and investigate the evolution of the cosmological
parameters. The last section is devoted to conclusions and discussions.

##  II GHDE in flat FRW cosmology with Hubble radius as IR cutoff

In the context of flat Friedmann-Robertson-Walker (FRW) cosmology, the
Friedmann equation can be written

$H^{2}=\frac{1}{3M_{p}^{2}}(\rho_{m}+\rho_{D}),$ (8)

where $\rho_{m}$ and $\rho_{D}$ are the energy densities of pressureless dark
matter and DE, respectively. By using the dimensionless energy densities

$\Omega_{m}=\frac{\rho_{m}}{3m_{p}^{2}H^{2}},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Omega_{D}=\frac{\rho_{D}}{3m_{p}^{2}H^{2}},$
(9)

the Friedmann equation (8) can be written as

$\Omega_{m}+\Omega_{D}=1.$ (10)

We shall assume there is no interaction between dark matter and GHDE.
Therefore, both components obey independent conservation equation. The
conservation equations for pressureless dark matter and DE, are given by

$\displaystyle\dot{\rho}_{m}+3H\rho_{m}=0,$ (11)
$\displaystyle\dot{\rho}_{D}+3H(1+\omega_{D})\rho_{D}=0,$ (12)

where $w_{D}=p_{D}/\rho_{D}$ is the EoS parameter of GHDE. In this section we
consider the Hubble radius as IR cutoff, $(L=H^{-1})$, thus the energy density
of GHDE model from (2) can be written as

$\rho_{D}=3M_{p}^{2}c^{2}(z)H^{2}.$ (13)

Using definition (9), the dimensionless energy density for the GHDE becomes

$\Omega_{D}(z)=c^{2}(z).$ (14)

Taking the time derivative of Eq. (13), we find

$\dot{\rho}_{D}=2\rho_{D}\left(\frac{\dot{c}(z)}{c(z)}+\frac{\dot{H}}{H}\right).$
(15)

Besides, if we take the time derivative of Friedmann equation (8), after using
Eqs. (10), (11) and (15), we find

$\frac{\dot{H}}{H^{2}}=-\frac{3}{2}+\frac{{c}^{\prime}(z)c(z)}{1-c^{2}(z)},$
(16)

where ${c}^{\prime}=\dot{c}/H$ and the prime represents derivative with
respect to $x=\ln a$. Combining Eqs. (15) and (16) with Eq. (12), one can
obtain the EoS parameter of GHDE as

$\omega_{D}=-\frac{2}{3}\frac{{c}^{\prime}(z)}{c(z)[1-c^{2}(z)]},$ (17)

Clearly, for constant $c$ parameter we have, ${c}^{\prime}(z)=0$ which leads
to a wrong equation of state, namely that for dust with $\omega_{D}=0$Hsu .
This implies that the HDE model with $L=H^{-1}$ as IR cutoff cannot describe
an accelerating universe Hsu . In contrast, taking the time varying $c^{2}$
term in the energy density of the HDE, it is quite possible to reproduce the
acceleration of the cosmic expansion in HDE model with the Hubble radius as IR
cutoff.

Another important cosmological parameter for studying the evolution of the
universe is the deceleration parameter which is given by

$q=1-\frac{\dot{H}}{H^{2}}.$ (18)

Substituting Eq. (16) into (18) yields

$q=\frac{1}{2}-\frac{{c}^{\prime}(z)c(z)}{1-c^{2}(z)}.$ (19)

Again for $c^{\prime}(z)=0$ the declaration parameter reduces to $q={1}/{2}>0$
which implies a decelerated universe.

We see from Eqs. (17) and (19) that the evolution of these cosmological
parameters depend on the functional form of $c(z)$. In what follow, we
consider four types of parametrization for $c(z)$ as given in Eqs. (4)-(7).

###  II.1 GHDE1: The CPL type

We start with the CPL type for $c(z)$, namely

$c(z)=c_{0}+c_{1}\frac{z}{(1+z)}.$ (20)

When $z\rightarrow\infty$ (in the early universe), we see that $c\rightarrow
c_{0}+c_{1}$ and as $z\rightarrow 0$ (at the present time), $c\rightarrow
c_{0}$. Thus holographic parameter $c$ varies slowly from $c_{0}+c_{1}$ to
$c_{0}$ from past to the present.

Using the fact that $a/a_{0}=(1+z)^{-1}$, where $z$ is the redshift parameter
and prime denotes the derivative with respect to $x=\ln a$, we arrive at
${c}^{\prime}=-(1+z){dc}/{dz}$. Taking derivatives of Eq. (20), it follows
that

${c}^{\prime}=-(1+z)\frac{dc}{dz}=-\frac{c_{1}}{1+z}.$ (21)

Substituting Eqs. (20) and (21) into (17), the EoS parameter is obtained as

$\omega_{D}(z)=\frac{2}{3}\frac{c_{1}(1+z)^{2}}{[c_{0}(1+z)+c_{1}z](1+z)^{2}-[c_{0}(1+z)+c_{1}z]^{3}}.$
(22)

Combining Eqs. (20) and (21) with Eq. (19), we can obtain the deceleration
parameter as

$q(z)=\frac{1}{2}+\frac{c_{0}c_{1}(1+z)+c_{1}^{2}z}{(1+z)^{2}-[c_{0}(1+z)+c_{1}z]^{2}}.$
(23)

The behavior of $\omega_{D}(z)$ and $q(z)$ are plotted for different values of
model parameters $c_{0}$ and $c_{1}$ in Fig.1. From these figures we see that
our Universe has a transition from deceleration to the acceleration phase
around $z\approx 0.6$ which is consistent with observations Daly ; Kom1 ; Kom2
; Planck1 ; Planck2 .

Figure 1: The evolution of EoS parameter $\omega_{D}$ and the deceleration
parameter $q$ versus redshift $z$ for GHDE1 model with $L=H^{-1}.$

###  II.2 GHDE2: The JBP model

The second parametrization of $c(z)$ is the JBP parametrization which is
written as

$c(z)=c_{0}+c_{1}\frac{z}{(1+z)^{2}}.$ (24)

From Eq. (24) we see that at the late time where $z\rightarrow 0$, we have
$c(z)\rightarrow c_{0}$, and as $z\rightarrow 1$, we have $c(z)\rightarrow
c_{0}+{c_{1}}/{4}$. Besides, in the early universe where $z\rightarrow\infty$,
we have $c(z)\rightarrow c_{0}$.

Taking derivative of Eq. (24) we find

${c}^{\prime}(z)=-c_{1}\frac{1-z}{(1+z)^{2}}.$ (25)

Substituting Eqs. (24) and (25) into Eq. (17), we obtain the EoS parameter for
GHDE2 model as

$\omega_{D}=\frac{2}{3}\frac{c_{1}(1-z)(1+z)^{4}}{[c_{0}(1+z)^{2}+c_{1}z][(1+z)^{4}-(c_{0}(1+z)^{2}+c_{1}z)^{2}]}.$
(26)

Combining Eqs. (24) and (25) with Eq. (19) yields

$q=\frac{1}{2}+\frac{c_{1}(1-z)[c_{0}(1+z)^{2}+c_{1}z]}{(1+z)^{4}-[c_{0}(1+z)^{2}+c_{1}z]^{2}}.$
(27)

The behavior of $\omega_{D}(z)$ and $q(z)$ are plotted for different values of
model parameters $c_{0}$ and $c_{1}$ in Fig. 2. Again, the universe has a
transition from deceleration to the acceleration phase around $z\approx 0.6$
and at the late time where $z\rightarrow 0$, the EoS parameter can cross the
phantom line $w_{D}=-1$.

Figure 2: The evolution of EoS parameter $\omega_{D}$ and the deceleration
parameter $q$ versus redshift $z$ for GHDE2 model with $L=H^{-1}.$

###  II.3 GHDE3: The Wetterich type

The third parametrization is Wetterich-type which assumes the following form
of $c(z)$:

$c(z)=\frac{c_{0}}{1+c_{1}\ln(1+z)}.$ (28)

In this model, at the late time where $z\rightarrow 0$, we have
$c(z)\rightarrow c_{0}$, while at the early universe where
$z\rightarrow\infty$, we get $c(z)\rightarrow 0$ and thus $\rho_{D}\rightarrow
0$. This implies that at the early universe the HDE did not have significant
contribution in the evolution of the universe. It follows directly that,

${c}^{\prime}(z)=-(1+z)\frac{dc(z)}{dz}=\frac{c_{0}c_{1}}{[1+c_{1}\ln(1+z)]^{2}}.$
(29)

Inserting Eqs. (28) and (29) into (17), one gets

$\omega_{D}(z)=-\frac{2}{3}\frac{c_{1}[1+c_{1}\ln(1+z)]}{[1+c_{1}\ln(1+z)]^{2}-c_{0}^{2}}.$
(30)

Combining Eqs. (19), (28) and (29), we find

$q(z)=\frac{1}{2}+\frac{c_{0}^{2}c_{1}}{[1+c_{1}\ln(1+z)]\Big{(}c_{0}^{2}-[1+c_{1}\ln(1+z)]^{2}\Big{)}}.$
(31)

We have plotted the evolutions of $\omega_{D}(z)$ and $q(z)$ in terms of the
redshift parameter $z$ in Fig. 3. From these figures it is obvious that the
present acceleration can be addressed in this model and the transition from
deceleration to the acceleration phase occurs for $0.3<z<0.7$.

Figure 3: The evolution of $\omega_{D}$ and $q$ versus redshift $z$ for GHDE3
with $L=H^{-1}.$

###  II.4 GHDE4: The Ma-Zhang type

The last choice for the parametrization of $c(z)$ was proposed in Ma  and can
be written as

$c(z)=c_{0}+c_{1}\left(\frac{\ln(2+z)}{1+z}-\ln 2\right).$ (32)

At the present time where $z\rightarrow 0$ we have $c(z)\rightarrow c_{0}$,
and at the early time where $z\rightarrow\infty$, one gets $c(z)\rightarrow
c_{0}-c_{1}\ln 2$. It is worth noting that for the previous choice of $c(z)$,
we could not investigate the future behavior of $c(z)$, because it diverges at
the future time where $z\rightarrow-1$. However, in case of the Ma-Zhang
parametrization we have $c(z)\rightarrow c_{0}+c_{1}(1-\ln 2)$ as
$z\rightarrow-1$. From Eq. (32) it is easy to show that

${c}^{\prime}(z)=-(1+z)\frac{dc(z)}{dz}=c_{1}\frac{(2+z)\ln(2+z)-(1+z)}{(1+z)(2+z)}.$
(33)

Combining Eqs. (32) and (33) with (17) and (19), we arrive at

$\omega_{D}(z)=\frac{2}{3}\frac{c_{1}(1+z)^{2}[(1+z)-(2+z)\ln(2+z)]}{(2+z)\Big{(}c_{0}(1+z)+c_{1}\ln(2+z)-c_{1}(1+z)\ln
2\Big{)}\Big{[}(1+z)^{2}-[c_{0}(1+z)+c_{1}\ln(2+z)-c_{1}(1+z)\ln
2]^{2}\Big{]}},$ (34)
$q(z)=\frac{1}{2}+\frac{c_{1}[c_{0}(1+z)+c_{1}\ln(2+z)-c_{1}(1+z)\ln
2][1+z-(2+z)\ln(2+z)]}{(2+z)\Big{(}(1+z)^{2}-[c_{0}(1+z)+c_{1}\ln(2+z)-c_{1}(1+z)\ln
2]^{2}\Big{)}}.$ (35)

In order to have an insight on the behaviour of these functions, we plot them
in terms of $z$ in Fig. 4. From these figures, we see that the behaviour is
similar to the previous parameterizations.

Figure 4: The evolution of EoS parameter $\omega_{D}$ and the deceleration
parameter $q$ versus redshift $z$ for GHDE4 with $L=H^{-1}.$

##  III GHDE in flat universe with GO cutoff

In this section we consider the GO cutoff as system’s IR cutoff, namely
$L=(\alpha H^{2}+\beta\dot{H})^{{-1}/{2}}$ which first proposed in Granda .
With this IR cutoff, the energy density (2) is written

$\rho_{D}=3M_{p}^{2}c^{2}(z)\left(\alpha H^{2}+\beta\dot{H}\right),$ (36)

where $\alpha$ and $\beta$ are constants that should be constrained by
observational data. Using the definition of density parameter (10) one can
obtain

$\Omega_{D}=c^{2}(z)\left(\alpha+\beta\frac{\dot{H}}{H^{2}}\right).$ (37)

Taking the time derivative of Eq. (37), we get

$\dot{\Omega}_{D}=2\Omega_{D}\left(\frac{\dot{c}(z)}{c(z)}-\frac{\dot{H}}{H}\right)+\frac{c^{2}(z)}{H^{2}}\left(2\alpha\dot{H}H+\beta\ddot{H}\right).$
(38)

Now, if we take the time derivative of both sides of Friedmann equation (8)
and after using Eq. (36), we arrive at

$c^{2}H^{-3}(2\alpha\dot{H}H+\beta\ddot{H})=2\frac{\dot{H}}{H^{2}}+3(1-\Omega_{D})-2\Omega_{D}\frac{\dot{c}(z)}{Hc(z)}.$
(39)

Combining Eqs. (38) and (39) the equation of motion for the dimensionless GHDE
density can be written as

$\dot{\Omega}_{D}=(1-\Omega_{D})\left(\dfrac{2\dot{H}}{H}+3H\right).$ (40)

By help of Eq. (37) and using the fact that
$\dot{\Omega}_{D}=H{\Omega}^{\prime}_{D}$, the evolution of dimensionless GHDE
density may be rewritten as

${\Omega}^{\prime}_{D}(z)=(1-\Omega_{D})\left(\frac{2\Omega_{D}}{\beta
c^{2}(z)}-\frac{2\alpha}{\beta}+3\right),$ (41)

where the prime denotes derivative with respect to $x=\ln{a}$. Taking the time
derivative of both sides of Eq. (8) and using Eqs. (12) and (37) we can obtain
the EoS parameter of GHDE model as follows

$\omega_{D}(z)=\frac{2\alpha}{3\beta\Omega_{D}}-\frac{2}{3\beta
c^{2}(z)}-\frac{1}{\Omega_{D}}.$ (42)

Substituting Eq. (37) into (18), we get

$q(z)=-1-\frac{\Omega_{D}}{\beta c^{2}(z)}+\frac{\alpha}{\beta}.$ (43)

Following the previous section, we shall consider four types of
parametrization of $c(z)$ listed in Eqs. (4)-(7), respectively.

Figure 5: The evolution of the dimensionless density parameter $\Omega_{D}$
versus redshift $z$ for GHDE1 model with GO cutoff.

###  III.1 GHDE1: The CPL type

Substituting Eq. (4) in (41), the equation of the evolutionary of
dimensionless GHDE1 density is obtained as

${\Omega}^{\prime}_{D}=(1-\Omega_{D})\left(\frac{2\Omega_{D}(1+z)^{2}}{\beta[c_{0}(1+z)+c_{1}z]^{2}}-\frac{2\alpha}{\beta}+3\right)$
(44)

The evolution of the dimensionless GHDE density parameter $\Omega_{D}$ as a
function of $1+z=a^{-1}$ is shown in Fig. 5. From this figure we see that at
the early universe where $z\rightarrow\infty$ we have $\Omega_{D}\rightarrow
0$, while at the late time where $z\rightarrow-1$, the DE dominated, namely
$\Omega_{D}\rightarrow 1$.
Using Eqs. (4), (42) and (43) we can obtain the EoS and the deceleration
parameters as

$\displaystyle\omega_{D}(z)$ $\displaystyle=$
$\displaystyle\frac{2\alpha}{3\beta\Omega_{D}}-\frac{2(1+z)^{2}}{3\beta[c_{0}(1+z)+c_{1}z]^{2}}-\frac{1}{\Omega_{D}},$
(45) $\displaystyle q(z)$ $\displaystyle=$
$\displaystyle-1-\frac{\Omega_{D}(1+z)^{2}}{\beta(c_{0}(1+z)+c_{1}z)^{2}}+\frac{\alpha}{\beta}.$
(46)

The evolution of EoS parameter $\omega_{D}(z)$ and deceleration parameter
$q(z)$ are shown numerically in Fig. 6. where we have fixed
$c_{0}=2.7,c_{1}=-0.7$ for different values of $\alpha$ and $\beta$. From
these figures we clearly see that we have a transition from a decelerated to
an accelerated universe around $z\approx 0.6$ which is compatible with
observations Planck1 ; Planck2 . However, not only at the present time
($z=0$), but also at the far future where $1+z\rightarrow 0$, the EoS
parameter, $\omega_{D}$, cannot cross the phantom line and we have always
$\omega_{D}>-1$.

Figure 6: The evolution of equation of state parameter $\omega_{D}$ and
deceleration parameter $q$ versus $1+z$ for GHDE1 model with GO cutoff.

###  III.2 GHDE2: The JBP model

Combining Eqs. (5) and (41) one can derive the evolution of dimensionless
GHDE2 density as

${\Omega}^{\prime}_{D}=(1-\Omega_{D})\left(\frac{2\Omega_{D}(1+z)^{4}}{\beta(c_{0}(1+z)^{2}+c_{1}z)^{2}}-\frac{2\alpha}{\beta}+3\right).$
(47)

The evolution of the dimensionless GHDE density parameter $\Omega_{D}$ as a
function of redshift $z$ is shown in Fig. 7.

Figure 7: The evolution of the dimensionless density parameter $\Omega_{D}$
versus $1+z$ for GHDE2 model with GO cutoff.

Using Eqs. (5), (42) and (43) we can obtain the equation of state and
deceleration parameters as

$\displaystyle\omega_{D}(z)$ $\displaystyle=$
$\displaystyle\frac{2\alpha}{3\beta\Omega_{D}}-\frac{2(1+z)^{4}}{3\beta[c_{0}(1+z)^{2}+c_{1}z]^{2}}-\frac{1}{\Omega_{D}},$
(48) $\displaystyle q(z)$ $\displaystyle=$
$\displaystyle-1-\frac{\Omega_{D}(1+z)^{4}}{\beta[c_{0}(1+z)^{2}+c_{1}z]^{2}}+\frac{\alpha}{\beta}.$
(49)

The behavior of the EoS parameter $\omega_{D}$ and the deceleration parameter
$q$ are plotted in Fig. 8. Again, our universe has a phase transition during
its history from a deceleration to an accelerated phase and $\omega_{D}$
cannot cross the phantom line even in the far future where $1+z\rightarrow 0$.

Figure 8: The evolution of equation of state parameter $\omega_{D}$ and
deceleration parameter $q$ versus $1+z$ for GHDE2 model with GO cutoff.

###  III.3 GHDE3: The Wetterich type

Using Eqs. (6) and (41) we obtain the evolution of dimensionless GHDE3 density
as

${\Omega}^{\prime}_{D}=(1-\Omega_{D})\left(\frac{2\Omega_{D}(1+c_{1}\ln(1+z))^{2}}{\beta
c_{0}^{2}}-\frac{2\alpha}{\beta}+3\right).$ (50)

The evolution of the dimensionless GHDE density parameter $\Omega_{D}$ as a
function of redshift $z$ is shown in Fig. (9).

Figure 9: The evolution of the dimensionless density parameter $\Omega_{D}$
versus redshift $1+z$ for GHDE3 model with GO cutoff.

Inserting Eq. (6) in Eqs. (42) and (43) we can obtain the EoS and the
deceleration parameters as follows

$\displaystyle\omega_{D}(z)$ $\displaystyle=$
$\displaystyle\frac{2\alpha}{3\beta\Omega_{D}}-\frac{2[1+c_{1}\ln(1+z)]^{2}}{3\beta
c_{0}^{2}}-\frac{1}{\Omega_{D}},$ (51) $\displaystyle q(z)$ $\displaystyle=$
$\displaystyle-1-\frac{\Omega_{D}[1+c_{1}\ln(1+z)]^{2}}{\beta
c_{0}^{2}}+\frac{\alpha}{\beta}.$ (52)

The behavior of the EoS parameter $\omega_{D}$ and deceleration parameter $q$
are plotted in Fig. 10. From these figures we see that, in contrast to the two
previous parametrization of $c(z)$, which $\omega_{D}$ cannot cross the
phantom line and $\omega_{D}>-1$ as $1+z\rightarrow 0$, here for the Wetterich
type of parametrization, the EoS parameter can cross the phantom line
($\omega_{D}<-1$) as $1+z\rightarrow 0$.

Figure 10: The evolution of equation of state parameter $\omega_{D}$ and
deceleration parameter $q$ versus $1+z$ for GHDE3 model with GO cutoff.

###  III.4 GHDE4: The Ma-Zhang type

Using Eqs. (7) and (41), it is easy to show that the evolution of the
dimensionless GHDE4 density can be obtained as

${\Omega}^{\prime}_{D}=(1-\Omega_{D})\left(\frac{2\Omega_{D}(1+z)^{2}}{\beta[c_{0}(1+z)+c_{1}\ln(2+z)-c_{1}(1+z)\ln
2]^{2}}-\frac{2\alpha}{\beta}+3\right).$ (53)

Fig. 9 shows the evolution of the dimensionless GHDE density parameter
$\Omega_{D}$. Again, at the early universe where $z\rightarrow\infty$ we have
$\Omega_{D}\rightarrow 0$, while at the late time where $z\rightarrow-1$, the
DE dominated, namely $\Omega_{D}\rightarrow 1$.

Figure 11: The evolution of the dimensionless density parameter $\Omega_{D}$
versus $1+z$ for GHDE4 model with GO cutoff

Inserting Eq. (7) in Eqs. (42) and (43) we can obtain the EoS and the
deceleration parameters as follows

$\displaystyle\omega_{D}(z)$ $\displaystyle=$
$\displaystyle\frac{2\alpha}{3\beta\Omega_{D}}-\frac{2(1+z)^{2}}{3\beta[c_{0}(1+z)+c_{1}\ln(2+z)-c_{1}(1+z)\ln
2]^{2}}-\frac{1}{\Omega_{D}},$ (54) $\displaystyle q(z)$ $\displaystyle=$
$\displaystyle-1-\frac{\Omega_{D}(1+z)^{2}}{\beta[c_{0}(1+z)+c_{1}\ln(2+z)-c_{1}(1+z)\ln
2]^{2}}+\frac{\alpha}{\beta}.$ (55)

The behavior of the EoS parameter $\omega_{D}(z)$ and the deceleration
parameter $q(z)$ are plotted in Figs. 12 and 13. From Fig. 12, we see that
when we fix $1\leq c_{0}\leq 2$ and $-1\leq c_{1}\leq 0$ and vary the
parameters $0\leq\alpha,\beta\leq 1$, the GHDE4 can behave as phantom DE model
at the future, while for $c_{0},c_{1}>0$ and $1<\alpha,\beta<2$ the EoS
parameter cannot cross the phantom line and is always larger than $-1$ (see
Fig. 13).

Figure 12: The evolution of equation of state parameter $\omega_{D}$ and
deceleration parameter $q$ versus $1+z$ for GHDE4 model with GO cutoff.

Figure 13: The evolution of equation of state parameter $\omega_{D}$ and
deceleration parameter $q$ versus $1+z$ for GHDE4 model with GO cutoff.

##  IV Conclusion and discussion

In this paper, we have studied HDE with time varying parameter $c^{2}$, the
so-called generalized holographic dark energy (GHDE), in a spatially flat
universe. It is important to note that, for the sake of simplicity, very often
the $c^{2}$ parameter in the HDE model is assumed constant. However, in
general, it can be regarded as a function of redshift parameter $z$ during the
history of the universe. By choosing four parameterizations of $c(z)$,
including the CPL type, JBP type, Wetterich type and Ma-Zhang type
parameterizations for $c(z)$, we have investigated the effects of varying
$c^{2}(z)$ term on the cosmological evolutions of GHDE model. As system’s IR
cutoff we have considered the Hubble radius $L=H^{-1}$ and the GO cutoff,
$L=(\alpha H^{2}+\beta\dot{H})^{{-1}/{2}}$. We have investigated the evolution
of EoS and deceleration parameters for all these parameterizations. We found
that in all GHDE models, with both Hubble and GO cutoffs, the universe has a
transition from a deceleration to an acceleration phase during its history. We
have found that for Hubble cutoff, the EoS parameter of GHDE can realize a
quintom behavior; namely, it evolves from a quintessence-like component to a
phantom-like component. While for the GO cutoff, not only at the present time
($z=0$), but also at the far future where $1+z\rightarrow 0$, the EoS
parameter, $\omega_{D}$, cannot cross the phantom line and we have always
$\omega_{D}>-1$ unless in a very special case.

It is worth mentioning that the simple and natural choice for the system’s IR
cutoff in the HDE model, is the Hubble radius $L=H^{-1}$. However, it was
argued that this choice for the IR cutoff leads to a wrong equation od state
for dark energy, namely $\omega_{D}=0$ Hsu , unless the interaction between
two dark components of the universe is taken into account Pavon . In this
paper, we demonstrated that by taking into account the time varying parameter
$c(z)$ can leads to an accelerated universe for $L=H^{-1}$ IR cutoff, even in
the absence of the interaction between two dark components of energy of the
universe. Besides, by suitably choosing of the parameter, not only the
accelerated universe can be achieved, but also the EoS parameter can cross the
phantom line $\omega_{D}=-1$, even in the absence of interaction. As far as we
know, this is a new result, which has not been reported already. In order to
more investigate the behavior of the EoS and deceleration parameters, we
plotted the evolution of these parameters versus redshift parameter $z$. From
these figures we see that has a decelerated phase at the early time
$(z\rightarrow\infty)$ and encounters a phase transition to an accelerated
phase around $z\approx 0.6$ which is consistent with recent observations Daly
; Kom1 ; Kom2 ; Planck1 ; Planck2 .

###### Acknowledgements.

We thank from the Research Council of Shiraz University. This work has been
supported financially by Research Institute for Astronomy & Astrophysics of
Maragha (RIAAM), Iran.

## References

  * (1) A.G. Riess, et al., Astron. J. 116 (1998) 1009;
S. Perlmutter, et al., Astrophys. J. 517 (1999) 565;
P. de Bernardis, et al., Nature 404 (2000) 955;
S. Perlmutter,et al., Astrophys. J. 598 (2003) 102.

  * (2) A. Cohen, D. Kaplan, A. Nelson, Phys. Rev. Lett. 82 (1999) 4971\.
  * (3) M. Li, Phys. Lett. B 603 (2004) 1.
  * (4) Q. G. Huang, M. Li, JCAP 0408 (2004) 013.
  * (5) S. D. H. Hsu, Phys. Lett. B 594 (2004) 13.
  * (6) E. Elizalde, S. Nojiri, S.D. Odintsov, P. Wang, Phys. Rev. D 71 (2005) 103504;
B. Guberina, R. Horvat, H. Stefancic, JCAP 0505 (2005) 001;
B. Guberina, R. Horvat, H. Nikolic, Phys. Lett. B 636 (2006) 80;
H. Li, Z. K. Guo, Y. Z. Zhang, Int. J. Mod. Phys. D 15 (2006) 869;
Q. G. Huang, Y. Gong, JCAP 0408 (2004) 006;
J. P. B. Almeida, J. G. Pereira, Phys. Lett. B 636 (2006) 75;
Y. Gong, Phys. Rev. D 70 (2004) 064029;
B. Wang, E. Abdalla, R. K. Su, Phys. Lett. B 611 (2005) 21.

  * (7) M. R. Setare, S. Shafei, JCAP 09 (2006) 011;
M. R. Setare, Phys. Lett. B 644 (2007) 99;
M. R. Setare, JCAP 0701 (2007) 023;
M. R. Setare, Phys. Lett. B 654 (2007) 1.

  * (8) B. Wang, Y. Gong and E. Abdalla, Phys. Lett. B 624 (2005) 141;
B. Wang, C. Y. Lin. D. Pavon and E. Abdalla, Phys. Lett. B 662 (2008) 1;
B. Wang, C. Y. Lin and E. Abdalla, Phys. Lett. B 637 (2005) 357.

  * (9) A. Sheykhi, Phys. Lett. B 681 (2009) 205 ;
A. Sheykhi, Class. Quantum Gravit. 27, (2010) 025007 ;
A. Sheykhi, Phys. Rev. D 84, (2011) 107302;
S. Ghaffari, A. Sheykhi and M. H. Dehghani, Phys. Rev. D 912, (2015) 023007;
S. Ghaffari, M. H. Dehghani and A. Sheykhi, Phys. Rev. D 89 (2014) 123009;
A. Sheykhi, M. Tavayef, Iran. J. Sci. Technol. Trans. Sci. DOI:
10.1007/s40995-016-0083-y.

  * (10) G. Hooft, gr-qc/9310026;
L. Susskind, J. Math. Phys. 36 (1995) 6377\.

  * (11) X. Zhang, F. Q. Wu, Phys. Rev. D 72 (2005) 043524;
X. Zhang, F. Q. Wu, Phys. Rev. D 76 (2007) 023502;
Q. G. Huang, Y.G. Gong, JCAP 0408 (2004) 006;
K. Enqvist, S. Hannestad, M. S. Sloth, JCAP 0502 (2005) 004;
J. Y. Shen, B. Wang, E. Abdalla, R.K. Su, Phys. Lett. B 609 (2005) 200.

  * (12) P. Horava and D. Minic, Phys. Rev. Lett. 85, 1610 (2000) [hep-th/0001145];
S. D. Thomas, Phys. Rev. Lett. 89, 081301 (2002).

  * (13) W. Fischler and L. Susskind, arXiv:hep-th/9806039;
R. Bousso, J. High Energy Phys. 9907, 004 (1999).

  * (14) S. Nojiri and S. D. Odintsov, Gen. Relativ. Gravit. 38, 1285 (2006).
  * (15) C. J. Gao, X. L. Chen, and Y. G. Shen, Phys. Rev. D 79, 043511 (2009);
R. G. Cai, B. Hu, and Y. Zhang, Commun. Theor. Phys. 51, 954 (2009).

  * (16) D. Pavon, W. Zimdahl, Phys. Lett. B 628 (2005) 206;
W. Zimdahl and D. Pavon, Class. Quantum Grav. 24 (2007) 5461\.

  * (17) N. Radicella, D. Pavon, JCAP 10, 005 (2010).
  * (18) M. Malekjani, R. Zarei, M. Honari-Jafarpour, Astrophys Space Sci. 343, 799 (2013), [arXiv:1209.1089],
  * (19) M. Chevallier and D. Polarski, Int. J. Mod. Phys. 10, 213 (2001);
E.V. Linder, Phys. Rev. Lett. 90, 091301 (2003).

  * (20) H. K. Jassal, J. S. Bagla and T. Padmanabhan, MNRAS. 356, 11 (2005).
  * (21) C. Wetterich, Phys. Lett. B 594, 17 (2004); Y. G. Gong, Class. Quantum Grav. 22, 2121(2005).
  * (22) J. Z. Ma and X. Zhang, Phys. Lett. B 699, 233 (2011); H. Li and X. Zhang, Phys. Lett. B 703,2 (2011).
  * (23) L. N. Granda and A. Oliveros, Phys. Lett. B 671, 199 (2009).
  * (24) R.A. Daly et al., Astrophys. J. 677, 1 (2008).
  * (25) E. Komatsu et al. [WMAP Collaboration], Astrophys. J. Suppl. 192, 18 (2011).
  * (26) V. Salvatelli, A. Marchini, L. L. Honorez and O. Mena, Phys. Rev. D 88, 023531 (2013).
  * (27) J.Q. Xia, H. Li, X. Zhang, Phys. Rev. D 88, 063501 (2013).
  * (28) W. Yang, L. Xu, Phys. Rev. D 89, 083517 (2014) ;
B. Novosyadlyj, O. Sergijenko, R. Durrer, V. Pelykh, DOI:
10.1088/1475-7516/2014/05/030, [arXiv:1312.6579].