An axisymmetric finite volume formulation for the solution of heat conduction problems using unstructured meshes

TECHNICAL PAPERS

An axisymmetric finite volume formulation for the solution of heat conduction problems using unstructured meshes

P. R. M. Lyra^I; R. de C. F. de Lima^I; D. K. E. de Carvalho^I; G. M. L. L. da Silva^II

^IDepartamento de Engenharia Mecânica – UFPE; Av. Acadêmico Hélio Ramos, S/N; 50740-530 Recife, PE. Brazil; [e-mail]; [e-mail]; [e-mail]

^IIDepartamento de Energia Nuclear – UFPE; Av. Prof. Luiz Freire, 1000; 50740-540 Recife, PE. Brazil; [e-mail]

ABSTRACT

In this work, a finite volume formulation developed for two-dimensional models is extended to deal with axisymmetric models of heat conduction applications. This formulation uses a vertex centered finite volume method and it was implemented using an edge-based data structure. The time and domain discretization using triangular meshes is described in details, including the treatment of boundary conditions, source terms, and domains with multiple materials. The proposed formulation is validated and proves to be effective and flexible through the solution of simple model problems.

Keywords: FVM, unstructured meshes, axisymmetric models, heat transfer applications

Introduction

Finite Volume Methods (FVM), are especially attractive for the solution of conservation laws problems. The construction of FVM methods capable to deal with unstructured meshes can be of utmost importance in order to handle applications that involve complex geometries (AGARD Report 787, 1992). The node-centered FVM using an edge-based data structure, which was presented in Lyra et al (2004a), is very flexible to deal with control volumes of different shapes associated to generic unstructured meshes.

When studying complex problems, bi-dimensional models may become quite limited when studying complex problems. On the other hand, a fully three dimensional model may lead to heavy computational requirements which may sometimes turn the analysis unfeasible. In many situations the tri-dimensional modeling may be avoided due to the axisymmetric behavior of the problems. Therefore, it is interesting to develop a FVM formulation suitable to deal with axisymmetric models, where the plane discretization is rotated 2p radians around the axis of symmetry. In the present paper a vertex centered finite volume formulation using median dual control volumes is developed and implemented for the simulation of axisymmetric heat conduction problems subjected to different types of boundary conditions (Dirichlet, Neumann, and Cauchy) and to some non-conventional loads, considering also systems with multi-materials. The temporal discretization is done by a simple first order Euler-forward finite difference method.

The use of an edge-based data structure is very effective in terms of reducing: indirect addressing to retrieve local information required by the solver; CPU time and memory requirements (Barth, 1992; Peraire et al., 1993; Sorensen et al., 2001; Lyra, 1994a, b), particularly when three-dimensional models are adopted.

In the next section we briefly describe the physical and mathematical models addressed. Then a detailed description of the finite volume formulation proposed for such class of models is fully developed and validated by considering the solution of simple model problems which allow to illustrate the versatility and effectiveness of the whole procedure. Finally, some conclusions and the present status of this research are presented.

Mathematical Model

As mentioned previously, an axisymmetric heat conduction model is adopted, where the z-axis is the axis of symmetry. Then it is convenient to formulate the problem using a cylindrical coordinate system. Since all coefficients are independent of q, the temperature distribution is a function of (r, z). The heat conduction equation is then written as (Ösizik, 1980):

where, r is the mass density, c is the specific heat, T is the temperature, and Q represents the source (or sink) terms. The spatial domain of the problem is represented by W, with r being the radial coordinate and z, the axial one. The time interval of integration is represented by , with tⁱ and t^f, being the initial and the final time stage. For simplicity, the medium is considered orthotropic.

The heat flux is a function of the temperature gradient, and is modeled by the Fourier's Law: , in which k_j is the thermal conductivity in the x_j direction, that represents the spatial independent variable. In an axisymmetric model using cylindrical coordinates, x_j represents the coordinates r and z.

The problem described by the Eq. (1) is subjected to boundary and initial conditions. The boundary conditions of interest may be of three different types:

a) Dirichlet boundary condition: prescribed temperature over a portion of the boundary G_D:

b) Neumann boundary condition: prescribed normal heat flux over G_N:

in which n_j are the outward normal direction cosines.

c) Cauchy or Robin boundary condition: mixed boundary condition, in other words, prescribed flux and convection heat transfer over G_C:

with h_S representing the convective heat transfer coefficient and T_h being the environmental fluid temperature.

The initial distribution of temperature is known for an initial time stage tⁱ, and the initial condition is expressed by:

An Axisymmetric Finite Volume Formulation

Equation (1) can be re-written in terms of the fluxes using an integral formulation over an axisymmetric volume W, as:

The infinitesimal volume in a tri-dimensional model using cylindrical coordinates is given by: , where dA is the infinitesimal cross-section area. For the axisymmetric case Dq = 2p, and the axisymmetric volume results in:

Substituting dW given by Eq. (7) in Eq. (6), excepted for the source term which will be considered in detail in the next session, we have:

The previous equation can also be written as:

Equation (9) can be re-written adopting an indicial notation, after the use of the divergence theorem, as:

in which S is the cross-sectional boundary of the volume W (See Figure 1).

To obtain the FV numerical formulation of Eq. (10), the computational domain is discretized into an unstructured assembly of elements. The first integral in Eq. (10) can be computed over the axisymmetric control volume associated to node I. The integral over the boundary present in the same equation is computed over the boundary face of the control volume associated to node I using an edge-based representation of the mesh. The last term, i.e. the source (or sink) thermal loads, will be detailed later in the next section. Following the approaches presented in Lyra et al. (2004a), the semi-discrete formulation of Eq. (10), for a node I, can be conveniently expressed as:

in which A_I is the cross-sectional area of the control volume associated to node I, r_C is the centroid radius of that control volume and represents the numerically calculated temperature at node I.

The centroid of the control volume is given by:

with r_i being the centroid of each sub-element that form the control volume and A_i is the area of the corresponding sub-element. It is important to observe that the control volume centroid coordinates do not necessarily coincide with the coordinates of node I.

Coefficients and in Eq. (11) represent the components in j direction of the area vector normal to the control volume surface. They must multiply the flux associated to edge IJ_L in order to obtain the heat flux contribution of this edge to the node I. The first summation in Eq. (11) extends over all edges in the domain which are connected to node I, while the second summation corresponds to the flux contributions over a boundary edge L connected to node I. Coefficients and , which correspond to the axisymmetric model proposed in the present work, are computed as:

in which, , with. In other words, r_k represents the mid point radial coordinate of interface K and L_K is the length of each interface K associated to edge IJ_L. Each interface connects the element centroid (C) to the mid point (MP) of one of the edges that belong to such element (for detail see Figure 1b). For each boundary edge, coefficient must be calculated with , in which L_L is the half-length of the boundary edge under consideration, and, for node I, and , for node J_L. The outward normal direction cosines components, in j direction, are given by and .

Figure 1 represents an example of a typical axisymmetric control volume and the coordinate system adopted. This figure shows the cross-section of an inner control volume and the geometrical parameters required for computing the weighting coefficients . For the two-dimensional case, a detailed description of these geometrical parameters can be found in Lyra et al. (2004a).

The approximation of the flux values , for the interior contribution of all edges and for the boundary edges contributions, are, respectively, given by the conventional FVM expressions as: .

In order to compute the edges fluxes, it is required to know the nodal fluxes values and, consequently, the values of the gradient of the nodal temperature. Using the divergence theorem and the approximation used to compute volume integrals over a control volume surrounding node I, we have,

where, and

On the other hand,

Using the same kind of approximation adopted to compute the boundary integral in Eq. (10) and considering Eqs. (14) and (15), we can express the nodal approximation of the gradient as:

Using the same approximations adopted for the calculation of edges fluxes, we have:

As stated in Lyra et al. (2004 a, b), the use of Eq. (16) to compute the gradients implies that the discretization of the diffusion term in Eq. (11) involves information from two layers of points surrounding the point I under consideration. Furthermore, if uniform structured quadrilateral (or hexahedral) meshes are adopted, the values computed at a given node are uncoupled from the values of those nodes directly connected to it. This fact may cause "checker-boarding" or "odd-even" oscillations (Lyra, 1994). When computing the diffusive term in non-uniform unstructured meshes, the adoption of an extended stencil and a weak coupling with the directly connected nodes may lead to some loss of robustness and reduction of convergence rate of the resulting scheme. In fact, it can be proved that this scheme is formally only first order in general triangular meshes (Svard and Nordstrom, 2003). In order to overcome such weaknesses, the gradients must be computed in an alternative way. Following the procedure suggested in the literature (Crumpton et al., 1997 and Sorensen, 2001) a better approach can be developed. Summing up, the procedure consists in a new evaluation of the domain edges fluxes, considering the flux contributions in the parallel and normal directions to the edge. The parallel gradient contribution is computed by a new approximation using a second order central difference (Lyra et al., 2004a), and the normal component is obtained as follows,

where each nodal gradient is calculated in a finite volume fashion (Eq. 16) and the edge gradient is given by the average between the gradients of nodes I and J_L. With the new approximation of the parallel component and the normal component given in Eq. (18) we have the gradient approximation and the corresponding fluxes , which is obtained using the Fourier's conductivity law.

With the new edge flux approximation (), Eq. (11) can be re-written as:

It is important to emphasize that, for each boundary edge, it must be calculated a different coefficient for each node J_L and I, as described just after Eq. (13).

Thermal Loads Discretization

The term Q represents the thermal sources that can act in different portions of the domain. The integral form of the thermal source Q, described in Eq. (19), can be approximated by:

where the superscripts P, C, R account for thermal sources (or "sinks") acting on a point, a curve or a region, respectively. In Equation (20), the summation is over the two edges L that approximate the curve C in the point I in which the source acts. The term considers the total value of heat source for a unity volume associated to a given node I. The term containing considers the heat source value for a strip of unity width over the plane and along the surface , which represents the transversal section along the line heat source. So , gives the heat source per unit of area over the surface for a given node I. The term represents the heat source per unit of volume.

It is important to emphasize that the expression proposed in Eq. (20) is of course only valid in the case of thermal sources acting on axisymmetric curve or region, or point sources applied over the symmetry axis.

Boundary Conditions Discretization

For the portion of the boundary subjected to Dirichlet boundary condition, the nodal value of temperature is known, being the prescribed temperature value .

For the imposition of Neumann boundary condition, the normal component of the flux (coming from FVM calculation) must be substituted by the prescribed flux .

The Robin boundary condition is implemented similarly to Neumann boundary condition, since in the explicit time formulation adopted the right hand side of Eq. (4) is known when we take the temperature value T of the previous time step. Further details are presented in Lyra et al. (2004a).

Multi-Materials Domain

When heat transfer problems involve different material properties, the discretization of the governing equations should guarantee the correct solution through the sub-regions interfaces. The mesh generator used (Lyra & Carvalho, 2000) has the flexibility to generate consistent meshes over multi-region domains.

For each edge at the interface of two regions, the edge coefficient is computed independently for each region with the identification of the region also kept in memory. Therefore, each interface edge has two weighting coefficients, defined in Eq. (13) (Lyra et al., 2004a).

In case of multi-materials, the discrete Equation (19) including the thermal loads terms, is then replaced by:

It is worthy to emphasize that the third summation of the right side is different of zero only when the node I is on an interface of two or more regions with different properties. The interface edge fluxes are computed in the same fashion as .

The gradient values and respective fluxes are obtained in three stages. At first, the summation over the boundary edges with the Neumman and Robin boundary conditions, already prescribed (second term of right hand side – RHS – of Eq. (21)). The second stage involves the summation over all domain edges (inner and boundary) (first term of RHS Eq. (21)). Finally a double summation over the interface edges is performed (third term of RHS Eq. (21)). During each summation the coefficients and material properties used are those correspondents to each region.

As the 2p constant appears in all equation terms, explicitly or implicitly as in the case of coefficients and , in practice it is removed from the computational implementation.

This procedure represents a good approximation for low ration of different material properties. For high ration, a more robust approximation is under development, see Carvalho (2005) for futher details.

Time Discretization

The time discretization is done through a simple explicit formulation ("Euler forward"), where the temperatures of node I are computed as functions of the neighborhood temperatures calculated in the previous time step (Maliska, 1995):

where RHS refers to the right hand side of Eq. (21), is the known temperature at time and the unknown temperature at time .

Numerical Results

In this section some simple, though representative academic examples are presented in order to validate and show some of the capabilities of the numerical scheme previously proposed. All steady-state solutions were obtained using the explicit transient algorithm, which is stopped when the time derivative of the temperature is not changing, i.e. the residual is below a pre-assigned tolerance, with the arbitrarily adopted value of 10^-5.

1^st Example: Steady-State Heat Transfer in a Solid Cylinder

The first example consists in the calculation of the steady state temperature profile in a stainless steel (AIST 302) solid cylinder of dimensions 0 < r < b e 0 < z < a. The heat is generated inside the cylinder at a constant rate of Q = 2.4 W/m³. The surfaces at z = 0 and z = a are insulated. The surface at r = b is under convection to an ambient temperature of T_¥ = 4.0°C, and the convective heat transfer coefficient is h = 10.0 W/m²°C. The thermal conductivity of the cylinder is k = 15.1 W/m°C, its mass density is r = 8055.0 kg/m³ and its specific heat is c = 480.0 J/kg K (Incropera & DeWitt, 1998). The cylinder dimensions are b = 5.0 m and a = 5.0 m. Figure 2 shows the boundary conditions of the problem. The trivial one-dimensional analytical solution that satisfies the boundary conditions of this problem can be found in Özisik (1980). For the numerical solution of the problem, three uniform isotropic unstructured triangular meshes were used, the first one with 52 elements and 37 nodes, the second one, with 102 elements and 66 nodes; and the third one with 1368 elements and 735 nodes. In Fig. 3, the radial temperature distributions for z = 0.0 using the three meshes, are plotted together with the analytical solution and Figure 4 shows the temperature contours for the different meshes. The maximum relative errors are 0.45%, 0.33% and 0.06% for the 1^st, 2^nd and 3^rd meshes, respectively. Even for the very coarse mesh we get good results and this simple mesh study was presented just to show the proper convergence of the solution.

2nd Example: Transient Heat Transfer in a Solid Cylinder

The second example presents the transient temperature profile in a cylinder of 5.0 m high and 5.0 m radius without heat generation and initially at T_o = 30.0°C. The physical properties are identical to the previous example. The surfaces z = 0.0 and z = 5.0 m are kept insulated, and the surface in r = 5.0 m is kept at =0. Figure 5 presents the boundary conditions of the problem. The analytical solution that satisfies the boundary conditions of this problem can be found in Özisik (1980). An uniform isotropic unstructured triangular mesh with 37 nodes and 52 elements was used. Figure 6 shows the temperature distribution versus time of a node at r = 3.0 m and z = 0.0. Figure (7) shows the temperature distributions for z = 0.0 at t = 110 h, t = 345 h and t = 555 h. Comparing analytical and numerical solutions, the maximum relative error for this coarse mesh was 4.9%.

3^rd Example: Steady State Heat Conduction in Nuclear Reactor Element

A solid cylindrical fuel element with 6 mm radius is made of UO₂ (uranium oxide). A 3-mm-thick cladding surrounds it. The cylinder is 6 mm high and the surfaces at z = 0.0 and z = 6 mm are kept insulated. The volumetric heat generation rate in the fuel element is q = 2.10⁸ W/m³. The coolant temperature is T_¥ = 27°C and the heat transfer coefficient is h = 2000 W/m²°C. The thermal conductivity of UO₂ is k_f = 2 W/m°C, its mass density is r_f = 10500 kg/m³ and its specific heat is c_f = 1033 J/kg K. The thermal conductivity of the cladding is k_c = 25 W/m°C, its mass density is r_c = 6504 kg/m³ and its specific heat is c_c = 1422 J/kg K. Figure 8 shows a cross-section of the fuel element with the cladding and the boundary conditions associated to the problem. This example is taken from Incropera & DeWitt (1998). The analytical solution that satisfies the boundary conditions of this problem can be found in El-Wakil (1971). We have used an uniform isotropic unstructured triangular mesh with 273 nodes and 484 elements. The radial temperature distribution is presented in Fig. 9. We can observe that the numerical solution presents the expected profile. For 0 < r < 6mm the profile is hyperbolic due to the heat generation in the fuel element, and the temperature in the cladding presents a linear profile. Table 1 shows the temperatures at the center of the fuel element, at the interface fuel-cladding, and at the surface of the cladding. We can see the good agreement achieved when we compare the numerical solution with the analytical one. The maximum relative error obtained was 0.08%. Such example validates the multi-materials problem approximation, while it extends the applicability of the presented formulation.

4^th Example: Steady-State Heat Conduction in a Radial Fin of Hyperbolic Profile

A stainless steel (AIST 302) radial fin is fixed on a circular vessel. The base width (z_b) of the fin is 4.0 mm and the inner (r_b) and outer (r_o) radii are 25.0 mm and 40.0 mm, respectively. The hyperbolic profile is given by. The surrounding fluid is at T_¥ = 20.0°C and the coefficient of heat transfer is h = 40.0 W/m²°C. The temperature of the base of the fin is T_b = 200.0°C. A negligible amount of heat is transferred out of the end of the fin. The thermal conductivity of stainless steel k =15.1 W/m°C, its mass density is r = 8055.0 kg/m³ and its specific heat is c = 480.0 J/kg K. Figure 10 presents a sketch of a radial fin with hyperbolic profile. Due to the symmetric nature of the problem around the central axis (z = 0), a symmetric model was adopted where it is necessary to analyze only a half of the domain. After a study of convergence, an uniform isotropic unstructured triangular mesh with 134 nodes and 199 elements was chosen.

The analytical solution of an one-dimensional model that approximates this problem can be obtained from the generalized differential equation (Kraus et al, 2000),

with q=T-T_¥and the following boundary conditions:

The resulting solution is:

where M=(2h/kz_br_b)^1/2 and I_n(x) are the modified Bessel function of the first kind and order "n".

Figure 11 presents the steady state temperature contour and the mesh utilized. As expected, the temperature profile varies along the z-axis, because the actual problem is axisymmetric. Despite this fact, excellent concordance of the axisymmetric numerical solution and the analytical solution of the approximate 1-D model problem is obtained along the central axis (z = 0) of the radial fin. Figure 12 shows the numerical solution of the temperature distribution. The maximum relative error obtained when the numerical solution is compared with the 1-D analytical one (Eq. (24)) is 0.091%. It must be said that the most correct solution (i.e. closer to the physical solution) in this case is the axisymmetric numerical solution. For a different configuration of the radial fin, such as when it becomes smaller (increasing Biot number), the 1-D analytical solution does not represent anymore a good approximation to the physical problem. These conclusions can be confirmed by the results shown in Fig. 13. The maximum relative error in this case is 1.5%.

Conclusions

An axisymmetric finite volume formulation, using unstructured meshes is presented in this paper. This formulation was used to solve some simple, though representative heat transfer problems. The axisymmetric finite volume formulation allows the treatment of certain three-dimensional problems at low cost and high accuracy. This formulation also allows for the utilization of several two-dimensional mesh adaptation tools developed in our research group with little effort (Araújo, et al, 2004) (work under development). Furthermore, the formulation can be extended to deal other classes of problems governed by conservation laws following with similar procedures as presented here.

In the near future, this formulation will be used to solve bioheat transfer applications such as the study of the temperature distribution in human eye with retinal implants (subretinal or epiretinal). These two kinds of implants are being used in humans who have retinitis pigmentosa or age-related macular degeneration, which are some of the leading causes of blindness in the population (Margalit et al., 2002). Lima et al. (2004) have performed a two-dimensional study of the temperature distribution in human eyes, both with and without retinal implants.

Acknowledgements

The authors would like to acknowledge the Brazilian Research Council (CNPq), the National Petroleum Agency (ANP) and (CAPES) for the financial support provided during the development of this research.

Paper accepted July, 2005.

Technical Editor: Atila P. Silva Freire.

• AGARD Report 787, 1992, "Special Course on Unstructured Grid Methods for Advection Dominated Flows", France.
• Araújo F. S., Lyra, P. R. M., Carvalho, D. K. E., 2004, "Procedimentos Adaptativos Aplicados À Simulação De Reservatórios De Petróleo", Iberian Latin American Congress In Computational Methods In Engineering, CILAMCE 2004, Recife - PE, Brasil, In CD-ROM (in Portuguese).
• Baliga, B.R., Patankar, S.V., 1983, "A Control Volume Finite-Element Method for Two-Dimensional Fluid Flow and Heat Transfer", Numerical Heat Transfer, Vol. 6, pp. 245-261.
• Barth, T. J., 1992, "Aspects of Unstructured Grids and Finite-Volume Solvers for The Euler and Navier-Stokes Equations", AGARD Report 787, pp. 6.1-6.61.
• Crumpton, P.I., Moinier, P. & Giles, M.B.T.J., 1997, "An Unstructured Algorithm for High Reynolds Number Flows on Highly Stretched Grids". In: Taylor, C. & Cross, J. T, ed., Numerical Methods In Laminar and Turbulent Flow, Pineridge Press, pp. 561-572.
• El-Wakil, M.M., 1971, "Nuclear Heat Transport", Ed. International Textbook Company, New York.
• Incropera, F.P., DeWitt, D.P., 1998, "Fundamentos da transferência de calor e massa", 4Ş ed., LTC, Rio de Janeiro (in Portuguese).
• Kraus, A.D., Aziz, A., Welty, J., 2000, "Extended Surface Heat Transfer", John Wiley & Sons Inc.
• Lima, R de C.F. de, Lyra, P.R.M., Silva, G.M.L.L. da & Carvalho, D.K.E. de, 2004, "Análise térmica dos tecidos oculares dotados de implantes retinianos através da utilização do método dos volumes finitos em malhas não estruturadas", Iberian Latin American Congress on Computational Methods in Engineering, CILAMCE 2004, Recife - PE, Brasil, in CD-ROM (in Portuguese).
• Lyra, P.R.M., 1994, "Unstructured Grid Adaptive Algorithms for Fluid Dynamics and Heat Conduction", Ph.D. thesis C/PH/182/94, University of Wales Swansea.
• Lyra, P.R.M., Carvalho, D.K.E. de, 2000, "A Flexible unstructured mesh generator for transient anisotropic remeshing", In: ECCOMAS 2000 European Congress on Comp. Meth. In Applied Sciences and Eng., Barcelona, Espanha, in CD-ROM.
• Lyra, P.R.M., Lima, R de C.F. de, Guimarães, C.S.C., Carvalho, D.K.E. de, 2002, "Uma formulação com estrutura de dados por arestas do método dos volumes finitos na solução de problemas de potencial", Anais do MECOM'2002 - First South American Congress on Computacional Mechanics, Parana - Santa Fé, Argentina, in CD-ROM (in Portuguese).
• Lyra, P.R.M., Lima, R de C.F. de, Guimarães, C.S.C., Carvalho, D.K.E. de, 2004a," An edge-based unstructured finite volume procedure for the numerical analysis of heat conduction applications", J. of the Braz. Soc. of Mech. Sci. & Eng., Vol. 26(2), pp. 160-169.
• Lyra, P.R.M., Lima, R de C.F. de, Silva, G.M.L.L. da, Carvalho, D.K.E. de, 2004b, "Uma formulação axissimétrica do MVF em malhas não-estruturadas para solução de problemas transientes de transferência de calor", Anais do CONEM2004 Congresso Nacional de Engenharia Mecânica, Belém - PA, Brasil, in CD-ROM (in Portuguese).
• Maliska, C.R., 1995, Transferência de Calor e Mecânica dos Fluidos Computacional; Fundamentos e Coordenadas Generalizadas, Ed. LTC, Rio de Janeiro (in Portuguese).
• Margalit, E., Maia, M., Weiland, J.D., Greenberg, R.J., Fujii, G.Y., Torres, G., Piyathaisere, D.V., O´Hearn, T.M., Liu, W., Lazzi, G., Dagnelie, G., Scribner, D.A., de Juan Jr, E., Humayun, M.S., 2002, "Retinal prosthesis for the blind", Survey of ophthalmology, vol. 47 (4), pp. 335-356.
• Mokheimer, E.M.A., 2002, "Performance of annular fins with different profiles subject to variable heat transfer coefficient", International Journal of Heat and Mass Transfer, vol. 45 (2002), pp.3631-3642.
• Özisik, M.N., 1980, "Heat Conduction", New York, Ed. John Wiley & Sons.
• Peraire, J., Peiró, J. & Morgan, K., 1993, "Finite Element Multigrid Solution of Euler Flows Past Installed Aero-Engines", J. Computational Mechanics, Vol.11, pp. 433-451.
• Sorensen, K.A., 2001, "A Multigrid Procedure for The Solution of Compressible Fluid Flows on Unstructured Hybrid Meshes", Ph.D. thesis C/PH/251/01, University of Wales Swansea.
• Svard, M. and Nordstrom, J., 2003, "A stable and accurate summation-by-parts finite volume formulation of the laplacian operator", Technical Report, 2003-03, Uppsala Universitet.
• Zienkiewicz, O.C. and Morgan, K., 1983, "Finite Element and Approximation", John Wiley & Sons, Inc.

Publication Dates

• Publication in this collection
  02 Jan 2006
• Date of issue
  Dec 2005

History

• Received
  July 2005
• Accepted
  July 2005