MATHEMATICAL FORMULAE AND NUMERICAL METHODS

The governing equations The continuity, momentum and energy conservation laws are used as the governing equations for fluid flows. A Newtonian fluid is considered in this study. The Boussinesq

approximation for the buoyancy due to the density variation of a fluid is employed. The governing equations for an incompressible Newtonian fluid are expressed in the following non-dimensional forms:

(1)

(2)

(3)

where u, p and are non-dimensional velocity, pressure and temperature, respectively. The dimensionless temperature is defined as , where is the wall/body temperature and is a reference temperature. In Eq. (2), f represents the virtual force only applied to solids. is the buoyancy with the direction of gravity being . In Eq. (3), represents the virtual heat source term only applied to solids to achieve desired temperature distribution of the solid object which will further imposes wanted thermal boundary conditions on fluids. The dimensionless parameters , and in Eqs. (2) and (3) are defined differently later depending on whether the convection is forced or natural. For forced convection without buoyancy, = 1/Re, = 0 and = 1/RePr where Re and Pr are Reynolds number and Prandtl number, respectively. For natural convection in which buoyancy is important, = Pr, = RaPr and = 1 where Ra

is Rayleigh number, or we have , , and , where is the

Rayleigh number.

Numerical procedures A staggered grid arrangement is utilized in this study. That is, pressure and temperature are located at the center of the computational cell while the velocities are placed at the faces of the cell. We use the time-splitting schemes to advance velocity and temperature in Eqs. (2) and (3). First, the velocity and temperature are advanced from the time level to the first intermediate level by solving the advection-diffusion equations without pressure and virtual force f in Eq. (2) and the virtual heat source q in Eq. (3). This step can be stated as the following

(4)

(5)

where and discretized by the second-order upwind scheme for

the convection term and central difference schemes for the diffusion term. Subsequently, we use the second-order Adam-Bashford method for the temporal discretization of those equations in such away that term and

can be approximated by

(6) and

(7)

Many conventional projection methods use Crank-Nicolson method for the diffusion term forming an implicit scheme, which has the advantage of allowing a larger time step, compared with the current explicit scheme, when Reynolds number is small. However, since the Reynolds number of many interesting heat convection problems is often moderate or high and therefore the time step constraint is dominated by CFL condition, here we use the explicit scheme stated above to save overhead cost of solving Helmholtz equation at each time step due to Crank-Nicolson scheme.

The first intermediate velocity in Eq. (4), in general, does not satisfy the divergence-free condition in Eq.

(1). At the second step, we advance the first intermediate velocity by the pressure term

(8) and apply the divergence to both sides. Eq. (8) becomes

(9) By imposing

(10) We achieve the pressure Poisson equation

(11)

Once Eq. (11) is solved, we can advance from to . Conventionally, this would be the end of projection method and actually . It is so indeed for a pure fluid cell ( ). However, for a cell fully or

11

solid that is known in advance. This is accomplished by the virtual force applied only to solid part and stated as follows

(12)

Likewise, we need to compensate the solid temperature to a desired distribution by the virtual heat source term and this gives the following step in advancing to :

(13)

Note that in Eqs. (12) and (13) and automatically for a pure fluid cell ( ). For a pure solid cell, ( ), and are known in advance and the virtual force and virtual heat source can be computed reciprocally. The integration of virtual force will be a good approximation of the resultant force exerted on a non-accelerated solid object via viscous and pressure gradient physically,

(14)

The thermal boundary conditions on a solid surface are and for hot and cold surfaces, respectively (see Dhole et al. 2006). Once the velocity and temperature fields are obtained, the local Nusselt number on the solid surface is evaluated using

(15)

where is the normal direction of the solid surface. Afterward, surface-average Nusselt number becomes (16)

where is the total surface area of the solid object .

In summary, the numerical procedure of the current method at each time step is given below:

1.Determine through the position and orientation of the solid rigid object.

2.Calculate and via Eqs. (4) and (5)).

3.Solve the Poisson equation (Eq. (11)) and advance to via Eq. (8).

4.Update the solid velocity to the prescribed value and compute the virtual force required.

5.Update the solid temperature to the prescribed value and compute the virtual heat source required.

The scheme is basically 2nd order accurate in time and space. However, an extra error may be introduced at a fluid/solid boundary and may degrade the total accuracy to be super-linear, unless we interpolate the fluid velocity to fit the exact boundary condition exerted by the solid object at the boundary, which will generate non-physical internal flow inside the solid object. We did not do the interpolation here to maintain the 2nd order accuracy in space for the following reason. Herein, by Eq. (12), we let the velocity in each solid cell (

=1) satisfy the actual rigid-body motion felt by the solid object and the velocity of a cut cell ( ) satisfy the weighted average (by ) of solid and fluid velocity to compute the virtual force in each cell, which summarizes to obtain the total force exerted on a solid object easily. The price of doing this is that it causes extra boundary error and sacrifices the 2nd order accuracy in space to be a super-linear one (see Fadlun et al.

[2000]). The same philosophy is applied to the temperature and the virtual heat sources in the energy equation, which sum to the physical Nusselt number.

The last step, Eq. (12), in the projection method to uphold the physical solid velocity due to the rigid-body motion may destroy the divergence-free condition in the cells in which the fluid/solid interface is located.

This implies non-physical mass sources or sinks may be generated in such cells. Fortunately, these minute sources and sinks come in pairs and end up canceling each other to maintain an overall mass conservation in a weak sense. This explains why the scheme still performs well in the benchmarking computations.