## Abstract

A novel high-temperature particle solar receiver is developed using a light trapping planar cavity configuration. As particles fall through the cavity, the concentrated solar radiation warms the boundaries of the receiver and in turn heats the particles. Particles flow through the system, forming a fluidized bed at the lower section, leaving the system from the bottom at a constant flowrate. Air is introduced to the system as the fluidizing medium to improve particle heat transfer and mixing. A laboratory scale cavity receiver is built by collaborators at the Colorado School of Mines and their data are used for model validation. In this experimental setup, near IR quartz lamp is used to provide flux to the vertical wall of the heat exchanger. The system is modeled using the discrete element method and a continuum two-fluid method. The computational model matches the experimental system size and the particle size distribution is assumed monodisperse. A new continuum conduction model that accounts for the effects of solid concentration is implemented, and the heat flux boundary condition matches the experimental setup. Radiative heat transfer is estimated using a widely used correlation during the post-processing step to determine an overall heat transfer coefficient. The model is validated against testing data and achieves less than 30% discrepancy and a heat transfer coefficient greater than 1000 W/m^{2} K.

## 1 Introduction

Concentrating solar power has been shown to be a promising source of continuous, dispatchable renewable energy, especially when combined with thermal energy storage (TES) systems [1,2]. TES traditionally works on the principle of heating a storage media in the charging phase and discharging that heat when needed to ensure a more consistent supply of energy [3,4]. Previous approaches for TES include using molten salts which circulate through the solar receiver where they absorb the concentrated solar radiation in the form of heat during the charging cycle and flow through heat exchangers for power generation in the discharging cycle [5]. Though molten salts have been shown to be an effective heat transfer medium, such salts can become chemically unstable at high temperatures and can also become corrosive [6]. In recent years, high operating temperature particle receivers have been designed. The main benefits of particle receivers are that the particulates are relatively inexpensive, and that the receivers can operate at high temperatures, enabling high thermodynamic efficiencies [7]. Particle receivers of various configurations such as cyclonic, rotary, and free falling as well as with fixed or fluidized bed configurations have been developed, each with their own unique set of benefits [8,9,3]. Lately, Brayton cycles have been shown to be more efficient than the traditionally used Rankine cycle. As a result of this, moving bed heat exchangers have become ubiquitous in systems where particles are used for thermal energy storage. Moving packed bed supercritical carbon dioxide heat exchangers have gained significant popularity due to the significant increase in efficiency of the $sCO2$ power cycles or recompression closed Brayton cycles [10,11]. Studies have shown that bubbling fluidized bed heat exchangers help obtain a significantly higher heat transfer coefficient as compared to moving packed beds [11]. In this study, we look at an inert solid particle receiver where the particles are used as the TES media and are directly heated in the receiver cavity by the concentrated solar radiation.

The National Renewable Energy Laboratory (NREL) is developing a high-temperature particle solar receiver using a light trapping planar cavity solar configuration. The receiver design consists of an array of light trapping receiver cavities arranged in the form of a cylinder. Light enters the receiver tower system, concentrated by a heliostat field, and gets trapped within the cavity. This geometry has been shown to have high thermal efficiencies with low reradiation losses. The planar cavity concept proposed for high-temperature operations and high thermal efficiencies of the receiver has been described by Ma and Martinek [12]. As mentioned, the particles are indirectly heated as they fall through the receiver cavity and are used as the TES medium. This heating process is integral to the performance of the receiver and hence significant testing and modeling efforts have been allocated to optimize the performance of this heat exchanger. The flow of the particles is gravity driven in the net downward direction at a fixed flowrate through the receiver cavity. A fluidized bed forms within the receiver cavity, which helps the wall to particle heat transfer by increasing the residence time and the mixing.

This study focuses on the development and validation of a computational fluid dynamics (CFD) model of the heat exchanger. The model is calibrated and validated against experimental data using a relatively novel approach. The aim of the model is to elucidate the heat transfer and fluid dynamics in the system and help optimize the performance of the receiver. There are two main motivations to conduct these simulations: (i) to study the particle and gas counterflow motion and study the effect of varying the particle or airflow rate on the heat transfer coefficient and (ii) to enhance the heat transfer.

This study follows the following structure:

Model details

- —
Two-fluid model (TFM)

- —
Discrete element modeling (DEM)

- —
Details of the experimental system modeled

- —
Batch mode: used for validation with HSP 40/70 Carbobead particles

- —
Continuous mode: counterflow between air and particles with silica sand

- —
Model validation: batch mode

- —
DEM validated against experimental data

- —
Nusselt number correlation by Morris et al. [13] against validated DEM model

- —
Two-fluid model using the Nusselt number correlation tuned to match experimental data

- —
Observations, results, and parametric study: continuous mode

## 2 Methods

### 2.1 Computational Details.

Multiphase flow with interphase exchange, MFiX, an open source CFD code for multiphase flows, is used as the modeling framework for the simulations. The TFM [14] is used to model the particle and air counterflow heat exchanger portion of the receiver, as well as the multiphase heat transfer in the system. The TFM is an Eulerian–Eulerian model, i.e., both the fluid as well as the particles are modeled as interpenetrating continuum fluids [14]. Since the solid particles are not individually tracked, the model is computationally inexpensive, but requires closures for the constitutive relations that dictate the solid stresses and heat flux [13]. The governing equations for mass, momentum, and energy that are solved numerically are shown in Table 1. The equations in Table 1 are obtained from the MFiX theory guide [14].

Continuity equations |
---|

Gas: |

$\u2202(\epsilon g\rho g)\u2202t+vg\u22c5\u2207(\epsilon g\rho g)=0$ |

Solid: |

$\u2202(\epsilon s\rho s)\u2202t+vs\u22c5\u2207(\epsilon s\rho s)=0$ |

Momentum equations |

Gas: |

$\u2202(\epsilon g\rho gvg)\u2202t+\u2207\u22c5(\epsilon g\rho gvgvg)=\u2212\epsilon g\u2207Pg+\u2207\u22c5\tau g\xaf\xaf\u2212\beta g,s(vs\u2212vg)\u2212\epsilon g\rho gg$ |

Solid: |

$\u2202(\epsilon s\rho svs)\u2202t+\u2207\u22c5(\epsilon s\rho svsvs)=\u2212\epsilon s\u2207Ps+\u2207\u22c5\tau s\xaf\xaf\u2212\beta g,s(vs\u2212vg)+\epsilon s\rho sg$ |

Energy equations |

Gas: |

$\u2202(\epsilon g\rho gCp,gTg)\u2202t+vg\u22c5\u2207(\epsilon g\rho gCp,gTg)=\u2212\u2207\u22c5qg\u2212Hg(Tg\u2212Ts)$ |

Solid: |

$\u2202(\epsilon s\rho sCp,sTs)\u2202t+vs\u22c5\u2207(\epsilon s\rho sCp,sTs)=\u2212\u2207\u22c5qs+Hg(Tg\u2212Ts)$ |

Interphase heat transfer |

$Hg=6\kappa g\epsilon sNugsdp2Nugs=(7\u221210\epsilon g+5\epsilon g2)(1+0.7Re0.2Pr1/3)+(1.33\u22122.4\epsilon g+1.2\epsilon g2)Re0.7Pr1/3$ [15] |

Continuity equations |
---|

Gas: |

$\u2202(\epsilon g\rho g)\u2202t+vg\u22c5\u2207(\epsilon g\rho g)=0$ |

Solid: |

$\u2202(\epsilon s\rho s)\u2202t+vs\u22c5\u2207(\epsilon s\rho s)=0$ |

Momentum equations |

Gas: |

$\u2202(\epsilon g\rho gvg)\u2202t+\u2207\u22c5(\epsilon g\rho gvgvg)=\u2212\epsilon g\u2207Pg+\u2207\u22c5\tau g\xaf\xaf\u2212\beta g,s(vs\u2212vg)\u2212\epsilon g\rho gg$ |

Solid: |

$\u2202(\epsilon s\rho svs)\u2202t+\u2207\u22c5(\epsilon s\rho svsvs)=\u2212\epsilon s\u2207Ps+\u2207\u22c5\tau s\xaf\xaf\u2212\beta g,s(vs\u2212vg)+\epsilon s\rho sg$ |

Energy equations |

Gas: |

$\u2202(\epsilon g\rho gCp,gTg)\u2202t+vg\u22c5\u2207(\epsilon g\rho gCp,gTg)=\u2212\u2207\u22c5qg\u2212Hg(Tg\u2212Ts)$ |

Solid: |

$\u2202(\epsilon s\rho sCp,sTs)\u2202t+vs\u22c5\u2207(\epsilon s\rho sCp,sTs)=\u2212\u2207\u22c5qs+Hg(Tg\u2212Ts)$ |

Interphase heat transfer |

$Hg=6\kappa g\epsilon sNugsdp2Nugs=(7\u221210\epsilon g+5\epsilon g2)(1+0.7Re0.2Pr1/3)+(1.33\u22122.4\epsilon g+1.2\epsilon g2)Re0.7Pr1/3$ [15] |

A model for solid shear developed by Johnson and Jackson [16] is used to describe the “solid shear” which transitions between the plastic and viscous regimes at a critical packing fraction. The Gidaspow drag model [17] is used to model fluid–solid momentum interactions.

The variables $ew$ and $ep$ are the emissivity coefficients of the wall and particle cluster, respectively. A median value from the acceptable range [20] is chosen, and this value affects the heat transfer coefficient by $\xb130W/m2K$.

The continuum conduction model developed by Morris et al. [13] was developed for moving beds and gravity-assisted chute flows. To test the validity of this Nusselt number correlation for fluidized beds and heat transfer from a vertical wall, DEM is performed for bubbling fluidized beds. DEM is an Eulerian–Lagrangian method where the fluid is modeled as a continuum media similar to the TFM while the solid phase is studied as discrete particles and their motions and momenta are individually tracked. MFiX-DEM is used as a solver for the simulations. The solid phase collisions are resolved using the soft sphere approach which is a kind of a spring-dashpot model [21]. The particle spring constant is set to 1000 N/m and the particle–particle and particle wall friction are set to 0. The gas and solid phases are fully coupled using interphase coupling equations, the Gidaspow drag model for momentum coupling [17], and the Gunn model for energy coupling [15]. The gas phase equations are solved similar to the continuum model while the solid phase equations are solved by resolving the collisions using Newton’s laws of motion. The Rong and Horio [22] model is used for conduction between particles and between particles and walls, through direct contact or more dominantly, via fluid film. Soft sphere corrections are applied by setting Young’s modulus to 70 GPa and the Poisson’s ratio to 0.17. This thickness of the fluid film around the particle through which conduction occurs is called the lens radius. The minimum conduction distance ($s$) is a model parameter that can be approximated well as the particle surface roughness. Both the lens radius and the minimum conduction distance have a significant effect on the heat transfer to the particles as shown by Morris et al. [13].

### 2.2 Test System and Modeling Domain.

A heat exchanger test station designed by the Colorado School of Mines and their data are used for the validation of the MFiX models. In this heat exchanger, there is a counterflow in the motion between air and particles. Particles enter the system from the top, form a bed, and traverse in a net downward motion to leave the system from the bottom at a constant flowrate. Air enters the system through a distributor near the bottom of the cavity, and lightly fluidizes the system to aid with the mixing and heat transfer. The air leaves the system from outlets present on the side walls near the top. A near IR quartz lamp provides heat flux to one of the vertical walls. Hence the energy input to the system is through a heated wall. Figure 1 shows the test station along with particle and air flow directions as well as location of thermocouples. Table 2 summarizes the particle properties and model inputs for both batch mode and continuous flow simulations. Particles do not exit the system during batch tests, whereas particles flow at a constant flowrate through the system for continuous mode tests. The experimental data used for model validation was obtained from the researchers at the School of Mines. The data are submitted and in review for publication.

Batch mode | HSP 40/70 (Carbobead) | Continuous mode | Silica sand |
---|---|---|---|

$dp$ | $360\mu m$ | $dp$ | $150\mu m$ |

$\rho p$ | $3620kg/m3$ | $\rho p$ | $2700kg/m3$ |

$Cp,p$ | $1150J/kgK$ | $Cp,p$ | $1155J/kgK$ |

$ep$ | 0.85 | $ep$ | 0.85 |

Batch mode | HSP 40/70 (Carbobead) | Continuous mode | Silica sand |
---|---|---|---|

$dp$ | $360\mu m$ | $dp$ | $150\mu m$ |

$\rho p$ | $3620kg/m3$ | $\rho p$ | $2700kg/m3$ |

$Cp,p$ | $1150J/kgK$ | $Cp,p$ | $1155J/kgK$ |

$ep$ | 0.85 | $ep$ | 0.85 |

### 2.3 Batch Mode Testing.

For batch mode testing, there is a fixed mass of particles, so the system acts as a fluidized bed with heat addition. Heat flux is incident on one of the vertical walls (front) and the other wall (back) is cooled with fans to keep the overall temperature of the bed quasi-steady and with small thermal gradients. The TFM model is calibrated for one value of constant bed temperature and the calibration is validated using a different bed temperature. The particles are preheated to this bed temperature before starting the experiment and hence in the models, the particles are initialized to have the bed temperature. Particles tested in batch mode are $360\mu m$ HSP 40/70 Carbobeads. The inlet air has the same temperature as the bed and a wide range of fluidization conditions are tested. The experimental data from School of Mines are used for model validation.

A DEM model is setup for batch mode which matches the test station dimensions in the direction of the heat transfer. The particles are monodisperse spheres with diameters equal to the Sauter mean diameter, $360\mu m$, and the other physical properties are matched. The heated wall is modeled as a constant temperature wall ($362\u2218C$ for a $300\u2218C$ bed bulk temperature) and all other walls are considered to be adiabatic.

Figure 1 shows a schematic of the test station in batch mode. The batch mode test data are used to validate and calibrate both the DEM and TFM conduction models. For the validation of the TFM model, the actual size of the test station ($0.1\xd70.25\xd7$$0.012m$) is simulated to eliminate challenges and uncertainties associated with system scaling. The heat fluxes are matched using the measured flux value as the boundary condition ($25kW/m2$). Convection cooling is modeled using a negative heat flux boundary condition.

### 2.4 Continuous Mode Simulations.

After the model is validated by the batch mode experimental results, TFM simulations with particle flow through the system are performed which mimic the working of the actual receiver cavity. The same dimensions as the School of Mines experimental test station are used. The particle and air counterflow is studied and parametric studies are conducted. Figure 2 shows the computational domain for the counterflow system. The front view is on the far left and shows the initial particle bed and freeboard region. The side view shows the particle bed, freeboard region, and the air outlet (top block on the side view). On the far right, the bottom view is shown followed by the top view of the domain. The bottom includes the air inlet (middle block, bottom view) as well as particle outlets (top and bottom blocks, bottom view) and the top shows the particle inlet (top view). The particles enter the domain from the top ($+y$) surface at a set mass flowrate and leave the domain from a portion of the bottom surface ($\u2212y$). Air enters the domain from the bottom ($\u2212y$), between the two particle outlets, and leaves from the top regions in the $+/\u2212x$ directions. The side walls ($+/\u2212x$) and the back surface ($\u2212$*z*) are no-slip adiabatic walls.

The front wall (+*z*) acts as the receiver surface and receives a constant heat flux. To aid with maximizing the heat transfer to the particles, smaller $150\mu m$ silica sand particles are tested because reducing the particle size has been shown in the literature to enhance the heat transfer [23]. Note that the continuous simulations are performed with only the two-fluid model.

## 3 Results and Discussion

### 3.1 Model Validation

#### 3.1.1 Validation of the Effective Conductivity Model.

Since the Nusselt number correlation by Morris et al. [13] was obtained from chute flows, the validity of the model for fluidized beds is tested. For this, the DEM results are validated against the experimental batch mode results and this validated DEM is used to obtain Nusselt numbers to validate the correlation.

The DEM model is setup to match the experimental domain in the direction of heat transfer and is scaled down in the other dimensions to reduce computational expense. The heated wall is set to a constant temperature, which is an average of the observed wall temperature and all other walls are setup as adiabatic boundaries. The time scales of the DEM simulations are quite small so the bed is mostly homogeneous at a temperature equal to the initial bed temperature, which in this case is $300\u2218C$. DEM is performed in batch mode for the $360\mu m$ Carbobead particles.

To validate the DEM models, the simulated heat transfer coefficient is compared to the experimentally obtained values. The heat transfer coefficient ($htotal$) is calculated as the sum of the conductive and convective ($hc$) and radiative ($hr$) heat transfer coefficients. This heat transfer coefficient is plotted against the ratio of the superficial velocity ($u$) and minimum fluidization velocity ($umf$). This minimum fluidization velocity is the minimum superficial velocity at which the weight of the particle bed is completely supported by the pressure drop of the fluid across the bed and is calculated using the Ergun equation [24].

Figure 3 compares the DEM heat transfer coefficients which are calculated as shown in Eq. (9) with the experiment, validating the model. This validated model is then used to obtain the Nusselt number near the heated wall. The Nusselt number obtained is scaled by the maximum Nusselt number that is obtained similar to the method followed by Morris et al. and plotted against the solid fraction of the cell. The results obtained are compared with the correlation by Morris et al. [13]. This also helps validate the inflow of heat to the system from the heated wall.

Figure 4 shows the comparison of the Nusselt number obtained from the validated DEM simulations and the Nusselt number correlation by Morris et al. [13]. The *y* axis on the left shows the ratio of the Nusselt number and the maximum Nusselt number and this is plotted against the solid fraction of each near-wall cell. These plots obtained from DEM agree well with the correlation by Morris et al. [13] at more fluidized conditions with small deviations at lesser fluidized cases. This helps validate the solid conductivity model for bubbling beds, which is then used in the two-fluid model.

#### 3.1.2 Validation of the Two-Fluid Model.

It is important to note that the calibrated model does not perform well in simulating heat transfer in the dense packed bed regime. It is suspected that the reason for this is the hydrodynamic predictions which are a known issue with the continuum model [26,27] and further studies are needed to ensure the validity of the model in those regimes.

### 3.2 Continuous System With Particle/Air Counterflow.

With the complete receiver system, the aim of the simulations is to model the receiver closer to its physical operating conditions and obtain an optimum regime to maximize the heat transfer to the particles. To this effect, smaller particles, silica sand with mean diameter of $150\mu m$, are tested because the literature shows that reducing particle size generally tends to increase the heat transfer coefficient [23]. These smaller particles are likely what will be used in the actual receiver and hence it is important to study the heat transfer through these particles. It is noted that the model predicted heat transfer coefficient values are likely larger than the actual case because of the assumption of smooth spheres. The actual values would likely be less depending on the surface roughness. Although not shown, simulations were run with the assumption that silica sand has the same surface roughness as the Carbo beads. Heat transfer coefficient values greater than $1000W/m2K$ are obtained for nearly all cases.

In the continuous system, heating takes place only from one wall which acts like the receiver wall. Due to this, the region of the bed near the hot wall is generally warmer than the rest of the domain. In these simulations, especially with smaller particles and lower regimes of fluidization, we observe that the gas velocity is higher near the heated wall which indicates a weak channel forms near the wall as shown in Fig. 7. Figure 7 shows the gas velocity contours on the side wall of the domain, where the scale is such that red is a faster $+y$ gas velocity. This channeling reduces the particle concentration near the wall which in turn reduces heat transfer from the wall to the particles. The channeling effect has been attributed to density reduction of air in the warmer regions as compared to the cooler regions. The reduction in density results in a higher gas velocity. Small changes in concentration have a significant effect on the thermal contact, especially as the concentration approaches a packed state. Hence, a proposed solution to eliminate the channeling is to slightly tilt the entire domain toward the heated wall, so that the particles are assisted by gravity to have better thermal contact with the hot wall. The right axis of Fig. 4 also shows how the solid conductivity is affected by the solid concentration.

Regular operating conditions would generally place our solid volume fraction between 0.45 and 0.58 where even slight changes have a significant effect on the effective solid conductivity and hence the heat transfer.

To study the effect of small inclination angles using MFiX simulations, $150\mu m$ silica sand is used, the particle flowrate is kept constant, and the particles enter the system at $750\u2218C$. Air enters the system at a velocity of $0.035m/s$ and $750\u2218C$ as well. In Fig. 8, it can be seen that even a small inclination significantly helps overcome the channel and improves the heat transfer coefficient.

A parametric study is performed to obtain the optimal particle and airflow rates. To study the air flowrate, $150\mu m$ silica sand is used and once again the particle flowrate is kept constant at $0.024kg/m3$ and the particles enter the system at 750 $\u2218C$. For simplicity and to better compare to prior simulations and experimental data, the domain is not inclined. Air enters at $750\u2218C$ and the air flowrate is varied. The inlet flux, $q\u2033$, is kept constant at $200kW/m2$. Figure 9 shows the results. The minimum fluidization velocity, $umf$, is defined for the particles at $750\u2218C$ using the Ergun equation [24]. This does not account for the downward motion of air which would change the effective superficial velocity slightly. It can be observed that fluidization significantly improves the mixing in the bed and hence has a large effect on the heat transfer coefficient. The minimum fluidization velocity of the bed is $0.0156m/s$. As the air flowrate increases, a slight peak is observed, similar to the batch mode results. Figure 10 shows that the heat transfer coefficient reduces as the mass flowrate of particles through the system is increased. This reduction of the heat transfer coefficient for increasing flowrate occurs because the increase in mass flowrate would reduce the residence time of the particles in the heat exchanger. Overall, the configuration of the heat exchanger seems to have good performance with heat transfer coefficients of over a $1000W/m2K$ which puts it well above the industry average and will hence be very effective in improving the receiver efficiency.

## 4 Conclusion

A high-temperature particle solar receiver is being developed which uses a planar cavity light trapping configuration. A novel narrow channel fluidized bed heat exchanger [11] is being tested to be implemented in the receiver cavity. The heat exchanger has a net gravity-aided downward motion of the particles which are in counterflow against air which acts as a fluidizing medium. This serves the benefit of increasing the residence time of the particles within the heat exchanger as well as mixing the particles for better heat transfer.

A Lagrangian–Eularian (DEM) model was built and validated and then used to validate the solid conduction model. A continuum model is built to simulate this heat exchanger and calibrate against test results. The calibration is conducted by comparing the model heat transfer coefficient to experimentally obtained values and the surface roughness of the particles is varied from completely smooth to a more realistic value. The calibrated model is then validated with a different set of experimental results. The calibration and validation are performed for batch mode with no particle flow and larger HSP 40/70 particles. Smaller silica sand particles are simulated in the continuous mode to enhance heat transfer to the particles. A few important conclusions can be obtained from the continuous simulations. Even slight differences in solid concentration near the heated wall can significantly affect the heat transfer to the particles. The simulations have a heated wall and hence a region where the density of air reduces causing an increase in air velocity which leads to the formation of a narrow channel near the heated wall. This channel reduces the thermal contact of particles with the heated wall. Using an incline to gravity assist particles toward the heated wall helps significantly increase the heat transfer from the wall to the particles.

Fluidization of the bed helps in better heat transfer but increasing the amount of fluidization does not make much difference to the heat transfer until the slugging stage is reached when the heat transfer reduces due to large bubbles. A lower particle flowrate leads to more heated particles as a direct result of the residence time of the particles within the system.

## Acknowledgment

The authors acknowledge individuals, institutions, or companies that supported them in preparing the work. Those mentioned might have provided technical support, insightful comments, or conversations, materials used in the work, or access to facilities.

## Funding Data

This work was sponsored by the National Renewable Energy Laboratory which is granted by the U.S. Department of Energy award No. DE-EE0038896.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## Nomenclature

### Symbols

- $e$ =
thermal emissivity

- $g\u2192$ =
gravity $(m/s2)$

- $h$ =
heat transfer coefficient $(W/m2K)$

- $q\u2192$ =
heat transfer rate $(W)$

- $t$ =
time $(s)$

- $v\u2192$ =
velocity vector $(m/s)$

- $D$ =
diameter (m)

- $H$ =
thermal conductance $(W/K)$

- $P$ =
pressure $(atm)$

- $T$ =
temperature $(K$ or $\u2218C)$

- $Cp$ =
specific heat $(J/kgK)$

- $q\u2033$ =
heat flux $(W/m2)$

- Nu =
Nusselt number

- $\beta $ =
gas–solids drag coefficient

- $\epsilon $ =
volume fraction

- $\kappa $ =
thermal conductivity $(W/mK)$

- $\rho $ =
density $(kg/m3)$

- $\tau $ =
shear stress $(N/m2)$

### Subscripts

## References

_{2}O and CO

_{2}

_{2}Heat Exchangers in Next Generation CPS Plants

*The Distinct Element Method as a Tool for Research in Granular Media: Report to the National Science Foundation Concerning NSF Grant ENG75-20711*, Vol. 2, Department of Civil and Mineral Engineering, Institute of Technology, University of Minnesota, Minneapolis, MN.