Yi Zhong and Marek-Jerzy Pindera

University of Virginia, Charlottesville, Virginia

Efficient Reformulation of HOTFGM:

Heat Conduction With Variable

Thermal Conductivity

NASA/CR—2002-211910

November 2002
The NASA STI Program Office . . . in Profile

Since its founding, NASA has been dedicated to

the advancement of aeronautics and space

science. The NASA Scientific and Technical

Information (STI) Program Office plays a key part

in helping NASA maintain this important role.

The NASA STI Program Office is operated by

Langley Research Center, the Lead Center for

NASA’s scientific and technical information. The

NASA STI Program Office provides access to the

NASA STI Database, the largest collection of

aeronautical and space science STI in the world.

The Program Office is also NASA’s institutional

mechanism for disseminating the results of its

research and development activities. These results

are published by NASA in the NASA STI Report

Series, which includes the following report types:

• TECHNICAL PUBLICATION. Reports of

completed research or a major significant

phase of research that present the results of

NASA programs and include extensive data

or theoretical analysis. Includes compilations

of significant scientific and technical data and

information deemed to be of continuing

reference value. NASA’s counterpart of peerreviewed formal professional papers but

has less stringent limitations on manuscript

length and extent of graphic presentations.

• TECHNICAL MEMORANDUM. Scientific

and technical findings that are preliminary or

of specialized interest, e.g., quick release

reports, working papers, and bibliographies

that contain minimal annotation. Does not

contain extensive analysis.

• CONTRACTOR REPORT. Scientific and

technical findings by NASA-sponsored

contractors and grantees.

• CONFERENCE PUBLICATION. Collected

papers from scientific and technical

conferences, symposia, seminars, or other

meetings sponsored or cosponsored by

NASA.

• SPECIAL PUBLICATION. Scientific,

technical, or historical information from

NASA programs, projects, and missions,

often concerned with subjects having

substantial public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific

and technical material pertinent to NASA’s

mission.

Specialized services that complement the STI

Program Office’s diverse offerings include

creating custom thesauri, building customized

databases, organizing and publishing research

results . . . even providing videos.

For more information about the NASA STI

Program Office, see the following:

• Access the NASA STI Program Home Page

at http://www.sti.nasa.gov

• E-mail your question via the Internet to

help@sti.nasa.gov

• Fax your question to the NASA Access

Help Desk at 301–621–0134

• Telephone the NASA Access Help Desk at

301–621–0390

• Write to:

 NASA Access Help Desk

 NASA Center for AeroSpace Information

 7121 Standard Drive

 Hanover, MD 21076
Yi Zhong and Marek-Jerzy Pindera

University of Virginia, Charlottesville, Virginia

Efficient Reformulation of HOTFGM:

Heat Conduction With Variable

Thermal Conductivity

NASA/CR—2002-211910

November 2002

National Aeronautics and

Space Administration

Glenn Research Center

Prepared under Grant NAG3–2524
Acknowledgments

This research was conducted under funding through the NASA Glenn Grant NAG3–2524. The NASA technical

monitor for this project was Dr. Steven M. Arnold whose support and encouragement for the undertaken

investigation, the authors gratefully acknowledge.

Available from

NASA Center for Aerospace Information

7121 Standard Drive

Hanover, MD 21076

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22100

Available electronically at http://gltrs.grc.nasa.gov
Abstract

Functionally graded materials (FGMs) have become one of the major research topics in the

mechanics of materials community during the past fifteen years. FGMs are heterogeneous materials, characterized by spatially variable microstructures, and thus spatially variable macroscopic properties, introduced to enhance material or structural performance. The spatially

variable material properties make FGMs challenging to analyze. The review of the various

techniques employed to analyze the thermomechanical response of FGMs reveals two distinct

and fundamentally different computational strategies, called uncoupled macromechanical and

coupled micro-macromechanical approaches by some investigators. The uncoupled macromechanical approaches ignore the effect of microstructural gradation by employing specific spatial

variations of material properties, which are either assumed or obtained by local homogenization,

thereby resulting in erroneous results under certain circumstances. In contrast, the coupled approaches explicitly account for the micro-macrostructural interaction, albeit at a significantly

higher computational cost. The higher-order theory for functionally graded materials (HOTFGM) developed by Aboudi et al. (1999) is representative of the coupled approach. However,

despite its demonstrated utility in applications where micro-macrostructural coupling effects are

important, the theoryís full potential is yet to be realized because the original formulation of

HOTFGM is computationally intensive. This, in turn, limits the size of problems that can be

solved due to the large number of equations required to mimic realistic material microstructures. Therefore, a basis for an efficient reformulation of HOTFGM, referred to as user-friendly

formulation, is developed herein, and subsequently employed in the construction of the efficient

reformulation using the local/global conductivity matrix approach. In order to extend HOTFGMís range of applicability, spatially variable thermal conductivity capability at the local level

is incorporated into the efficient reformulation. Analytical solutions to validate both the userfriendly and efficient reformulations are also developed. Volume discretization sensitivity and

validation studies, as well as a practical application of the developed efficient reformulation, are

subsequently carried out. The presented results illustrate the accuracy and implementability of

both the user-friendly formulation and the efficient reformulation of HOTFGM .

1 Introduction

Functionally graded materials (FGMs) have become one of the major research topics in the mechanics and materials communities during the past fifteen years. FGMs are heterogeneous materials

which are composed of two or more phases with different properties, and therefore they belong to

Efficient Reformulation of HOTFGM: Heat Conduction

 With Variable Thermal Conductivity

Yi Zhong and Marek-Jerzy Pindera

University of Virginia

Charlottesville, Virginia 22903

NASA/CR—2002-211910 1
the broad class of materials known as composites. Unlike traditional composites, however, the spatial distribution of phases in FGMs is not statistically homogeneous. As the name suggests, FGMs

are materials with macroscopic properties characterized by spatial gradients. This is accomplished

by using variable spacings between individual inclusions or by using inclusions with different properties, sizes and shapes. Research results obtained thus far have demonstrated that FGMs have

great potential for improving material/structural performance in many engineering applications

precisely because of their spatially graded heterogeneous microstructure. In particular, grading or

tailoring the internal microstructure of a composite material or structural component enables the

designer to integrate both the material and structural aspects into the final design and product.

This offers a number of advantages over traditional methods of tailoring composite materials and

structural components through the elimination of sharp material discontinuities which often lead

to undesirable stress concentrations, thereby degrading the structural componentís strength and

durability.

The concept of modern man-made FGMs was first proposed in 1984 by material scientists in

the Sendai area of Japan as a means of developing thermal protection materials. In particular, this

concept was originally conceived for thermal protection of fuselage exterior and engine materials

intended for use in high-speed civil transport aircraft under development in the mid 1980ís through

mid 1990ís, and then rapidly spread to other fields including physics, chemistry, agriculture, biology

and medicine. After almost twenty years of development, the term FGM is now widely used in

various fields to describe this class of materials. Strictly speaking, however, graded materials are not

new. They can be found in nature, and they also have been utilized extensively by human beings

for many years without being explicitly called FGMs. The examples include graded materials

developed long ago, such as case-hardened steel, which are still in common use today, Kaysser and

IIschner (1995). Other naturally occurring biological FGMs include bamboo, which has been used

for a long time for structural components, corn, barley, etc., Figure 1. Nevertheless, the new and

exciting thing about modern FGMs is the realization that the property gradient can be designed

at the microstructural level to tailor the materialís specific function and performance required by

engineering needs.

Figure 1. Selected FGMs found in nature. (Source: Amada, 1995).

As stated in the foregoing, FGMs are excellent candidates for applications involving severe

thermal gradient environments. Some practical examples of the successful application of the FGM

NASA/CR—2002-211910 2
concept include zirconia coatings deposited on a nickel-chromium alloy substrate by plasma spraying, where the composition profiles of the graded interlayers are carefully optimized to minimize

thermally-induced stresses under service loading conditions. By appropriately adjusting the microstructural transition, that is by placing the ceramic-rich part of a functionally graded material

in the high temperature environment and the metallic-rich part in the low temperature environment, with a gradual microstructural transition in the direction of the temperature gradient, it

is possible to achieve optimum temperature, deformation and stress distributions. An example

of a graded thermal barrier coating is given in Figure 2. In this case, the requirements for heatresistant and oxidation-resistant properties are fulfilled, as well as those for adequate mechanical

toughness. In another successful application, Honda Engineering Co. Ltd. developed a functionally

graded cemented carbide with improved wear resistance and tool life, owing to a coarse grain size

and smaller cobalt content in the surface region than in the core region with a graded interlayer,

Ichikawa (2001). Thus the highest hardness was achieved in the surface layer while the toughness

increased towards the core region.

Figure 2. Graded ceramic-metal thermal barrier coating (Source: Dollmeier et al. 1995).

There is no doubt that FGMs offer substantial advantages over traditional composites in the

design of structural components which require different types of materials to simultaneously satisfy

various requirements such as load bearing capability, surface protection, and perhaps decorative

function. In particular, it has been demonstrated that microstructural gradients not only can overcome problems associated with traditional composites by smoothing out property discontinuities,

but can also produce unique and desirable functions, such as focusing light in fibers, channeling heat in computer chips, and providing biocompatible implant-tissue transitions in biomedical

engineering applications, for instance.

However, there are many challenges that must be overcome in order to make full use of the FGM

potential, Ilschner (1997). Materials researchers agree that the FGM manufacturing capabilities

and material development are not in a dynamic growth phase comparable to that which occurred

during the 1960-1980 period for the then-emerging materials such as plastics. The reason is that

with the present state of knowledge, it is very difficult and expensive to fabricate graded materials

with various designed functions in large quantities for industrial use. Moreover, at the present

time, mature design standards and test methods for FGMs are insufficient or unavailable. Just as

importantly, the multiscale approach required in the analysis, design and optimization of FGMs

is still not well developed. The present investigation, therefore, was undertaken to address this

last point, as will be discussed in more detail later. The literature review of the following section

provides the motivation for the undertaken study.

NASA/CR—2002-211910 3
1.1 Literature Review

In this section, a brief overview of the FGM development is given for completeness, followed by

a discussion of theoretical modeling methods used to analyze FGMs. The emphasis is put on the

latter in view of the focus of this investigation which addresses one aspect of multiscale modeling

of FGMs, as discussed in the next section.

With the realization of the important role that FGMs will play in many fields, many researchers

have devoted a substantial amount of time and energy in recent years to the broad area of FGMs,

ranging from the different fabrication techniques and characterization methods to the development

of theoretical models for the response of these materials. The development and exploration of

FGMs was initiated in Japan under a national project during the period 1987-1991. The aim of

this government-sponsored project was to develop FGMs for future aerospace applications. Following this, similar projects and forums for FGM research and development were initiated in other

countries, in particular in Germany, Finland, China, the United States and Russia. Research on

functionally graded materials around the world increased rapidly after the First International Symposium devoted solely to these materials, which was held in Sendai, Japan in 1990. The second

International Symposium on FGMs was held in San Francisco in 1992, followed by the Third International Symposium held in Lausanne, Switzerland in 1994. These symposia, which were organized

by the International Forum on FGMs, are now held every two years on a regular basis with the

next one scheduled for the fall of 2002 in Beijing, China.

As a result of the above activities, the overall research and development of FGMs has expanded

to many areas. Consequently, many papers on different aspects of FGMs have been published in

standard journals, special issues of journals devoted solely to FGMs, conference proceedings, and

monographs, cf. Pindera et al. (1994a, 1995a, 1997), Needleman and Suresh (1996), Yamanouchi et

al. (1990), Ilschner and Cherradi (1994), Shiota and Miyamoto (1997), and Suresh and Mortensen

(1998). Therefore, a comprehensive review of the different research activities is outside the scope

of this report. Therefore, in keeping with the investigationís focus, an overview is provided below

of the different approaches employed by researchers to model the thermomechanical response of

FGMs.

The review of the various techniques employed to analyze the response of FGMs under thermal

and mechanical loading reveals two distinct and fundamentally different computational strategies.

These have been termed uncoupled and coupled micro-macromechanical approaches by Aboudi

et al. (1993). The uncoupled approaches do not explicitly couple the materialís heterogeneous

microstructure with the structural global analysis. Local effective or macroscopic properties at a

given point within the FGM are first obtained through homogenization based on a chosen micromechanics scheme, and subsequently used in a global thermomechanical analysis, Figure 3. The

response of an FGM under specified loading at a homogenized material point with equivalent properties is obtained by solving the governing differential equations (partial or ordinary) with variable

coefficients which represent the spatially macroscopic material property variation obtained from

homogenization. Sometimes, the spatial property variation is represented as piece-wise uniform

which is accomplished by discretizing the microstructure into layers or strips with uniform properties that change gradually from layer to layer. Clearly, by employing only effective or homogenized

properties in the global FGM analysis, the microstructural details can only affect the FGM in an

effective or average way. It is in this sense that the approach is uncoupled.

The key assumption made in the uncoupled approach is that the heterogeneous microstructure

of an FGM can be replaced by an equivalent continuum with a set of macroscopic properties that

vary with spatial coordinates in a manner commensurate with the materialís heterogeneity. The

various micromechanical approaches that may be used to calculate homogenized material properties

NASA/CR—2002-211910 4
Figure 3. Homogenization-based analysis of FGMs.

include Reuss and Voigt hypotheses, self-consistent schemes, Mori-Tanaka method, concentric cylinder models, and the Generalized Method of Cells, as discussed by Aboudi (1991). Therefore, the

use of the uncoupled micro-macromechanical approach in analyzing the response of FGMs is implicitly based on the assumption of the applicability of the representative volume element (RVE)

concept. However, the definition of an RVE remains to be established theoretically for materials

with spatially variable heterogeneous microstructures. In addition, ignoring the possibility of coupling between local and global effects itself may result in potentially erroneous predictions. This

is especially relevant when the thermal and mechanical field variable gradients are relatively large

compared with the characteristic dimension of the inclusion phase or the number of inclusions is

small, or the dimension of the inclusion phase is large relative to the global dimension of the composite, Aboudi et al. (1993, 1994). Therefore, this type of approach should be used in limited

circumstances such as when the size of the inclusion is relatively small compared with the overall

size of composite or the number of the inclusions is large, Pindera et al. (1995b).

The alternative to the uncoupled approach is the coupled approach wherein the micro-macrostructural

interaction is explicitly taken into account. This implies that the local micromechanics analysis

must be combined with the global or structural analysis in a way that connects the two scales

directly. One way to do this is to solve the entire problem directly by accounting for the presence of

every local heterogeneity. In the presence of many inclusions, the solution of the given boundaryvalue problem in an exact sense is prohibitive and often impossible, necessitating development of

approximate computational strategies.

NASA/CR—2002-211910 5
1.1.1 The uncoupled macrostructural approach

The uncoupled approach is the simpler of the two approaches in the analytical sense since it

treats the micro and macrostructural problems separately. As a result, many particular types

of problems involving FGM applications have been solved using this approach. One of the most

important applications involves thermal barrier coating (TBC) problems. As discussed in the

preceding section, internal stresses can be reduced and fracture toughness can be enhanced by

appropriate choice of microstructural gradation. For instance, Fukumoto et al. (1994) addressed

the influence of coating configuration on cyclic thermal shock fracture behavior of plasma sprayed

coatings. An extensive review of the different approaches employed to evaluate temperature and

stress fields in TBCs has been recently provided by Noda (1999). In the presence of inelastic

effects, the approaches range from one-dimensional finite-difference analyses, Kokini and Choules

(1995), to two-dimensional finite-element analyses based on layer-wise discretization of the coating

microstructure with piece-wise uniform properties calculated using simple rule-of-mixtures or MoriTanaka estimates, such as those by Jian et al. (1995) or Delfosse et al. (1997). Elastic analyses

often use continuously varying properties of particular functional form, thereby facilitating either

closed-form analytical solutions, or numerical solutions based on the Greenís function approach,

such as those presented by Nomura and Sheahen (1997) and Sutradhar et al. (2002).

Another important application of FGMs involves joining structural components made of dissimilar materials using graded transition regions to reduce the large interlaminar stresses that arise

at the free edge due to the material property mismatch. Analyses of these problems have been

conducted by Drake et al. (1993) and Williamson et al. (1993) using the finite-element approach

to study the effect of material gradation on the free-edge peel stresses in layered cylinders with

gradual material variation achieved by using slightly different material properties in each cylindrical layer. Plasticity was demonstrated to be an important factor in influencing the free-edge stress

fields. Suresh et al. (1994) demonstrated similar effects in analyzing the response of elastoplastic

biomaterial strips subjected to cyclic thermal loading.

Fracture analysis is also an important area of FGMs. In a sequence of papers, Erdogen and coworkers extended the elasticity techniques for crack-analysis problems of inhomogeneous materials

to functionally graded materials, Delale and Erdogan (1983), Erdogan (1985), Erdogan et al. (1991).

The focus of these studies were the different failure modes that occur in graded TBCs which are

initiated by free-edge cracks, and cracks parallel or perpendicular to the coating surface, as discussed

by Kokini et al. (1996) or Kawasaki and Watanabe (1997). These were subsequently analyzed by

Erdogan (1995), Erdogan and Wu (1996), Lee and Erdogan (1998), and most recently by Jin and

Paulino (2001).

Other applications of FGMs include the important problem of identifying optimum microstructural designs which produce the best temperature and stress distributions under given loading

conditions. Various analytical solutions have been developed for specific boundary-value problems

(typically involving axisymmetric fields and out-of-plane shearing) which can be incorporated into

optimization algorithms for design purposes. These analytical solutions are also useful as benchmark solutions for the validation of finite-element, finite-difference or boundary-element techniques.

For instance, Horgan and Chan (1999) developed elastic solutions for problems involving pressurized hollow cylinders, rotating disks, and bars under torsion with continuously graded isotropic

constituents. In order to enable fiber-matrix interface/interphase design and optimization in metal

matrix composites, Pindera et al. (1994b) extended the previously-developed thermoplastic solution to the problem of an arbitrary layered concentric cylinder under combined thermomechanical

axisymmetric loading, Pindera et al. (1993), by incorporating the capability to treat the interfacial

layers as two-phase composite materials. The results provided insight into designing/engineering

NASA/CR—2002-211910 6
homogenized and functionally graded interfaces in advanced metal-matrix composites with the objective of generating preferred fabrication-induced residual stress distributions. Salzar and Barton

(1994) subsequently incorporated this solution strategy into a commercial optimization code and

demonstrated the efficiency with which optimum interfacial architectures could be identified.

1.1.2 The coupled micro-macrostructural approach

The coupled micro-macromechanical approach circumvents the limitations of the standard uncoupled approach in FGM applications by explicitly coupling the local (microstructural) and global

(macrostructural) effects. The coupled analysis can be carried out using the finite-element technique. However, this can be prohibitively expensive due to the refined volume discretization required

for convergence when dealing with even simple FGM microstructures. For this reason, alternative

approximate computational techniques are needed which are sufficiently accurate to compute the

micro-macrostructural interaction without requiring excessively large storge.

The representative of such a coupled micro-macromechanical approach is the higher-order theory

for functionally graded materials (HOTFGM) developed by Aboudi and co-workers in a sequence of

papers and reports originating in 1993, Aboudi et al. (1993). Summaries of the different stages of

the higher-order theory have been provided by Pindera et al. (1995c, 1998), and most recently by

Aboudi et al. (1999). The theoretical framework is based on volumetric averaging of the various field

quantities, satisfaction of the field equations in a volumetric sense, and imposition in an average

sense of boundary and interfacial conditions between the subvolumes used to characterize the

compositeís functionally graded microstructure. The theory has been validated through comparison

with the results obtained from finite-element analysis, Pindera and Dunn (1997), and subsequently

applied to the following technologically important problems:

2 Investigation of the effect of microstructure on thermal and stress fields in MMC plates and

cylinders

2 Investigation of the use of functionally graded architectures in reducing edge effects in MMC

plates

2 Optimization of functionally graded microstructures in MMC plates and cylinders

2 Development of guidelines for the design of special coatings in exhaust nozzle applications

2 Investigation of the microstructural effects in functionally graded TBCs

The results obtained so far have demonstrated that this theory is an accurate and easily implementable tool in the analysis and design of FGMs. Furthermore, comparison of the results obtained

from the standard micromechanics approach with those of HOTFGM has demonstrated the need

for a theory like HOTFGM, which explicitly couples the micro (local) and macro (global) effects

within a unified analysis.

The recent developments of HOTFGM include extension to cylindrical coordinates to enable

analysis, design, and optimization of structural components found in aircraft engine applications,

Pindera and Aboudi (2000).

1.2 Objectives of the Investigation and Outline of the Report

Despite the demonstrated utility of the higher-order theory in applications where the micromacrostructural coupling effects are important, the theoryís full potential has not yet been attained.

NASA/CR—2002-211910 7
This is because the theory as originally formulated is computationally intensive. This, in turn, limits the types of problems that can be solved due to the large number of equations required to mimic

realistic material microstructures. Therefore, the objective of this investigation was to construct

a basis for an efficient reformulation of the higher-order theory, and subsequently implement it in

order to reduce the number of equations required in the solution of a given thermal boundary-value

problem. The basis for the efficient reformulation is the demonstration that the same results are

obtained when the volume discretization used in the original version is simplified by eliminating

the need to differentiate between generic cells and subcells. This is called a user-friendly formulation which, in turn, provides the basis for an efficient reformulation based on the local/global

conductivity matrix approach. The efficient reformulation admits spatially variable thermal conductivity within individual subcells, thereby making it possible to treat microstructures with two

fundamentally different length scales.

The theoretical framework of the original two-dimensional formulation of HOTFGM, limited to

thermal boundary-value problems, is briefly described in the first part of Section 2, with the second part devoted to the development of the user-friendly formulation. Based on the demonstrated

feasibility of the simplified user-friendly formulation, an efficient reformulation of HOTFGM with

variable thermal conductivity is developed in Section 3 using the local/global conductivity matrix

approach, which is an extension of the local/global stiffness matrix approach described by Pindera

(1991) in the context of purely mechanical boundary-value problems. Analytical solutions used

to validate both the user-friendly and efficient reformulations are developed in Section 4. Volume

discretization sensitivity and validation studies, as well as a practical application of the developed

efficient reformulation, are provided in Section 5. Section 6 summarizes the present accomplishments and discusses future desirable extensions of the developed computational capability.

2 Higher-Order Theory: Thermal Problem

2.1 Outline of the Original Formulation

The original formulation of the two-directional higher-order theory for functionally graded materials (HOTFGMó2D) was outlined by Aboudi et al. (1996). The basic idea was to develop a

two-dimensional framework for the response of composites which are functionally graded in two

directions, accounting for micro-macrostructural coupling effects. As indicated in Figure 4, the

heterogeneous composite model has finite length L and width H in the x2 ° x3 plane and extends

to infinity in the x1 direction. The thermal loading applied on the boundary can be specified in

terms of temperature or heat flux distribution. The microstructure of the heterogeneous composite

is discretized into Nq and Nr generic cells (q; r) in the intervals 0 ! x2 ! H; 0 ! x3 ! L in the

x2 °x3 plane. The indices q and r represent the location of the generic cell in the x2 °x3 plane and

they can vary from q = 1; 2::: to Nq and r = 1; 2; ::: to Nr. Each generic cell is further discretized

into four subcells denoted by the indices (Ø°), Figure 4, where the range of each index Ø; ° is from

1 to 2. These indices indicate the relative position of a given subcell within the generic cell. The

dimension of the generic cell along the x1 direction is infinite whereas the dimensions along the x2

and x3 axes, h(q)

1 ; h(q)2 and l

(r)

1 ; l(r)2 ; respectively, can vary from cell to cell such that

H = XNq

q=1

(h(q)

1 + h(q)2 ) , L = XNr

r=1

(l

(r)

1 + l

(r)

2 )

NASA/CR—2002-211910 8
Figure 4. Schematic of the HOTFGM analysis geometry.

An approximate solution for the temperature field under steady-state thermal loading is constructed on the basis of volumetric averaging of the heat conduction equation, together with the

imposition of boundary and continuity conditions in an average sense between the subvolumes used

to characterize the materialís microstructure. The solution to the heat equation is accomplished

by approximating the temperature field in each subcell of a generic cell using a quadratic expansion in the local coordinates x1

(β)

2 ; x1

(γ)

3 centered at subcellís midpoint. The higher-order terms in

the temperature field expansion account for the local effects caused by the thermal field gradients,

the microstructure of the composite and the finite dimensions in the functionally graded directions, in contrast with the linear expansion used in constructing the micromechanics model for

periodic multiphase materials called the Generalized Method of Cells (Paley and Aboudi, 1992).

The higher-order theory is an extension of this micromechanics model to materials without periodic

microstructures, and employs similar volume discretization and related notation.

The unknown coefficients associated with each term in the temperature field expansion are obtained by constructing a system of equations that satisfies the requirements of a standard boundaryvalue problem for the given field variable approximations. That is, the zeroth, first, and second

moments of the heat conduction equation are satisfied in a volumetric sense. The thermal and

heat flux interfacial continuity conditions within a given cell, as well as between a given cell and

adjacent cells, are imposed in an average sense together with the applied boundary conditions.

The following two subsections outline the thermal problem formulation and its solution based

on the volume discretization employed in the original higher-order theory in order to prepare the

stage for a generic cell-free, user-friendly formulation outlined in the second part of this section. A

detailed exposition was provided by Aboudi et al. (1999).

2.1.1 Thermal analysis

Heat conduction equation Assume that the functionally graded material model is subjected

to steady-state temperature or heat flux distribution on its bounding surfaces in the x2 ° x3 plane.

NASA/CR—2002-211910 9
Under steady-state heat conduction, the heat flux field in the material occupying subcell (Ø°) of

the (q; r)th generic cell, in the region Ø

Ø

Øx1

(β)

2

Ø

Ø

Ø ∙ 1

2h(q)

β ;

Ø

Ø

Øx1

(γ)

3

Ø

Ø

Ø ∙ 1

2 l

(r)

γ ; must satisfy:

@q(βγ)

2

@x1

(β)

2

+ @q(βγ)

3

@x1

(γ)

3

= 0 (1)

The components q

(βγ)

i of the heat flux vector in this subcell are obtained from the Fourierís heat

conduction law for anisotropic materials,

q

(βγ)

i = °k(βγ)ij

@T(βγ)

@x1

(·)

i

; (i; j = 2; 3; no sum) (2)

where T(βγ) is the temperature distribution in the subcell (Ø°) of the (q; r)th generic cell measured

with respect to a reference temperature, and k(βγ)

ij are the coefficients of thermal conductivity of

the material in the subcell (Ø°), with k(βγ)

ij = k(βγ)i ±ij (no sum on i) for orthotropic materials. No

sum is implied by the repeated Greek letters in the above and henceforth.

The solution to the heat conduction equation in each subcell of every generic cell is obtained

in conjuction with satisfaction of the continuity and boundary conditions in the manner described

next.

Heat flux continuity conditions The continuity of the heat flux vector q(βγ) at the interfaces separating adjacent subcells within the generic cell (q; r) is fulfilled by imposing the integral

conditions Z l

(r) γ /2

−l

(r) γ /2

q

(1γ)

2 j

(q,r)

xØ(1)

2 =h(q)1 /2

dx1

(γ)

3 =

Z l

(r) γ /2

−l

(r) γ /2

q

(2γ)

2 j

(q,r)

xØ(2)

2 =−h(q)2 /2

dx1

(γ)

3 (3)

Z h(q)

β /2

−h(q)

β /2

q

(β1)

3 j

(q,r)

xØ(1)

3 =l

(r)

1 /2

dx1

(β)

2 =

Z h(q)

β /2

−h(q)

β /2

q

(β2)

3 j

(q,r)

xØ(2)

3 =−l

(r)

2 /2

dx1

(β)

2 (4)

In addition to the above surface-averaged continuity conditions within the (q; r)th generic cell, the

heat flux continuity conditions at the interfaces between adjacent cells are given by

Z l

(r) γ /2

−l

(r) γ /2

q

(2γ)

2 j

(q,r)

xØ(2)

2 =h(q)2 /2

dx1

(γ)

3 =

Z l

(r) γ /2

−l

(r) γ /2

q

(1γ)

2 j

(q+1,r)

xØ(1)

2 =−h(q+1)1 /2

dx1

(γ)

3 (5)

Z h(q)

β /2

−h(q)

β /2

q

(β2)

3 j

(q,r)

xØ(2)

3 =l

(r)

2 /2

dx1

(β)

2 =

Z h(q)

β /2

−h(q)

β /2

q

(β1)

3 j

(q,r+1)

xØ(1)

3 =−l

(r+1)

1 /2

dx1

(β)

2 (6)

Thermal continuity conditions The thermal continuity conditions at the interfaces separating

adjacent subcells within the generic cell (q; r) are similar to the corresponding heat flux continuity

conditions, Z l

(r) γ /2

−l

(r) γ /2

T(1γ) j

(q,r)

xØ(1)

2 =h(q)1 /2

dx1

(γ)

3 =

Z l

(r) γ /2

−l

(r) γ /2

T(2γ) j

(q,r)

xØ(2)

2 =−h(q)2 /2

dx1

(γ)

3 (7)

Z h(q)

β /2

−h(q)

β /2

T(β1) j

(q,r)

xØ(1)

3 =l

(r)

1 /2

dx1

(β)

2 =

Z h(q)

β /2

−h(q)

β /2

T(β2) j

(q,r)

xØ(2)

3 =−l

(r)

2 /2

dx1

(β)

2 (8)

NASA/CR—2002-211910 10
Similarly, the thermal continuity conditions at the interfaces between neighboring cells are given

by

Z l

(r) γ /2

−l

(r) γ /2

T(2γ) j

(q,r)

xØ(2)

2 =h(q)2 /2

dx1

(γ)

3 =

Z l

(r) γ /2

−l

(r) γ /2

T(1γ) j

(q+1,r)

xØ(1)

2 =−h(q+1)1 /2

dx1

(γ)

3 (9)

Z h(q)

β /2

−h(q)

β /2

T(β2) j

(q,r)

xØ(2)

3 =l

(r)

2 /2

dx1

(β)

2 =

Z h(q)

β /2

−h(q)

β /2

T(β1) j

(q,r+1)

xØ(1)

3 =−l

(r+1)

1 /2

dx1

(β)

2 (10)

Boundary conditions The solution for the temperature field must satisfy the boundary conditions at top and bottom, left and right surfaces of the graded material, i.e., on the external surfaces

of the boundary cells q = 1; Nq; and r = 1; Nr. When q = 1 and q = Nq; the temperature in

the cells (1; r) and (Nq; r) at the x2 = 0 and x2 = H (i.e., left and right) surfaces, respectively,

must equal the applied temperature. As in the case of the interfacial continuity conditions, these

boundary conditions are imposed in an average sense,

Z l

(r) γ /2

−l

(r) γ /2

T(1γ) j

(1,r)

xØ(1)

2 =− 1

2h(1)

1

dx1

(γ)

3 =

Z l

(r) γ /2

−l

(r) γ /2

T(γ)

lef t dx1

(γ)

3 (11)

Z l

(r) γ /2

−l

(r) γ /2

T(2γ) j

(Nq,r)

xØ(2)

2 =+ 1

2h

(Nq)

2

dx1

(γ)

3 =

Z l

(r) γ /2

−l

(r) γ /2

T(γ)

right dx1

(γ)

3 (12)

where r = 1; ::::::; Nr:

Similarly, when r = 1 and r = Nr, the temperature in the cell (q; 1) and (q; Nr) at the x3 = 0

and x3 = L (top and bottom) surfaces, respectively, must equal the applied temperature,

Z h(q)

β /2

−h(q)

β /2

T(β1) j

(q,1)

xØ(1)

3 =− 1

2 l

(1)

1

dx1

(β)

2 =

Z h(q)

β /2

−h(q)

β /2

T(γ)

bot dx1

(β)

2 (13)

Z h(q)

β /2

−h(q)

β /2

T(β2) j

(q,Nr)

xØ(2)

3 =+ 1

2 l

(Nr)

2

dx1

(β)

2 =

Z h(q)

β /2

−h(q)

β /2

T(γ)

top dx1

(β)

2 (14)

where q = 1; :::::; Nq:

Alternatively, we can also apply mixed-boundary conditions involving temperature and heat

flux on different portions of the boundary.

2.1.2 Solution of the thermal problem

The temperature distribution in the subcell (Ø°) of the (q; r)th generic cell, measured with respect

to a reference temperature Tref ; is denoted by T(βγ)

: The temperature field is approximated by a

second-order expansion in the local coordinates x1

(β)

2 ; x1

(γ)

3 as follows (omitting the indices (q; r) for

notational simplicity):

T(βγ) = T(βγ)

(00) + 1x(β)

2 T(βγ)

(10) + 1x(γ)

3 T(βγ)

(01) +

1

2

(31x(β)2

2 °

h(q)2

β

4 )T(βγ)

(20) +

1

2

(31x(γ)2

3 ° l

(r)2

γ

4 )T(βγ)

(02) (15)

In the two-directional version of HOTFGM, the temperature field varies only in the x2°x3 plane, so

the temperature field expansion terms in the x1 direction are not included. In the above temperature

NASA/CR—2002-211910 11
field approximation, we can see clearly that there are five unknowns, where T(βγ)

(00) is the volumeaveraged temperature within the subcell, and T(βγ)

(mn) (m; n = 0; 1; 2 with m+n ∙ 2) are higher-order

terms which are determined from the conditions outlined previously. The reason for the absence of

the product terms x(.)

i x(.)j (i 6= j) is due to the averaging procedure employed in the higher-order

theory.

In summary, in each generic cell there are four subcells and in every subcell the five unknown

coefficients,

[T(00); T(10); T(01); T(20); T(02)]

(βγ)

describe the temperature field. Therefore, there are altogether 20NqNr unknown quantities which

must be determined for a composite with Nq rows and Nr columns of cells containing arbitrarily

specified materials within their subcells. These unknown coefficients are determined by first satisfying the heat conduction equation, as well as the first and second moment of this equation in

each subcell in a volumetric sense. Then the temperature and heat flux continuity conditions are

applied in an average sense at the interfaces separating adjacent subcells as well as adjacent cells.

Finally, the boundary conditions are imposed in an average sense. Following this procedure, we

obtain 20NqNr equations for the 20NqNr unknown coefficients in the temperature field expansion

as summarized below.

2.1.3 Heat conduction equation

In the course of satisfying the steady-state heat equation in a volumetric sense, the heat flux

quantities given below are defined for convenience

Q(βγ)

i(m,n) = 1

A(qr)

(βγ)

Z +h(q)

(β)/2

−h(q)

(β)/2

Z +l

(r)

(γ)/2

−l

(r)

(γ)/2

(1x(β)

2 )

m(1x(γ)

3 )

nq

(βγ)

i dx1

(β)

2 dx1

(γ)

3 (16)

where m; n = 0; 1; or 2 with m+n ∙ 2, and A(qr)

(βγ) = h(q)

β l

(r)

γ is the cross-sectional area of the subcell

(Ø°) in the (q; r)th generic cell. For example, when m = n = 0; Q(βγ)

i(00) is the volume-averaged

value of the heat flux component q

(βγ)

i in the subcell, while for other values of (m; n), higher-order

volume integrals of the heat flux components are obtained. These flux quantities can be evaluated

explicitly in terms of the coefficients T(βγ)

(mn) by performing the required volume integration using

Eqs. (2) and (15) in Eq. (16). This yields the following non-vanishing zeroth-order and first-order

heat fluxes expressed in terms of the unknown coefficients in the temperature field expansion.

Q(βγ)

2(00) = °k(βγ)

2 T(βγ)

(10) (17)

Q(βγ)

2(10) = °k(βγ)

2

h(q)2

β

4

T(βγ)

(20) (18)

Q(βγ)

3(00) = °k(βγ)

3 T(βγ)

(01) (19)

Q(βγ)

3(01) = °k(βγ)

3

l

(r)2

γ

4

T(βγ)

(02) (20)

Satisfaction of the zeroth, first and second moment of the steady-state heat equation, Eq. (1), in

volumetric sense results in the following four relationships among the first-order heat fluxes Q(βγ)

i(m,n)

NASA/CR—2002-211910 12
in the different subcells (Ø°) of the (q; r)th generic cell, where (Ø°) assumes all permutations of

the integers 1 and 2,

[Q(βγ)

2(10)=h2

β + Q(βγ)

3(01)=l2

γ]

(q,r) = 0 (21)

Heat flux continuity equations The equations that ensure heat flux continuity at the subcell

interfaces, as well as between individual cells, associated with the x2 and x3 directions, Eqs. (3)-(6)

are given by

[°12Q(1γ)

2(1,0)=h1 + Q(2γ)2(0,0) ° 6Q(2γ)2(1,0)=h2]

(q,r) ° [Q(2γ)

2(0,0) + 6Q(2γ)2(1,0)=h2]

(q−1,r) = 0 (22)

[°Q(1γ)

2(0,0) +

1

2

Q(2γ)

2(0,0) ° 3Q(2γ)2(1,0)=h2]

(q,r) +

1

2

[Q(2γ)

2(0,0) + 6Q(2γ)2(1,0)=h2]

(q−1,r) = 0 (23)

[°12Q(β1)

3(0,1)=l1 + Q(β2)3(0,0) ° 6Q(β2)3(0,1)=l2]

(q,r) ° [Q(β2)

3(0,0) + 6Q(β2)3(0,1)=l2]

(q,r−1) = 0 (24)

[°Q(β1)

3(0,0) +

1

2

Q(β2)

3(0,0) ° 3Q(β2)3(0,1)=l2]

(q,r) +

1

2

[Q(β2)

3(0,0) + 6Q(β2)3(0,1)=l2]

(q,r−1) = 0 (25)

Equations (22)-(25) provide eight relations among the zeroth-order and first-order heat fluxes. It

is observed that Eqs. (22)-(25) together with Eq. (21) can be expressed in terms of T(βγ) mn by using

Eqs. (17)-(20).

Thermal continuity equations The equations that ensure the continuity of temperature at the

interfaces between the subcells, as well as between the neighboring cells in the x2 and x3 directions,

obtained from Eqs. (7)-(8) and Eqs. (9)-(10), are as follows:

[T(1γ)

(00) +

1

2

h1T(1γ)

(10) +

1

4

h2

1T(1γ)

(20) ]

(q,r) = [T(2γ)

(00) ° 1

2

h2T(2γ)

(10) +

1

4

h2

2T(2γ)

(20) ]

(q,r) (26)

[T(2γ)

(00) +

1

2

h2T(2γ)

(10) +

1

4

h2

2T(2γ)

(20) ]

(q,r) = [T(1γ)

(00) ° 1

2

h1T(1γ)

(10) +

1

4

h2

1T(1γ)

(20) ]

(q+1,r) (27)

[T(β1)

(00) +

1

2

l1T(β1)

(01) +

1

4

l

2

1T(β1)

(02) ]

(q,r) = [T(β2)

(00) ° 1

2

l2T(β2)

(01) +

1

4

l

2

2T(β2)

(02) ]

(q,r) (28)

[T(β2)

(00) +

1

2

l2T(β2)

(01) +

1

4

l

2

2T(β2)

(02) ]

(q,r) = [T(β1)

(00) ° 1

2

l1T(β1)

(01) +

1

4

l

2

1T(β1)

(02) ]

(q,r+1) (29)

Equations (26)-(29) provide us with the remaining eight relations governing the temperature field

within interior generic cells.

Governing equations The heat conduction equation, Eq. (21), together with the heat flux

and thermal continuity conditions, Eqs. (22)-(25) and (26)-(29), respectively, form 20NqNr linear

algebraic equations which govern the 20NqNr microvariables T(βγ) mn in the four subcells (Ø°) of an

interior cell (q; r); q = 2; :::Nq ° 1; r = 2; :::Nr ° 1: For the boundary cells q = 1; Eq. (21) as well

as Eqs. (24)-(29) are still applicable. However, the heat flux continuity conditions between a given

cell and the preceding cell given by Eqs. (22) and (23) are not applicable. They are replaced by

the applied temperature or heat flux conditions at the left boundary x2 = 0; i.e., Eq. (11), and

the condition that the heat flux at the interface between subcell (1°) and (2°) of the cell (1; r) is

NASA/CR—2002-211910 13
continuous. For the boundary cells q = Nq; the previous equations are applicable except for Eq.

(27), which is replaced by the specific temperature applied at the right boundary x2 = H; i.e., Eq.

(12). Similar situation holds for the boundary cells r = 1 and r = Nr:

The governing equations for the interior and boundary cells form a system of 20NqNr algebraic

equations in the unknown coefficients T(βγ)

(mn)

: The solution of these equations determines the temperature distribution within the functionally graded composite subjected to the boundary conditions

specified by Eqs. (11)-(14). The final form of this system of equations is symbolically represented

below:

kT1 =t (30)

where the structural thermal conductivity matrix k1 contains information on the geometry and

thermal conductivities of the individual subcells (Ø°) in the NqNr generic cells (q; r). The thermal

coefficient vector T contains the unknown coefficients that describe the temperature fields in each

subcell. It is given by

T = [T(11)

11, ::::::::::::::::::; T(22)

NqNr ]

where

T(βγ)

qr = [T(00); T(10); T(01); T(20); T(02)]

(βγ)

(qr)

Finally, the thermal force vector t contains information on the boundary conditions.

2.2 User-Friendly Formulation of HOTFGM-2D

The motivation for the development of user-friendly thermal formulation of HOTFGM-2D is to

ultimately construct a basis for an efficient reformulation using the local/global conductivity matrix

approach, which will be outlined in the next section. User-friendly formulation of HOTFGM2D basically follows the solution technique of the original formulation which was outlined in the

preceding section. The major difference between these two formulations is simplification of the

equations through the elimination of the generic cell concept which also simplifies the heat flux

continuity conditions.

In the user-friendly formulation, therefore, the microstructure of the heterogeneous composite

is discretized into NβN ° subcells in the intervals 0 ∙ x2 ∙ H; 0 ∙ x3 ∙ L, i.e., the x2 ° x3 plane,

Figure 5. The arbitrary subcell (Ø°) has length hβ, and width lγ, where the Ø = 1; 2; :::::; Nβ;

° = 1; 2; :::::; Nγ: In the original formulation, there are five unknown quantities associated with each

subcell and four subcells within each generic cell, therefore, 20NqNr unknown coefficients must be

determined by establishing a 20NqNr system of equations involving field equations and continuity

conditions, together with the imposed thermal or heat flux boundary conditions. In the user-friendly

formulation, there are still five unknown quantities in the temperature field expansion associated

with every subcell, because the same temperature expansion is employed. Since the total number

of subcells is NβNγ, 5NβNγ equations must be constructed for the determination of the unknown

quantities. Although the number of unknowns in each subcell and the total number of equations

remain the same in the user-friendly formulation (since 2Nq = Nβ and 2Nr = Nγ); the simplified

notation leads to a more systematic method of assembling the final system of equations for the

determination of the unknown coefficients in the temperature field expansion. More importantly,

it also provides a basis for the efficient reformulation described in the next section.

NASA/CR—2002-211910 14
Figure 5. Schematic of the user-friendly HOTFGM formulation geometry.

2.2.1 Thermal analysis

As before, the temperature distribution T(βγ) in the subcell (Ø°) is approximated by

T(βγ) = T(βγ)

(00) + 1x(β)

2 T(βγ)

(10) + 1x(γ)

3 T(βγ)

(01) +

1

2

(31x(β)2

2 °

h2

β

4 )T(βγ)

(20) +

1

2

(31x(γ)2

3 ° l

2

γ

4 )T(βγ)

(02) (31)

Heat conduction equation As in the original formulation of HOTFGM-2D, under steady-state

heat conduction, the heat flux field in the material occupying subcell (Ø°) in the region jx11j < 1; Ø

Ø

Øx1

(β)

2

Ø

Ø

Ø ∙ 1

2hβ;

Ø

Ø

Øx1

(γ)

3

Ø

Ø

Ø ∙ 1

2 lγ must satisfy Eq. (1) in a volumetric sense,

Z

V (βγ)

(

@q(βγ)

2

@x1

(β)

2

+ @q(βγ)

3

@x1

(γ)

3

)dV (βγ) = 0 (32)

The heat flux components q

(βγ)

2 and q

(βγ)

3 can be expressed in terms of the unknown coefficients

in the temperature field expansion, Eg. (31), using Eq. (2),

q

(βγ)

2 = °k(βγ)2

@T(βγ)

@x1

(β)

2

= °k(βγ)

2 (T(βγ)

(10) + 31x(β)

2 T(βγ)

(20) ) (33)

q

(βγ)

3 = °k(βγ)3

@T(βγ)

@x1

(γ)

3

= °k(βγ)

3 (T(βγ)

(01) + 31x(γ)

3 T(βγ)

(02) ) (34)

NASA/CR—2002-211910 15
Therefore, substituting Eqs. (33) and (34) into the steady-state heat conduction equation, Eq.

(32), yields,

k(βγ)

2 T(βγ)

(20) + k(βγ)

3 T(βγ)

(02) = 0 (35)

This is the same equation as Eq.(21) in the original formulation at the subcell level, after replacing

the temperature microvariables with the volume-averaged heat fluxes using Eqs. (18) and (20).

Heat flux continuity conditions between adjacent subcells Using the new generic cell-free

notation, the heat flux continuity condition between adjacent subcellís along the x2 direction is

Z lγ/2

−lγ/2

q

(β,γ)

2 jxØ2=hβ/2 dx1

(γ)

3 =

Z lγ/2

−lγ/2

q

(β+1,γ)

2 jxØ2=−hβ+1/2 dx1

(γ)

3 (36)

Substituting for the heat flux components in the x2 direction using Eq. (33), the above equation

becomes:

°k(β,γ)

2 (T(β,γ)

(10) lγ +

3hβ

2

T(β,γ)

(20) lγ) = °k(β+1,γ)

2 (T(β+1,γ)

(10) lγ ° 3hβ+1

2

T(β+1,γ)

(20) lγ) (37)

Similarly, the heat flux continuity condition along the x3 direction is

Z hβ/2

−hβ/2

q

(β,γ)

3 jxØ3=lγ/2 dx1

(β)

2 =

Z hβ/2

−hβ/2

q

(β,γ+1)

3 jxØ3=−lγ+1/2 dx1

(β)

2 (38)

which, upon using Eq. (34), becomes

°k(β,γ)

3 (T(β,γ)

(01) hβ +

3lγ

2

T(β,γ)

(02) hβ) = °k(β,γ+1)

3 (T(β,γ+1)

(01) hβ ° 3lγ+1

2

T(β,γ+1)

(02) hβ) (39)

Thermal continuity conditions between adjacent subcells Using the new notation, the

thermal continuity condition at adjacent subcells interfaces along the x2 direction is,

Z lγ/2

−lγ/2

T(βγ) jxØ2=hβ /2 dx1

(γ)

3 =

Z lγ/2

−lγ/2

T(β+1,γ) jxØ2=−hβ+1 /2 dx1

(γ)

3 (40)

Upon substituting for the temperature field defined by Eq. (31) and performing the required

integration, the above condition reduces to

T(βγ)

(00) lγ +

hβlγ

2

T(βγ)

(10) +

h2

βlγ

4

T(βγ)

(20) = T(β+1,γ)(00) lγ ° hβ+1lγ

2

T(β+1,γ)

(10) +

h

2

β+1lγ

4

T(β+1,γ)

(20) (41)

Similarly, in the x3 direction, the thermal continuity condition is,

Z hβ/2

−hβ/2

T(βγ) jxØ3=lγ /2 dx1

(β)

2 =

Z hβ/2

−hβ/2

T(β,γ+1) jxØ3=−lγ+1 /2 dx1

(β)

2 (42)

which, upon using Eq. (31), becomes

T(βγ)

(00) hβ +

lγhβ

2

T(βγ)

(01) +

l

2

γhβ

4

T(βγ)

(02) = T(β,γ+1)(00) hβ ° lγ+1 hβ

2

T(β,γ+1)

(01) +

l

2

γ+1hβ

4

T(β,γ+1)

(02) (43)

NASA/CR—2002-211910 16
Boundary conditions The temperature or heat flux at the exterior surfaces of the boundary

subcells (Ø; 1) and (Ø;Nγ) at x3 = 0 and x3 = L, Figure 6, respectively, must equal the applied

boundary condition (given temperature or given heat flux) in an average sense,

Figure 6. Boundary subcells.

Z hβ/2

−hβ/2

T(β,1)jxØ3=−l1/2dx1

(β)

2 = T(β,1)

(00) ° T(β,1)(01) l1=2 + T(β,1)(02) l

2

1=4 = Z hβ/2

−hβ/2

T(β)

bot dx1

(β)

2 = T1(β)

bot (44)

or

Z hβ/2

−hβ/2

q

(β,1)

3 jxØ3=−l1/2dx1

(β)

2 = °k(β,1)3 (T(β,1)

(01) ° 3l1T(β,1)(02) =2) = Z hβ/2

−hβ/2

q

(β)

bot dx1

(β)

2 = 1q

(β)

bot (45)

and

Z hβ/2

−hβ/2

T(β,Nγ)jxØ3=lNγ /2dx1

(β)

2 = T(β,Nγ)

(00) + T(β,Nγ)(01) lNγ =2 + T(β,Nγ)(02) l

2

Nγ =4 = Z hβ/2

−hβ/2

T(β)

top dx1

(β)

2 = T1(β)

top

(46)

or

Z hβ/2

−hβ/2

q

(β,Nγ)

3 jxØ3=lNγ /2dx1

(β)

2 = °k(β,Nγ)

3 (T(β,Nγ)

(01) + 3lNγ T(β,Nγ)(02) =2) = Z hβ/2

−hβ/2

q

(β)

top dx1

(β)

2 = 1q

(β)

top (47)

Similarly the temperature or heat flux in the boundary subcells (1; °) and (Nβ; °) at x2 = 0

and x2 = H must be equal to the applied boundary condition (given temperature or given heat

flux) in an average sense,

Z lγ/2

−lγ/2

T(1,γ)jxØ2=−h1/2dx1

(γ)

3 = T(1,γ)

(00) ° T(1,γ)(10) h1=2 + T(1,γ)(20) h2

1=4 = Z lγ/2

−lγ/2

T(γ)

lef tdx1

(γ)

3 = T1(γ)

lef t (48)

or

NASA/CR—2002-211910 17
Z lγ/2

−lγ/2

q

(1,γ)

2 jxØ2=−h1/2dx1

(γ)

3 = °k(1,γ)2 (T(1,γ)

(10) ° 3h1T(1,γ)(20) =2) = Z lγ/2

−lγ/2

q

(γ)

lef tdx1

(γ)

3 = 1q

(γ)

lef t (49)

and

Z lγ/2

−lγ/2

T(Nβ,γ)jxØ2=hNβ /2 dx1

(γ)

3 = T(Nβ,γ)

(00) +T(Nβ,γ)(10) hNβ =2 +T(Nβ,γ)(20) h2

Nβ =4 = Z lγ/2

−lγ/2

T(γ)

rightdx1

(γ)

3 = T1(γ)

right

(50)

or

Z lγ/2

−lγ/2

q

(Nβ,γ)

2 jxØ2=hNβ /2dx1

(γ)

3 = °k

(Nβ,γ)

2 (T(Nβ,γ)

(10) +3hNβ T(Nβ,γ)(20) =2) = Z lγ/2

−lγ/2

q

(γ)

rightdx1

(γ)

3 = 1q

(γ)

right (51)

2.2.2 Assembly of the governing equations

As stated before, there are five unknown quantities in the temperature field expansion associated

with an arbitrary subcell (Ø°). They can be represented by the vector T(βγ) defined as follows,

T(βγ) = [T(00); T(10); T(01); T(20); T(02)]

(βγ) (52)

Therefore, we need to construct five equations for every subcell to solve for the unknown coefficients.

Two equations involve temperature continuity in the x2 and x3 directions, another two involve heat

flux continuity in the x2 and x3 directions, and the remaining equation is the steady-state heat

conduction equation. These will be arranged into an organized system of equations along each row

sequence and column sequence of subcells representing the functionally graded material.

Let the assembly of the system of equations proceed first from left to right (in the x2 direction)

along each row of subcells starting from the first row, Figure 7. Subsequently, we proceed from

bottom to top (in the x3 direction) along each column of subcells starting from the first column.

Finally, we take care of the steady-state heat equation within each subcell along either a row

sequence or a column sequence. In summary, we obtain two equations at each interior interface

along either direction after rearrangement, one involving temperature continuity and the other

associated with heat flux continuity. We also obtain one steady-state heat equation for each subcell

along either a row sequence or a column sequence for a total of five equations.

Figure 7. Governing equation assembly.

NASA/CR—2002-211910 18
Continuity and boundary conditions along x2 direction For a given °; the left boundary

condition, Eq. (48) or (49), can be rewritten in matrix form as follows,

£

KB

1γ

§ h

T(1γ)

i

= T1(γ)

lef t or q1

(γ)

lef t (53)

where KB

1γ depends on the applied boundary condition, and T1(γ)

lef t and q1

(γ)

lef t are the surface integrals

of the temperature and heat flux on the external surface of the boundary subcell (1; °) given by

the right side of Eqs. (48) and (49). If the applied boundary condition is temperature, then

KB

1γ = £

1 ° h1=2 0 h2

1=4 0§

Alternatively, if the applied boundary condition is heat flux, then

KB

1γ =

h

0 ° k(1,γ)

2 0 ° 3=2k(1,γ)2 h1 0

i

The definition of T(βγ) is given by Eq. (52).

Proceeding into the interior, we can rewrite the temperature and heat flux continuity conditions

at the interface separating two adjacent subcells (Ø; °) and (Ø + 1; °) in the interior, Eqs. (37) and

(41), along the x2 direction as follows,

" KT +

βγ KT −β+1,γ

KH+

βγ KH−β+1,γ

# ∙ T(βγ)

T(β+1,γ)

 ̧

=

∙ 0

0

 ̧

(54)

where

KT +

βγ =

h

1 hβ=2 0 h2

β=4 0i

KT −

β+1,γ =

h

°1 hβ+1=2 0 ° h2

β+1=4 0i

KH+

βγ =

h

0 ° k(βγ)

2 0 ° 3=2hβk(βγ)2 0

i

KH−

β+1,γ =

h

0 k(β+1,γ)

2 0 ° 3=2hβ+1k(β+1,γ)2 0

i

and Ø = 1; 2; :::::; Nβ ° 1:

Finally, the right boundary condition, Eq. (50) and (51), can be rewritten in matrix form as,

h

KB

Nβγ

i hT(Nβ,γ)i

= T1(γ)

right or q1

(γ)

right (55)

where KNβ,γ depends on the applied boundary condition, and T1(γ)

right and q1

(γ)

right are the surface

integrals of the temperature and heat flux on the external surface of the boundary subcell (Nβ; °)

given by the right side of Eqs. (50) and (51). If the boundary condition is temperature, then

KNβ,γ =

h

1 hNβ =2 0 h2

Nβ =4 0i

Alternatively, if the applied boundary condition is heat flux, then

KNβ,γ =

h

0 ° k

(Nβ,γ)

2 0 ° 3=2k

(Nβ,γ)

2 hNβ 0

i

NASA/CR—2002-211910 19
Therefore, for an arbitrary row °, let Ø = 1; :::; Nβ, Figure 7, and apply the above equations

from the left boundary to the right boundary, respectively. This procedure produces a system of

equations defined by the matrix K(#

β)

γ ; which contains information on the left boundary condition,

temperature and heat flux continuity conditions at the interfaces between adjacent subcells, and

the right boundary condition, respectively, along the row sequence. The structure of K(#

β)

γ , whose

size is 2Nβ £ 5Nβ, is shown below,

K(#

β)

γ =

2

6

6

6

6

6

6

6

6

6

6

6

6

4

KB

1γ 00 0 0 0 0

:::::: :::::: 0 0 0 00

0 0KT +

βγ KT −β+1,γ 0 00

0 0KH+

βγ KH−β+1,γ 0 00

0 0 0KT +

β+1,γ

KT −

β+2,γ 0 0

0 0 0KH+

β+1,γ

KH−

β+2,γ 0 0

00 0 0 0 ::::::: ::::::

0 0 0 0 0 0KB

Nβ,γ

3

7

7

7

7

7

7

7

7

7

7

7

7

5

Following the same procedure applied to every row, ° = 1; 2; ::::Nγ; we obtain a system of

equations defined by the matrix A; consisting of the individual matrices K(#

β)

γ ; that has the form,

A =

2

6

6

6

6

6

6

4

K(#

β)

1 00 0 0

0 ::::: 000

0 0 K(#

β)

γ 0 0

000 :::::: 0

000 0 K(#

β)

Nγ

3

7

7

7

7

7

7

5

Since the matrix A consists of Nγ matrices K(#

β)

γ placed along the main diagonal, its size is 2NβN °£

5NβNγ.

Continuity and boundary conditions along x3 direction Repeating the preceding procedure

along the x3 direction, the bottom boundary condition, Eq. (44) or (45), for a given Ø; can be

rewritten in matrix form as follows,

h

K^ B

β1

i hT(β1)i

= T1(β)

bot or q1

(β)

bot (56)

As in the case of the boundary matrix KB

1γ in A; K^ B

β1 also can be obtained in two ways. When

temperature is applied to the boundary, then we get

K^ B

β1 =

h

1 0 ° l1=2 0 l

2

1=4

i

Alternatively, when the boundary condition involves heat flux, we get

K^ B

β1 =

h

0 0 ° k(β1)

3 0 3=2k(β1)3 l1

i

T1(β)

bot and q1

(β)

bot are the surface integrals of the temperature and heat flux on the external surface of

the boundary subcell (Ø; 1) given by the right side of Eqs. (44) and (45).

NASA/CR—2002-211910 20
Proceeding into the interior; we can assemble the temperature and heat flux continuity conditions at the interface separating two adjacent subcells (Ø; °) and (Ø; ° +1) in the interior, Eqs. (39)

and (43), along the x3 direction as follows,

" K^ T +

βγ K^ T −β,γ+1

K^ H+

βγ K^ H−β,γ+1 # ∙ T(βγ)

T(β,γ+1)  ̧

=

∙ 0

0

 ̧

(57)

where

K^ T +

βγ = £

1 0 lγ=2 0 l

2

γ=4

§

K^ T −

β,γ+1 = £

°1 0 lγ+1=2 0 ° l

2

γ+1=4

§

K^ H+

βγ =

h

0 0 ° k(βγ)

3 0 ° 3=2k(βγ)3 lγ

i

K^ H−

β,γ+1 =

h

0 0 k(β,γ+1)

3 0 ° 3=2k(β,γ+1)3 lγ+1i

and ° = 1; ::::; Nγ ° 1:

Finally, we rewrite the top boundary condition, Eq. (46) or (47), in matrix form as,

h

K^ B

β,Nγ

i hT(β,Nγ)i

= T1(β)

top or q1

(β)

top (58)

where T1(β)

top and q1

(β)

top are the surface integrals of the temperature and heat flux on the external

surface of the boundary subcell (Ø;Nγ) given by the right side of Eqs. (46) and (47). If the applied

boundary condition is temperature defined as T1(β)

top , then

K^ B

β,Nγ =

h

1 0 lNγ =2 0 l

2

Nγ =4

i

Similarly, if the boundary is subjected to heat flux given by the form q1

(β)

top , then we obtain

K^ B

β,Nγ =

h

0 0 ° k(β,Nγ)

3 0 3=2k(β,Nγ)3 lNγ

i

Following these steps, for an arbitrary column Ø, letting ° = 1 ! Nγ, Figure 7, we construct a

system of equations defined by the matrix Kà (#γ)

β which contains information on the bottom boundary,

interfacial temperature and heat flux continuity conditions, and the top boundary, respectively,

along the column sequence (i.e., in the x3 direction).

K^ (#γ)

β =

2

6

6

6

6

6

6

6

6

6

6

6

6

4

: K^ B

β1 :: 0 :: 0 :: 0 :: 0 :::

: ¢ :: 0 :: 0 :: 0 :: 0 :::

: 0 :: K^ T +

βγ :: K^ T −β,γ+1 :: 0 :: 0 :::

: 0 :: K^ H+

βγ :: K^ H−β,γ+1 :: 0 :: 0 :::

: 0 :: 0 :: K^ T +

β,γ+1 :: K^ T −β,γ+2 :: 0 :::

: 0 :: 0 :: K^ H+

β,γ+1 :: K^ H−β,γ+2 :: 0 :::

: 0 :: 0 :: 0 :: 0 :: ¢ :::

: 0 :: 0 :: 0 :: 0 :: K^ B

β,Nγ :::

3

7

7

7

7

7

7

7

7

7

7

7

7

5

NASA/CR—2002-211910 21
In the above matrix, the designation : indicates Ø ° 1 0 matrices, the designation :: indicates

Nβ ° 1 0 matrices, and the designation ::: indicates Nβ ° Ø 0 matrices, with the zero matrices

having the dimensions [1 £ 5]. The size of the K^ (#γ)

β matrices is 2Nγ £ 5NβNγ.

Extending Kà (#γ)

β to all Ø, Ø = 1; 2; ::::; Nβ; we obtain the matrix B that has the form

B =

2

6

6

6

6

6

6

4

K(#γ)

1

......

K(#γ)

β

......

K(#γ)

Nβ

3

7

7

7

7

7

7

5

with B having the dimensions 2NβNγ £ 5NβNγ:

Heat conduction equation in x2 or x3 direction Finally, we write the steady-state heat

conduction equation, Eq. (35), for an arbitrary (Ø°) subcell as follows

h

K∗(β,γ)

γ

i hT(βγ)i

= 0

where

K∗(β,γ)

γ = [0 0 0 ° 3k(βγ)

2 ° 3k(βγ)3 ]

Assuming that the assembly proceeds from left to right along a row sequence, we obtain the matrix

K∗(#

β)

γ with dimensions Nβ £ 5Nβ for a given ° row of subcells

K∗(#

β)

γ =

2

6

6

6

6

6

4

K∗(1,γ)

γ 000 0

0 :::: 00 0

0 0K∗(β,γ)

γ 0 0

000 :::: 0

0 0 0 0K∗(Nβ,γ)

γ

3

7

7

7

7

7

5

where Ø = 1; 2; :::::Nβ along the diagonal, respectively.

Extending the procedure to all °, we obtain the final system of equations that satisfies the heat

conduction equation in every subcell. This system of equation is represented by the matrix C,

which contains the individual matrices K∗(#

β)

γ along its diagonal as shown below.

C =

2

6

6

6

6

6

6

4

K∗(#

β)

1 00 0 0

0 ::::: 000

0 0 K∗(#

β)

γ 0 0

000 :::::: 0

000 0 K∗(#

β)

Nγ

3

7

7

7

7

7

7

5

where ° = 1; 2; ::::::Nγ, respectively, along the diagonal. Since the matrix C consists of Nγ matrices

K∗(#

β)

γ placed along the diagonal, its size is NβNγ £ 5NβNγ:

NASA/CR—2002-211910 22
Governing equations The governing equations for the interior and boundary subcells form a

system of 5NβN ° algebraic equations in the 5NβNγ unknown coefficients T(βγ)

(mn)

. This system of

equations can be written symbolically as

kT = t (59)

where k is the structural thermal conductivity matrix. It contains information on the geometry

and thermal conductivity of every subcell. The matrix k is made up of three submatrices as follows

k =

2

4

A

B

C

3

5 (60)

T is the thermal coefficient vector. It contains the unknown coefficients that describe the temperature field in all the subcells arranged as follows

T = [T1; :::::;Tγ; ::::::TNγ ] (61)

where

Tγ = [T(1γ)

; ::::::; T(βγ); ::::::T(Nβγ)] (62)

with T(βγ) defined previously in Eq. (52). Finally, t is the thermal force vector. It contains

information on the boundary conditions specified either in terms of temperatures or heat fluxes.

In summary, the elimination of generic cells in the user-friendly formulation has resulted in a

system of equations which is easier to visualize and implement in a computer code. Most importantly, the new formulation based on subcells as the basic building blocks sets the stage for an

efficient reformulation presented in the next section.

3 Reformulation

The user-friendly formulation of HOTFGM-2D, which was presented in the preceding section,

illustrates the feasibility of eliminating the distinction between generic cells and subcells. It should

be noted that the simplified geometric model played an important role in the entire analytical

process. In other words, the basic building block defined by the subcell itself without a generic

cell, which was used to construct the heterogeneous composite, made the indices simpler and the

derivation process more straightforward. More importantly, formulating the higher-order theory

solely in terms of subcells, which play the role of generic subvolumes, reduces it to a standard

thermal boundary-value problem on a discretized domain, thereby facilitating the search for a

more efficient solution method.

As discussed by Pindera (1991), the local/global stiffness matrix method is an efficient algorithm

for solving mixed boundary-value problems in composite mechanics. This method is particularly

well suited to problems involving layered composite materials or structures that require satisfaction

of both continuity of tractions and displacements in mechanical problems, or temperature and

heat flux in thermal problems, along common interfaces. The method thus far has been applied

to layered media composed of continuous layers in Cartesian coordinates, or circular concentric

rings/cylinders in polar coordinates. Examples of the methodís application include axisymmetric

and axial shear problems involving multiple concentric cylinders, Pindera et al. (1993) and Williams

and Pindera (1997), composite tube problems, Salzar et al. (1996), a layered composite cylinder

subjected to diametral loading, Davison et al. (1994), and contact problems, Pindera and Lane

(1993) and Urquhart and Pindera (1994). The local/global stiffness matrix method is based on

NASA/CR—2002-211910 23
a systematic construction of a global stiffness matrix (mechanical problem) or global conductivity

matrix (thermal problem) for the entire composite structure in terms of local stiffness matrices

(mechanical problem) or local conductivity matrices (thermal problem) of individual layers. Herein,

we extend this method to HOTFGM-2D where the volume discretization involves subvolumes of

finite dimensions in the x2 ° x3 plane defined by the subcells (Ø°). Since we only discuss the

thermal problem here, we will call it henceforth the local/global conductivity matrix method.

In the local/global conductivity matrix method, the local conductivity matrix provides a connection between the surface-averaged heat flux and the corresponding surface-averaged temperature

components at top, bottom, left and right surfaces of a given subcell. The assembly of local conductivity matrices into the global conductivity matrix is based on the direct enforcement of continuity

conditions along the subcell interfaces, and the interfacial temperatures are considered as the basic

unknown variables in the reformulation. Therefore, for a volume discretized into NβN ° subcells,

the total number of interfacial surface-averaged temperatures will be (Nβ + 1)Nγ + (Nγ + 1)Nβ;

which includes both the interior subcell interfaces and the external surfaces of boundary subcells.

For a given (Ø°) subcell, the local conductivity matrix can be expressed as follows,

2

6

6

4

°Q1−

2

Q1+

2

°Q1−

3

Q1+

3

3

7

7

5

(βγ)

=

2

6

6

4

∙11 ∙12 ∙13 ∙14

∙21 ∙22 ∙23 ∙24

∙31 ∙32 ∙33 ∙34

∙41 ∙42 ∙43 ∙44

3

7

7

5

(βγ) 2

6

6

4

T1−

2

T1+

2

T1−

3

T1+

3

3

7

7

5

(βγ)

(63)

It relates the surface-averaged heat flux components °(Q1−

2 )(βγ)

;(Q1+

2 )(βγ)

; °(Q1−

3 )(βγ)

;(Q1+

3 )(βγ) on

the left, right, top and bottom surfaces of a given subcell to their corresponding surface-averaged

temperatures (T1−

2 )(βγ)

;(T1+

2 )(βγ)

;(T1−

3 )(βγ)

;(T1+

3 )(βγ)

; Figure 8. The superscripts + and - refer to the

corresponding heat flux or temperature acting on the right/top surface and left/bottom surface of

a given subcell. The subscripts 2 and 3 indicate the x2 and x3 directions, respectively. The minus

signs in front of the heat flux components Q1−

2 ; Q1−3 indicate that the direction of the heat flux vector

coincides with the direction of the unit normal vector associated with the left and bottom surfaces

of the (Ø°) subcell.

Figure 8. Definition of surface-averaged temperatures and heat fluxes.

NASA/CR—2002-211910 24
Using the Fourierís heat conduction law for orthotropic materials, Eqs. (33)-(34),

q

(βγ)

2 = °k(βγ)2

@T(βγ)

@x1

(β)

2

; q(βγ)

3 = °k(βγ)3

@T(βγ)

@x1

(γ)

3

the definitions of surface-averaged heat fluxes and temperatures in Eq. (63) are:

(Q1+

2 )

(βγ) = 1

l° Z lγ/2

−lγ/2

°k(βγ)

2

@T(βγ)

@x1

(β)

2

j

xØ(β)

2 =hβ/2 dx1

(γ)

3 (64)

(Q1−

2 )

(βγ) = 1

l° Z lγ/2

−lγ/2

°k(βγ)

2

@T(βγ)

@x1

(β)

2

j

xØ(β)

2 =−hβ/2 dx1

(γ)

3 (65)

(Q1+

3 )

(βγ) = 1

hβ

Z hβ/2

−hβ/2

°k(βγ)

3

@T(βγ)

@x1

(γ)

3

j

xØ(γ)

3 =lγ/2 dx1

(β)

2 (66)

(Q1−

3 )

(βγ) = 1

hβ

Z hβ/2

−hβ/2

°k(βγ)

3

@T(βγ)

@x1

(γ)

3

j

xØ(γ)

3 =−lγ/2 dx1

(β)

2 (67)

and

(T1+

2 )

(βγ) = 1

lγ

Z lγ/2

−lγ/2

T(βγ)

2 j

xØ(β)

2 =hβ/2 dx1

(γ)

3 (68)

(T1−

2 )

(βγ) = 1

lγ

Z lγ/2

−lγ/2

T(βγ)

2 j

xØ(β)

2 =−hβ/2 dx1

(γ)

3 (69)

(T1+

3 )

(βγ) = 1

hβ

Z hβ/2

−hβ/2

T(βγ)

3 j

xØ(γ)

3 =lγ/2 dx1

(β)

2 (70)

(T1−

3 )

(βγ) = 1

hβ

Z hβ/2

−hβ/2

T(βγ)

3 j

xØ(γ)

3 =lγ/2 dx1

(β)

2 (71)

The detailed derivation of the elements of the local conductivity matrix in Eq. (63) will be outlined

in the following section.

Imposition of continuity of heat flux and temperature along common subcell interfaces, together

with the external boundary conditions, gives rise to a system of equations in the unknown interfacial

temperatures along the x2 and x3 directions, respectively. The continuity of interfacial heat flux is

ensured by requiring that the sum of the heat fluxes acting on the interface separating (Ø; °) and

(Ø + 1; °) subcells and the interface separating (Ø; °) and (Ø; ° + 1) subcells be zero, i.e.,

(Q1+

2 )

(β,γ) + (°Q1−

2 )

(β+1,γ) = 0; (72)

(Q1+

3 )

(β,γ) + (°Q1−

3 )

(β,γ+1) = 0 (73)

where Ø = 1; 2; ::::Nβ ° 1 and ° = 1; 2; ::::Nγ ° 1: The thermal continuity is directly enforced by

setting the surface-averaged temperatures along interfaces separating adjacent subcells to common

values,

(T1+

2 )

(β,γ) = (T1−

2 )

(β+1,γ) = T1(β+1,γ)

2 (74)

(T1+

3 )

(β,γ) = (T1−

3 )

(β,γ+1) = T1(β,γ+1)

3 (75)

NASA/CR—2002-211910 25
The system of equations is constructed by applying Eqs. (72)-(75) in conjunction with Eq. (63)

to each interface, starting first with the left boundary in the first row of subcells where the boundary

condition is prescribed and then progressing in the x2 direction to each interface separating adjacent

interior subcells, and then finally the right boundary with the prescribed boundary condition. This

procedure is repeated along every row of subcells. We then repeat the same procedure in the x3

direction, starting at the bottom of the first column and proceeding up, and then repeating it

for each subsequent column. At each interface, the thermal continuity expressed in terms of the

common interfacial temperatures defined by Eqs. (74) and (75) is directly enforced in the heat

flux continuity conditions given by Eqs. (72) and (73) through the use of Eq. (63). Actually, the

assembly of the global conductivity matrix for the entire structure is carried out by superposing

local conductivity matrices of individual subcells along x2 and x3 directions, respectively, in the

manner that will be described in Section 3.3.

Since the continuity of interfacial temperatures is enforced directly through Eqs. (74) and

(75) along x2 and x3 directions, respectively, the above procedure results in the elimination of a

substantial number of redundant interfacial continuity equations which are retained in the original

formulation. Thus the reformulation simplifies the computational procedure for thermal boundaryvalue problems by reducing the size of the global conductivity matrix, which can be large for a large

number of subcells. Furthermore, the global conductivity matrix can be generated systematically

and automatically, thereby facilitating numerical implementation for any number of subcells.

The above presentation illustrates the general formulation process of the local/global conductivity matrix method in solving thermal boundary-value problems. In the following, the detailed

derivation of how to construct the local conductivity matrix based on the generic cell-free formulation of HOTFGM-2D problem is outlined. It is should be noted that herein the thermal conductivity

coefficients k(βγ)

2 and k(βγ)3 are assumed to vary with local subcell coordinates. The same problem

but with constant thermal conductivity coefficients was investigated by Bansal (2002).

3.1 Local Conductivity Matrix: Variable Thermal Conductivity

In the present work, the thermal conductivity coefficients k(βγ)

2 and k(βγ)3 in each (Ø°) subcell are

assumed to vary linearly with the local subcell coordinates x1

(β)

2 and x1

(γ)

3 as follows

k(βγ)

2 = k(βγ)20 + k(βγ)22 x1

(β)

2 + k(βγ)23 x1

(γ)

3 (76)

k(βγ)

3 = k(βγ)30 + k(βγ)32 x1

(β)

2 + k(βγ)33 x1

(γ)

3 (77)

We note that k(βγ)

23 and k(βγ)32 cannot be completely arbitrary as they do not explicitly appear in

the surface-averaged heat fluxes and the volume-averaged heat conduction equation derived in the

sequel. This is a direct result of the averaging technique employed in the higher-order theory.

Therefore, we must impose the following constraint for consistency: k(βγ)

23 = k(βγ)32 = 0. In other

words, we obtain one-dimensional variation of k(βγ)

2 and k(βγ)3 along the respective directions, if

linear thermal conductivity variation is assumed.

As before, the temperature field in the subcell (Ø°) is approximated by a second-order expansion

in the local coordinates x1

(β)

2 ; x1

(γ)

3 in the following way (see Eq. (31))

T(βγ) = T(βγ)

(00) + 1x(β)

2 T(βγ)

(10) + 1x(γ)

3 T(βγ)

(01) +

1

2

(31x(β)2

2 °

h2

β

4 )T(βγ)

(20) +

1

2

(31x(γ)2

3 ° l

2

γ

4 )T(βγ)

(02)

NASA/CR—2002-211910 26
Using Eqs. (68)-(71), the surface-averaged temperatures (T1+

2 )(βγ)

;(T1−

2 )(βγ)

;(T1+

3 )(βγ)

;(T1−

3 )(βγ) are

expressed in terms of the subcell temperature microvariables T(βγ) mn as follows

(T1+

2 )(βγ) = T(βγ)

(00) +

hβ

2

T(βγ)

(10) +

h2

β

4

T(βγ)

(20)

(T1−

2 )(βγ) = T(βγ)

(00) ° hβ

2

T(βγ)

(10) +

h2

β

4

T(βγ)

(20)

(T1+

3 )(βγ) = T(βγ)

(00) +

lγ

2

T(βγ)

(01) +

l

2

γ

4

T(βγ)

(02)

(T1−

3 )(βγ) = T(βγ)

(00) ° lγ

2

T(βγ)

(01) +

l

2

γ

4

T(βγ)

(02)

Adding and subtracting the above first two equations involving (T1+

2 )(βγ) and (T1−2 )(βγ)

; we obtain

expressions for T(βγ)

(10) and T(βγ)(20) in terms of (T1+

2 )(βγ)

, (T1−

2 )(βγ) and T(βγ)

(00) . The matrix form of these

two expressions become

∙ T(10)

T(20)  ̧(βγ)

=

2

6

6

4

° 1

hβ

1

hβ

2

h2

β

2

h2

β

3

7

7

5

∙ T1−

2

T1+

2

 ̧(βγ)

° 4

h2

β

∙ 0

1

 ̧

T(βγ)

(00) (78)

Similarly, T(βγ)

(01) and T(βγ)(02) can be expressed in terms of (T1−

3 )(βγ)

;(T1+

3 )(βγ) and T(βγ)

(00) by applying

the same technique to the above last two equations,

∙ T(01)

T(02)  ̧(βγ)

=

2

6

6

4

° 1

lγ

1

lγ

2

l2

γ

2

l2

γ

3

7

7

5

∙ T1−

3

T1+

3

 ̧(βγ)

° 4

l2

γ

∙ 0

1

 ̧

T(βγ)

(00) (79)

The next step is to evaluate the surface-averaged heat flux components in terms of the microvariables T(βγ)

(mn)

. Integrating Eqs. (64)-(67) with k(βγ)

2 and k(βγ)3 given by Eqs. (76) and (77),

the surface-averaged heat fluxes (Q1+

2 )(βγ)

;(Q1−

2 )(βγ)

;(Q1+

3 )(βγ)

;(Q1−

3 )(βγ) are obtained in terms of the

temperature microvariables T(βγ)

(mn) as follows

(Q1+

2 )(βγ) = T(βγ)

(10) (°k(βγ)

20 ° hβ

2

k(βγ)

22 ) + T(βγ)

(20) (

°3hβ

2

k(βγ)

20 °

3h2

β

4 k(βγ)

22 )

(Q1−

2 )(βγ) = T(βγ)

(10) (°k(βγ)

20 +

hβ

2

k(βγ)

22 ) + T(βγ)

(20) (

3hβ

2

k(βγ)

20 °

3h2

β

4 k(βγ)

22 )

(Q1+

3 )(βγ) = T(βγ)

(01) (°k(βγ)

30 ° lγ

2

k(βγ)

33 ) + T(βγ)

(02) (

°3lγ

2

k(βγ)

30 ° 3l

2

γ

4 k(βγ)

33 )

(Q1−

3 )(βγ) = T(βγ)

(01) (°k(βγ)

30 +

lγ

2

k(βγ)

33 ) + T(βγ)

(02) (

3lγ

2

k(βγ)

30 ° 3l

2

γ

4 k(βγ)

33 )

The above four equations can be rewritten in matrix form as follows

∙ °Q1−

2

Q1+

2

 ̧(βγ)

=

2

6

6

4

k20 ° k22

hβ

2 °3hβk202 +

3h2

β

4 k22

°k20 ° k22

hβ

2 °3hβk202 °

3h2

β

4 k22

3

7

7

5

(βγ)

∙ T(10)

T(20)  ̧(βγ)

(80)

NASA/CR—2002-211910 27
and

∙ °Q1−

3

Q1+

3

 ̧(βγ)

=

2

6

4

k30 ° k33

lγ

2 °3lγk302 +

3l

2

γ

4 k33

°k30 ° k33

lγ

2 °3lγk302 ° 3l

2

γ

4 k33

3

7

5

(βγ)

∙ T(01)

T(02)  ̧(βγ)

(81)

Substituting for T(βγ)

(10) ; T(βγ)(20) and T(βγ)(01) ; T(βγ)(02) using Eqs. (78) and (79), the above equations become

∙ °Q1−

2

Q1+

2

 ̧(βγ)

=

2

6

6

4

2k22 ° 4k20

hβ

k22 ° 2k20

hβ

°k22 ° 2k20

hβ

°2k22 ° 4k20

hβ

3

7

7

5

(βγ)

∙ T1−

2

T1+

2

 ̧(βγ)

+

2

6

6

4

6k20

hβ

° 3k22

6k20

hβ

+ 3k22

3

7

7

5

(βγ)

T(βγ)

(00)

(82)

and

∙ °Q1−

3

Q1+

3

 ̧(βγ)

=

2

6

6

4

2k33 ° 4k30

lγ

k33 ° 4k30

lγ

°k33 ° 4k30

lγ

°2k33 ° 4k30

lγ

3

7

7

5

(βγ)

∙ T1−

3

T1+

3

 ̧(βγ)

+

2

6

6

4

6k30

lγ

° 3k33

6k30

lγ

° 3k33

3

7

7

5

(βγ)

T(βγ)

(00)

(83)

The last step is to develop a relationship between the surface-averaged temperatures and the

volume-averaged temperature T(βγ)

(00) . This relationship is obtained from the volume-averaged heat

conduction equation

Z lγ/2

−lγ/2

Z hβ/2

−hβ/2

Ã

@q(βγ)

2

@x1

(β)

2

+ @q(βγ)

3

@x1

(γ)

3

!

dx1

(β)

2 dx1

(γ)

3 = 0 (84)

where q

(βγ)

2 and q

(βγ)

3 are given by Eqs. (33) and (34). Since the thermal conductivity coefficients

k(βγ)

2 and k(βγ)3 are functions of the coordinates x1

(β)

2 and x1

(γ)

3 , the above equation can be developed

by using the chain rule as

lγ/2 Z

−lγ/2

hβ/2 Z

−hβ/2

(

@k(βγ)

2

@x1

(β)

2

@T(βγ)

@x1

(β)

2

+

@2T(βγ)

@x1

2 (β)

2

k(βγ)

2 +

@k(βγ)

3

@x1

(γ)

3

@T(βγ)

@x1

(γ)

3

+

@2T(βγ)

@x1

2 (γ)

3

k(βγ)

3 )dx1

(β)

2 dx1

(γ)

3 = 0 (85)

Upon determining the derivatives of temperature and thermal conductivity, respectively, the above

equation becomes

lγ/2 Z

−lγ/2

hβ/2 Z

−hβ/2

[k(βγ)

22 T(βγ)

(10) + 3T(βγ)(20) (2k(βγ)

22 x1

(β)

2 + k(βγ)20 ) +

k(βγ)

33 T(βγ)

(01) + 3T(βγ)(02) (2k(βγ)

33 x1

(γ)

3 + k(βγ)30 )]dx1

(β)

2 dx1

(γ)

3 = 0 (86)

Substituting for T(10)(βγ)

, T(20)(βγ); and T(01)(βγ); T(02)(βγ) using Eqs. (78) and (79), the above

equation becomes

NASA/CR—2002-211910 28
k(βγ)

22

hβ

[(T1+

2 )

(βγ) ° (T1−

2 )

(βγ)

] + 3k(βγ)

20

h2

β

[2(T1−

2 )

(βγ) + 2(T1+

2 )

(βγ) ° 4T(βγ)

(00) ] +

k(βγ)

33

lγ

[(T1+

3 )

(βγ) ° (T1−

3 )

(βγ)

] + 3k(βγ)

30

l2

γ

[2(T1−

3 )

(βγ) + 2(T1+

3 )

(βγ) ° 4T(βγ)

(00) ]=0 (87)

Therefore, the expression for T(βγ)

(00) can be obtained from the above equation in terms of (T1−

2 )(βγ)

;

(T1+

2 )(βγ)

; (T1−

3 )(βγ)

; (T1+

3 )(βγ) as follows

T(βγ)

(00) = [(T1−

2 )

(βγ)

(

6k(βγ)

20

h2

β

° k(βγ)

22

hβ

)+(T1+

2 )

(βγ)

(

k(βγ)

22

hβ

+

6k(βγ)

20

h2

β

) +

(T1−

3 )

(βγ)

(

6k(βγ)

30

l2

γ

° k(βγ)

33

lγ

)+(T1+

3 )

(βγ)

(

k(βγ)

33

lγ

+

6k(βγ)

30

l2

γ

)]=12(k(βγ)

20

h2

β

+

k(βγ)

30

l2

γ

) (88)

Upon substituting the above expression for T(βγ)

(00) into the expressions for (Q1+

2 )(βγ)

; (Q1−

2 )(βγ)

;

(Q1+

3 )(βγ)

; (Q1−

3 )(βγ) given by Eqs. (82) and (83), the local conductivity matrix κ(βγ) is obtained

which directly relates the surface-averaged heat fluxes to the corresponding surface-averaged temperatures. These expressions can be written in matrix form as follows

2

6

6

4

°Q1−

2

Q1+

2

°Q1−

3

Q1+

3

3

7

7

5

(βγ)

=

2

6

6

4

∙11 ∙12 ∙13 ∙14

∙21 ∙22 ∙23 ∙24

∙31 ∙32 ∙33 ∙34

∙41 ∙42 ∙43 ∙44

3

7

7

5

(βγ) 2

6

6

4

T1−

2

T1+

2

T1−

3

T1+

3

3

7

7

5

(βγ)

(89)

with the individual local conductivity matrix elements given below in terms of the subcell thermal

conductivity coefficients and geometry, upon setting  ̧ = 12(k(βγ)

20 =h2

β + k(βγ)

30 =l2

γ) for convenience.

∙(βγ)

11 = 2k(βγ)22 ° 4k(βγ)

20

hβ

+ (6k(βγ)

20

hβ

° 3k(βγ)

22 )(°k(βγ)

22

hβ

+

6k(βγ)

20

h2

β

)= ̧

∙(βγ)

12 = k(βγ)22 ° 2k(βγ)

20

hβ

+ (6k(βγ)

20

hβ

° 3k(βγ)

22 )(k(βγ)

22

hβ

+

6k(βγ)

20

h2

β

)= ̧

∙(βγ)

13 = (6k(βγ)

20

hβ

° 3k(βγ)

22 )(°k(βγ)

33

lγ

+

6k(βγ)

30

l2

γ

)= ̧

∙(βγ)

14 = (6k(βγ)

20

hβ

° 3k(βγ)

22 )(k(βγ)

33

lγ

+

6k(βγ)

30

l2

γ

)= ̧

∙(βγ)

21 = °k(βγ)22 ° 2k(βγ)

20

hβ

+ (6k(βγ)

20

hβ

+ 3k(βγ)

22 )(°k(βγ)

22

hβ

+

6k(βγ)

20

h2

β

)= ̧

∙(βγ)

22 = °2k(βγ)22 ° 4k(βγ)

20

hβ

+ (6k(βγ)

20

hβ

+ 3k(βγ)

22 )(k(βγ)

22

hβ

+

6k(βγ)

20

h2

β

)= ̧

∙(βγ)

23 = (6k(βγ)

20

hβ

+ 3k(βγ)

22 )(°k(βγ)

33

lγ

+

6k(βγ)

30

l2

γ

)= ̧

∙(βγ)

24 = (6k(βγ)

20

hβ

+ 3k(βγ)

22 )(k(βγ)

33

lγ

+

6k(βγ)

30

l2

γ

)= ̧

NASA/CR—2002-211910 29
∙(βγ)

31 = 2k(βγ)33 ° 4k(βγ)

30

lγ

+ (6k(βγ)

30

lγ

° 3k(βγ)

33 )(°k(βγ)

33

lγ

+

6k(βγ)

30

l2

γ

)= ̧

∙(βγ)

32 = k(βγ)33 ° 2k(βγ)

30

lγ

+ (6k(βγ)

30

lγ

° 3k(βγ)

33 )(k(βγ)

33

lγ

+

6k(βγ)

30

l2

γ

)= ̧

∙(βγ)

33 = (6k(βγ)

30

lγ

° 3k(βγ)

33 )(°k(βγ)

22

hβ

+

6k(βγ)

20

h2

β

)= ̧

∙(βγ)

34 = (6k(βγ)

30

lγ

° 3k(βγ)

33 )(k(βγ)

22

hβ

+

6k(βγ)

20

h2

β

)= ̧

∙(βγ)

41 = °k(βγ)33 ° 2k(βγ)

30

lγ

+ (6k(βγ)

30

lγ

+ 3k(βγ)

33 )(°k(βγ)

33

lγ

+

6k(βγ)

30

l2

γ

)= ̧

∙(βγ)

42 = °2k(βγ)33 ° 4k(βγ)

30

lγ

+ (6k(βγ)

30

lγ

+ 3k(βγ)

33 )(k(βγ)

33

lγ

+

6k(βγ)

30

l2

γ

)= ̧

∙(βγ)

43 = (6k(βγ)

30

lγ

+ 3k(βγ)

33 )(°k(βγ)

22

hβ

+

6k(βγ)

20

h2

β

)= ̧

∙(βγ)

44 = (6k(βγ)

30

lγ

+ 3k(βγ)

33 )(k(βγ)

22

hβ

+

6k(βγ)

20

h2

β

)= ̧

3.2 Constant Thermal Conductivity

When the coefficients associated with the local coordinates x1

(β)

2 and x1

(γ)

3 in the expressions for k(βγ)2

and k(βγ)

3 are set to zero, i.e., let k(βγ)22 = k(βγ)23 = k(βγ)32 = k(βγ)33 = 0 in Eqs. (76) and (77), constant

thermal conductivities are obtained, i.e., k(βγ)

2 = k(βγ)20 and k(βγ)3 = k(βγ)30 :

Replacing k(βγ)

20 with k(βγ)2 ; k(βγ)30 with k(βγ)3 and letting k(βγ)22 = k(βγ)33 = 0 in the local conductivity

matrix κ(βγ) with variable thermal conductivity in Eq. (89), the new local conductivity matrix with

constant thermal conductivities is obtained.

2

6

6

4

°Q1−

2

Q1+

2

°Q1−

3

Q1+

3

3

7

7

5

(βγ)

=

2

6

6

4

∙11 ∙12 ∙13 ∙14

∙21 ∙22 ∙23 ∙24

∙31 ∙32 ∙33 ∙34

∙41 ∙42 ∙43 ∙44

3

7

7

5

(βγ) 2

6

6

4

T1−

2

T1+

2

T1−

3

T1+

3

3

7

7

5

(βγ)

(90)

where

∙(βγ)

11 = ∙(βγ)22 = °2k(βγ)

2

hβ

(2 ° k(βγ)

2

21k(βγ)

2

)

∙(βγ)

12 = ∙(βγ)21 = °2k(βγ)

2

hβ

(1 ° k(βγ)

2

21k(βγ)

2

)

∙(βγ)

13 = ∙(βγ)14 = ∙(βγ)23 = ∙(βγ)24 = k(βγ)

2 k(βγ)3 hβ

1

k(βγ)

2 l2

γ

∙(βγ)

31 = ∙(βγ)32 = ∙(βγ)41 = ∙(βγ)42 = k(βγ)

2 k(βγ)3 lγ

1

k(βγ)

2 h2

β

∙(βγ)

33 = ∙(βγ)44 = °2k(βγ)

3

lγ

(2 ° k(βγ)

3

21k(βγ)

3

)

NASA/CR—2002-211910 30
∙(βγ)

34 = ∙(βγ)43 = °2k(βγ)

3

lγ

(1 ° k(βγ)

3

21k(βγ)

3

)

and where

1

k(βγ)

2 = h2

β(

k(βγ)

2

h2

β

+

k(βγ)

3

l2

γ

)=3

1

k(βγ)

3 = l

2

γ(

k(βγ)

2

h2

β

+

k(βγ)

3

l2

γ

)=3

This is the same local conductivity matrix as that obtained by Bansal (2002), which confirms

that the present development reduces correctly.

3.3 Global Conductivity Matrix Assembly

As mentioned above, the assembly of the global conductivity matrix is carried out by applying

heat flux continuity conditions, Eq. (72) and (73), at each interface along x2 and x3 directions,

respectively, together with the temperature continuity condition defined by Eqs. (74) and (75),

starting at the boundary with prescribed boundary conditions.

3.3.1 Temperature continuity conditions

From the temperature continuity, we know that the surface-averaged temperature (T1+

2 )(βγ) of the

(Ø; °) subcell is the same as (T1−

2 )(β+1,γ) of the (Ø+1; °) subcell defined in Eq. (74). The same holds

true in the x3 direction, Eq. (75). These conditions are applied Ø = 1; :::; Nβ°1 and ° = 1; :::; Nγ°1

times at the subcell interfaces in the x2 and x3 directions, respectively, in the heat flux continuity

conditions described below. Thus, altogether, the number of unknown surface-averaged interfacial

temperatures in the interior is (Nβ ° 1)Nγ + (Nγ ° 1)Nβ.

3.3.2 Heat flux continuity conditions

The continuity of interfacial surfaced-averaged heat flux along each interface requires that the sum

of the heat fluxes across the interfaces separating adjacent subcells be zero, i.e., Eq. (72) in the x2

direction and Eq. (73) in the x3 direction, respectively. Expressing these two equations in terms

of the surface-averaged temperatures by using Eq. (63) yields for the interface between (Ø; °) and

(Ø + 1; °) subcells,

∙(β,γ)

21 (T1−

2 )

(βγ) + ∙(β,γ)

22 (T1+

2 )

(βγ) + ∙(β,γ)

23 (T1−

3 )

(βγ) + ∙(β,γ)

24 (T1+

3 )

(βγ) + ∙(β+1,γ)

11 (T1−

2 )

(β+1,γ)

+∙(β+1,γ)

12 (T1+

2 )

(β+1,γ) + ∙(β+1,γ)

13 (T1−

3 )

(β+1,γ) + ∙(β+1,γ)

14 (T1+

3 )

(β+1,γ) = 0 (91)

while for the interface between (Ø; °) and (Ø; ° + 1) subcells we get,

∙(β,γ)

41 (T1−

2 )

(βγ) + ∙(β,γ)

42 (T1+

2 )

(βγ) + ∙(β,γ)

43 (T1−

3 )

(βγ) + ∙(β,γ)

44 (T1+

3 )

(βγ) + ∙(β,γ+1)

31 (T1−

2 )

(β,γ+1)

+∙(β,γ+1)

32 (T1+

2 )

(β,γ+1) + ∙(β,γ+1)

33 (T1−

3 )

(β,γ+1) + ∙(β,γ+1)

34 (T1+

3 )

(β,γ+1) = 0 (92)

Employing the temperature continuity conditions, Eqs. (74) and (75), the above equations become

∙(β,γ)

21 T1(β,γ)2 + (∙(β,γ)22 + ∙(β+1,γ)11 )T1(β+1,γ)2 + ∙(β+1,γ)12 T1(β+2,γ)2 + ∙(β,γ)23 T1(β,γ)3 (93)

+∙(β,γ)

24 T1(β,γ+1)3 + ∙(β+1,γ)13 T1(β+1,γ)3 + ∙(β+1,γ)14 T1(β+1,γ+1)3 = 0

NASA/CR—2002-211910 31
∙(β,γ)

41 T1(β,γ)2 + ∙(β,γ)42 T1(β+1,γ)2 + ∙(β,γ+1)31 T1(β,γ+1)2 + ∙(β,γ+1)32 T1(β+1,γ+1)2 + ∙(β,γ)43 T1(β,γ)3

+(∙(β,γ)

44 + ∙(β,γ+1)33 )T1(β,γ+1)3 + ∙(β,γ+1)34 T1(β,γ+2)3 = 0 (94)

Thus, these two equations provide us with (Nβ ° 1)Nγ + (Nγ ° 1)Nβ equations for the unknown

(Nβ ° 1)Nγ + (Nγ ° 1)Nβ common interfacial surface-averaged temperatures.

3.3.3 Boundary conditions

At external boundaries, the remaining 2(Nβ +Nγ) equations are obtained from the imposed boundary conditions specified either in terms of applied temperature or heat flux as given below

T1(β,1)

3 jxØ3=−l1/2 = 1

hβ

Z hβ/2

−hβ/2

T(β)

bot dx1

(β)

2 (95)

or

Q1(β,1)

3 jxØ3=−l1/2 = 1

hβ

Z hβ/2

−hβ/2

q

(β)

bot dx1

(β)

2 (96)

and

T1(β,Nγ)

3 jxØ3=lNγ /2 = 1

hβ

Z hβ/2

−hβ/2

T(β)

top dx1

(β)

2 (97)

or

Q1(β,Nγ)

3 jxØ3=lNγ /2 = 1

hβ

Z hβ/2

−hβ/2

q

(β)

top dx1

(β)

2 (98)

at (Ø; 1) and (Ø;Nγ) subcells. Similar reasoning holds for subcells (1; °); and (Nβ; °) where the

applied boundary conditions are given by

T1(1,γ)

2 jxØ2=−h1/2 = 1

lγ

Z lγ/2

−lγ/2

T(γ)

lef tdx1

(γ)

3 (99)

or

Q1(1,γ)

2 jxØ2=−h1/2 = 1

lγ

Z lγ/2

−lγ/2

q

(γ)

lef tdx1

(γ)

3 (100)

and

T1(Nβ,γ)

2 jxØ2=hNβ /2 = 1

lγ

Z lγ/2

−lγ/2

T(γ)

rightdx1

(γ)

3 (101)

or

Q1(Nβ,γ)

2 jxØ2=hNβ /2 = 1

lγ

Z lγ/2

−lγ/2

q

(γ)

rightdx1

(γ)

3 (102)

3.3.4 Governing equations

The general form of the final system of equations can be written symbolically as

Q = κ^T (103)

which can be decomposed for convenience as follows

∙ Q1 2

Q1 3

 ̧

=

∙ κ^11 κ^12

κ^21 κ^22  ̧ ∙ T12

T13

 ̧

(104)

NASA/CR—2002-211910 32
where the heat flux vectors Q1 2 and Q1 3 are defined as

Q1 2= [Q1 (1)

2 ; :::; Q1 (γ)

2 ; :::; Q1 (Nγ)

2 ]

Q1 3 = [Q1 (1)

3 ; :::; Q1(β)3 ; :::; Q1 (Nβ)

3 ] (105)

where

Q1 (γ)

2 = [Q1(1,γ)2 ; 0; :::; 0; Q1(Nβ+1,γ)

2 ]

Q1 (β)

3 = [Q1(β,1)3 ; 0; :::; 0; Q1(β,Nγ+1)

3 ]

The vectors T12 and T13 contain the unknown surface-averaged temperatures at the subcell

interfaces and the outer edges of the composite and are given by

T12 = [T1(1)

2 ; :::; T1(γ)

2 ; :::; T1(Nγ)

2 ]

T13 = [T1(1)

3 ; :::; T1(β)3 ; :::; T1(Nβ)

3 ] (106)

where

T1(γ)

2 = [T1(1,γ)2 ; :::; T1(Nβ+1,γ)

2 ]

T1(β)

3 = [T1(β,1)3 ; :::; T1(β,Nγ+1)

3 ]

The global thermal conductivity matrix κ^ is obtained by assembling the local conductivity

matrices as explained above. It contains information on the geometry and thermal conductivities

of the individual NβNγ subcells. The general structure of the global conductivity submatrices

κ^11;κ^12;κ^21;κ^22 is given below

κ^11 =

2

6

6

6

6

6

6

6

4

A(1)

#

β

000 0

0 ::::: 00 0

0 0A(γ)

#

β

0 0

000 ::::: 0

0 0 0 0A(Nγ)

#

β

3

7

7

7

7

7

7

7

5

κ^12 =

2

6

6

6

6

6

6

4

B(1,1)

βØ B(2,1)βØ ::::: B(Nβ,1)

βØ

B(1,2)

βØ B(2,2)βØ ::::: B(Nβ,2)

βØ

::::: ::::: ::::: :::::

::::: ::::: ::::: :::::

B(1,Nγ)

βØ B(2,Nγ)βØ ::::: B(Nβ,Nγ)

βØ

3

7

7

7

7

7

7

5

κ^21 =

2

6

6

6

6

6

6

4

B(1,1)

#γ B(1,2)#γ ::::: B(1,Nγ)

#γ

B(2,1)

#γ B(2,2)#γ ::::: B(2,Nγ)

#γ

::::: ::::: ::::: :::::

::::: ::::: ::::: :::::

B(Nβ,1)

#γ B(Nβ,2)#γ ::::: B(Nβ,Nγ)#γ

3

7

7

7

7

7

7

5

κ^22 =

2

6

6

6

6

6

6

4

A(1)

#γ 0 00 0

0 ::::: 00 0

0 0A(β)

#γ 0 0

000 ::::: 0

0 0 00A(Nβ)

#γ

3

7

7

7

7

7

7

5

The elements of the submatrices A(γ)

#

β ; A(β)

#γ ; B(β,γ)

βØ ; and B(β,γ)

γØ are given below

A(γ)

#

β =

2

6

6

6

6

6

6

6

6

4

∙(1,γ)

11 ∙(1,γ)12 0 0 0 00 0

:::: :::: :::: 0 0 00 0

0 0 ∙(β,γ)

11 ∙(β,γ)12 0 00 0

0 0 ∙(β,γ)

21 ∙(β,γ)22 + ∙(β+1,γ)11 ∙(β+1,γ)12 00 0

000 :::: :::: :::: :::: ::::

000 0 0 0 ∙

(Nβ,γ)

21 ∙

(Nβ,γ)

22

3

7

7

7

7

7

7

7

7

5

NASA/CR—2002-211910 33
where ° = 1; 2; :::::Nγ:

A(β)

#γ =

2

6

6

6

6

6

6

6

6

4

∙(β,1)

33 ∙(β,1)34 0 0 0 00 0

:::: :::: :::: 0 0 00 0

0 0 ∙(β,γ)

33 ∙(β,γ)34 0 00 0

0 0 ∙(β,γ)

43 ∙(β,γ)44 + ∙(β,γ+1)33 ∙(β,γ+1)34 00 0

000 :::: :::: :::: :::: ::::

000 0 0 0 ∙(β,Nγ)

43 ∙(β,Nγ)44

3

7

7

7

7

7

7

7

7

5

where Ø = 1; 2; :::Nβ:

B(β,γ)

#

β =

2

6

6

6

6

6

6

6

6

6

6

6

6

6

4

γ ! 1 2 :::: γthcol: (γ + 1)

thcol: :::: Nγ+1

β #

1 0 0 :::: 0 0 :::: 0

2 0 0 :::: 0 0 :::: 0

:::: :::: :::: :::: :::: :::: :::: ::::

βthrow 0 0 :::: ∙(β,γ)

13 ∙(β,γ)14 :::: 0

(β + 1)

throw 0 0 :::: ∙(β,γ)

23 ∙(β,γ)24 :::: 0

:::: :::: :::: :::: :::: :::: :::: ::::

Nβ+1 0 0 :::: 0 0 :::: 0

3

7

7

7

7

7

7

7

7

7

7

7

7

7

5

B(β,γ)

#γ =

2

6

6

6

6

6

6

6

6

6

6

6

6

6

4

Ø! 1 2 :::: βthcol: (β + 1)

thcol: :::: Nβ+1

γ #

1 0 0 :::: 0 0 :::: 0

2 0 0 :::: 0 0 :::: 0

:::: :::: :::: :::: :::: :::: :::: ::::

γthrow 0 0 :::: ∙(β,γ)

31 ∙(β,γ)32 :::: 0

(γ + 1)

th

row 0 0 :::: ∙(β,γ)

41 ∙(β,γ)42 :::: 0

:::: :::: :::: :::: :::: :::: :::: ::::

Nγ+1 0 0 :::: 0 0 :::: 0

3

7

7

7

7

7

7

7

7

7

7

7

7

7

5

3.4 Reformulated vs Unreformulated HOTFGM

The motivation for the above outlined reformulation of HOTFGM is the reduction in the number of

simultaneous equations whose solution determines the temperature field in the individual subcells.

The reduction is a direct result of enforcing the thermal continuity at the common subcell interfaces

directly. This is possible because the fundamental unknowns in the reformulation are the surfaceaveraged interfacial temperatures. In the original formulation, the size of the global conductivity

matrix is 5NβNγ £ 5NβNγ, while in the reformulation version we have (2NβNγ + Nβ + Nγ) £

(2NβNγ + Nβ + Nγ): Figure 9 compares the size of the global thermal conductivity matrices in

the two versions. For discretized volumes with the same number of subcells in the two directions,

Nβ = Nγ, the reduction can be as much as about 60%.

NASA/CR—2002-211910 34
Figure 9. Reduction in the size of the global conductivity matrix due to reformulation.

NASA/CR—2002-211910 35
4 Analytical Solution

In this section, an analytical solution methodology is developed for determining temperature distributions in homogeneous and layered plates subjected to steady-state thermal loading. The solutions

will be employed in validating both the user-friendly and reformulated versions of the higher-order

theory. First, we develop a solution for a single strip or plate subjected to thermal boundary

conditions, and then generalize it to layered plates consisting of an arbitrary number of strips.

4.1 Single Plate/Strip Case

Consider the problem of determining the temperature distribution in a rectangular plate whose

length is a and height is b in the x2 ° x3 plane, Figure 10, subjected to steady-state thermal

boundary conditions to be specified later. The thermal conductivity of the plateís material is

constant. Therefore, the governing equation for the temperature distribution within the plate in

this problem is the Laplace equation

"2T = 0 (107)

We solve this equation by the separation-of-variables method, seeking a solution of the form

T(x2; x3) = X(x2)Y (x3) (108)

Figure 10. Problem definition.

Substituting Eq. (108) into Eq. (107) and separating the variables gives

X!!

X = °Y !!Y = k2 (109)

so that

X!! ° k2X = 0 (110)

Y !!

+ k2Y = 0 (111)

and

X(x2) = f A + Bx2; k = 0

C cosh kx2 + D sinh kx2; k 6= 0 (112)

NASA/CR—2002-211910 36
Y (x3) = f E + F x3; k = 0

G cos kx3 + H sin kx3; k = 0 6 (113)

Because the Laplace equation "2T = 0 is linear, we can use superposition and combine the k = 0

and k = 0 6 product solutions as follows

T(x2; x3)=(A + Bx2)(E + F x3)+(C cosh kx2 + D sinh kx2)(G cos kx3 + H sin kx3) (114)

The coefficients A; :::::; H depend on the particular boundary conditions applied at the plateís edges.

We consider two cases below:

Case (i)

The first case produces an internal temperature field that has no axes of symmetry when a =6 b,

and when a = b the temperature field is symmetric about the diagonal connecting the upper left

and the lower right corners. The corresponding boundary conditions are

T(0; x3) = 100; T(a; x3) = 200

T(x2; 0) = 200; T(x2; b) = 100

Application of the boundary conditions at the bottom face x3 = 0; where T = 200; yields

(A + Bx2)E + (C cosh kx2 + D sinh kx2)G = 200

from which we obtain

AE = 200; B = 0; G = 0

Therefore Eq. (114) becomes

T(x2; x3) = 200 + Jx3 + (C!

cosh kx2 + D!sinh kx2) sin kx3 (115)

where AF has been combined into J; CH into C!

, and DH into D!for convenience.

Next, application of the boundary condition at the top face x3 = b; where T = 100; yields

200 + Jb + (C!cosh kx2 + D!sinh kx2) sin kb = 100

from which we obtain

J = °100=b; and sin kb = 0

Hence,

k = mo=b (m = 1; 2; :::::)

Therefore, Eq.(115) becomes

T(x2; x3) = 200 ° 100x3

b +X∞

n=1

(Cm cosh mox2

b + Dm sinh mox2b ) sin mox3b (116)

Application of the boundary condition at the left face x2 = 0; where T = 100; yields

200 ° 100x3

b + X∞

m=1

Cm sin mox3

b = 100 (117)

Solving for the unknown coefficients Cm using the orthogonality properties of sin mox3

b ; we obtain

Cm = 2

b

Z b

0

(

100x3

b ° 100) sin mox3b dx3 (118)

NASA/CR—2002-211910 37
or

Cm = °200

mo (119)

Finally, application of the boundary condition at the right face x2 = a; where T = 200; yields

200 ° 100x3

b + X∞

m=1

(°200

mo

cosh moa

b + Dm sinh moab ) sin mox3b = 200 (120)

from which the solution for Dm can be obtained in the form

Dm = 2

b sinh(moa

b )

Z b

0

(

100x3

b sin mox3b +

200

mo

cosh moa

b sin2 mox3b )dx3 (121)

or

Dm =

200[cosh(moa

b ) ° cos mo]

mo sinh(moa

b ) (122)

Thus, the solution to Eq. (107) is given by Eq. (116), with the coefficients Cm and Dm computed

using Eqs. (119) and (122).

Case (ii)

If the boundary conditions are changed to produce temperature distribution which is symmetric

about the x3 = b=2 line, i.e.,

T(0; x3) = 100; T(a; x3) = 200

T(x2; 0) = 100; T(x2; b) = 100

then we obtain a simpler form of solution. First, applying the boundary condition T(x2; 0) = 100;

we get

T(x2; x3) = 100 + Jx3 + (C!

cosh kx2 + D!sinh kx2) sin kx3 (123)

where AE = 100; B = 0; G = 0; and CH has been combined into C!

;DH into D!; and AF into J:

Next, application of the boundary condition T(x2; b) = 100 in the above equation produces

T(x2; b) = 100 + Jx3 + (C!

cosh kx2 + D!sinh kx2) sin kb = 100 (124)

from which we obtain

J = 0; sin kb = 0

Therefore, the solution becomes

T(x2; x3) = 100 +X∞

n=1

(Cm cosh mox2

b + Dm sinh mox2b ) sin mox3b (125)

After application of the boundary conditions T(0; x3) = 100 and T(a; x3) = 200; the unknown

coefficients Cm; Dm are obtained, by using orthogonality properties of sin moy

b ; in the form

Cm = 0 (126)

Dm = 2

sinh moa

b

Z b

0

100 sin mox3

b dx3 = 400

mo sinh moa

b

(m = 1; 3; 5; :::) (127)

Thus, the solution in this case is given by Eqs. (125), (126), and (127).

NASA/CR—2002-211910 38
4.2 Nβ-strip Case

In the Nβ-strip problem, the rectangular plate, which has length a and height b, is divided into Nβ

strips along the x2 direction. Therefore, every strip has height b and width aβ; where Ø = 1; :::::; Nβ;

Figure 11. Since the boundary conditions on the top and bottom surface are given in terms of

constant temperature, the analytical solution can be developed first from Eqs. (116) or (125) of

the previous section depending on the applied boundary conditions.

Figure 11. Nβ-strip case problem definition.

Case (i)

Assume the boundary conditions on the top and bottom faces are Ttop = 200 and Tbot = 100:

Then the temperature field in the arbitrary Øth strip is given by Eq. (116) as follows,

T(β)

(x(β)

2 ; x3) = 200 ° 100x3

b + X∞

m=1

(C(β) m cosh mox(β)

2

b + D(β) m sinh mox(β)

2

b ) sin mox3b

However, it should be noted that this solution is obtained in the coordinate system with the

origin located at the left-bottom corner. In order to obtain the solution in the coordinate system

with the origin in the plateís center, we use the following coordinate transformation, Figure 11.

x1

(β)

2 = x(β)2 ° aβ=2 ) x(β)2 = 1x(β)2 + aβ=2

x13 = x3 ° b=2 ) x3 = 1x3 + b=2

Therefore, the solution of the Laplaceís equation in the translated coordinate system becomes

T(β)

(1x(β)

2 ; x13) = X∞

m=1

[C(β) m cosh mo(1x(β)

2 + aβ=2)

b + D(β) m sinh mo(1x(β)

2 + aβ=2)

b ] £

sin mo(1x3 + b=2)

b + 200 ° 100(1x3 + b=2)b (128)

Case (ii)

If the boundary conditions are changed to produce temperature distribution which is symmetric

about the x3 = b=2 line, i.e.,

Ttop = Tbottom = 100

NASA/CR—2002-211910 39
the temperature field in the arbitrary Øth strip is given by Eq. (125) as follows

T(β)

(x(β)

2 ; x3) = 100 + X∞

m=1

(C(β) m cosh mox(β)

2

b + D(β) m sinh mox(β)

2

b ) sin mox3b

Upon transforming the above solution from the coordinate system with the origin located at the

left-bottom corner to that with the origin in the plateís center, the solution of the Laplaceís equation

with the symmetry boundary condition becomes,

T(β)

(1x(β)

2 ; x13) = 100 + X∞

m=1

[C(β) m cosh mo(1x(β)

2 + aβ=2)

b + D(β) m sinh mo(1x(β)

2 + aβ=2)

b ] £

sin mo(1x3 + b=2)

b (129)

In both cases, the solution for the unknown coefficients C(β) m ; D(β) m is obtained by constructing

a system of equations which satisfies the heat flux and temperature continuity conditions at the

interfaces separating the individual layers and the remaining boundary conditions. In order to

make this as efficient as possible for large number of layers, we formulate the solution procedure

using the local/global conductivity matrix approach.

In the local/global conductivity matrix approach, the unknown coefficients C(β) m and D(β) m are

the key elements in constructing the local conductivity matrix. They are expressed and calculated

in terms of the harmonics of common interfacial temperatures between adjacent strips and subsequently related to the corresponding interfacial heat flux harmonics through the local conductivity

matrix. The detailed derivation is outlined in the following section.

4.2.1 Construction of the local conductivity matrix

In order to construct the local conductivity matrix, we start with the heat flux constitutive equation

in the x2 direction

q

(β)

2 = °kβ

@T

@x1

(β)

2

; Ø = 1; :::::; Nβ (130)

Substituting the expression for the temperature field given by Eq. (128) or (129), the above equation

becomes

q

(β)

2 = °kβ

X∞

m=1

mo

b

"

C(β) m sinh mo(1x(β)

2 + aβ=2)

b + D(β) m cosh mo(1x(β)

2 + aβ=2)

b

#

sin mo(1x3 + b=2)

b

(131)

The determination of the coefficients C(β) m and D(β) m for case (i) and case (ii) is outlined below.

Case (i)

In order to express the coefficients C(β) m and D(β) m in terms of the surface temperatures on

the vertical boundaries of the Øth strip, Eq. (128) is evaluated at (aβ=2; x13) and (°aβ=2; x13);

respectively, yielding

T(β)

(aβ=2; x13) = 200 ° 100(1x3 + b=2)

b + X∞

m=1

[C(β) m cosh moaβ

b + D(β) m sinh moaβb ] sin mo(1x3 + b=2)b

(132)

NASA/CR—2002-211910 40
or

[T(β)(aβ=2; x13)°200+ 100(1x3 + b=2)

b ]m = X∞

m=1

[C(β) m cosh moaβ

b +D(β) m sinh moaβb ] sin mo(1x3 + b=2)b

(133)

and

T(β)

(°aβ=2; x13) = 200 ° 100(1x3 + b=2)

b + X∞

m=1

C(β) m sin mo(1x3 + b=2)

b (134)

or

[T(β)(°aβ=2; x13) ° 200 + 100(1x3 + b=2)

b ]m = X∞

m=1

C(β) m sin mo(1x3 + b=2)

b (135)

The corresponding interfacial temperature harmonics at the left (°aβ=2; x13) and right (aβ=2; x13)

boundary for an arbitrary Øth strip can be obtained from Eqs. 133) and (135), respectively, using

orthogonality relations

(T +

2 )

(β) m =

b/2 Z

−b/2

∙

T(β)

(aβ=2; x13) ° 200 + 100(1x3 + b=2)

b

 ̧

sin mo(1x3 + b=2)

b dx13

=

b/2 Z

−b/2

X∞

m=1

[C(β) m cosh moaβ

b + D(β) m sinh moaβb ] sin2 mo(1x3 + b=2)b dx13

and

(T −

2 )

(β) m =

b/2 Z

−b/2

∙

T(β)

(°aβ=2; x13) ° 200 + 100(1x3 + b=2)

b

 ̧

sin mo(1x3 + b=2)

b dx13

=

b/2 Z

−b/2

X∞

m=1

C(β) m sin2 mo(1x3 + b=2)

b dx13

Performing the required integrations, the above two equations become

(T +

2 )

(β) m = b

2

(C(β) m cosh moaβ

b + D(β) m sinh moaβb ) (136)

(T −

2 )

(β) m = b

2

C(β) m (137)

Upon writing the above two equations in matrix form, we obtain

∙ T −

2

T +

2

 ̧(β)

m

=

2

6

4

b

2 0

b

2

cosh moaβ

b

b

2

sinh moaβ

b

3

7

5

∙ Cm

Dm

 ̧(β)

(138)

Inverting the above system, the coefficients C(β) m and D(β) m are expressed in terms of the interfacial

temperature harmonics as follows

NASA/CR—2002-211910 41
∙ Cm

Dm

 ̧(β)

=

2

6

4

b

2 0

b

2

cosh moaβ

b

b

2

sinh moaβ

b

3

7

5

−1

∙ T −

2

T +

2

 ̧(β)

m

(139)

Similarly, evaluating the heat flux at the left (°aβ=2; x13) and right (aβ=2; x13) boundary for the Øth

strip, respectively, yields

q

(β)

2 (aβ=2; x13) = °kβ

X∞

m=1

(

mo

b C(β) m sinh moaβb + D(β) m cosh moaβb ) sin mo(1x3 + b=2)b

°q

(β)

2 (°aβ=2; x13) = kβ

X∞

m=1

mo

b D(β) m sin mo(1x3 + b=2)b

The interfacial heat flux harmonics for an arbitrary Øth strip are obtained from the above two

equations by using orthogonality relations,

(Q+

2 )

(β) m =

b/2 Z

−b/2

q

(β)

2 (aβ=2; x13) sin mo(1x3 + b=2)

b dx13 (140)

=

b/2 Z

−b/2

°kβ

X∞

m=1

(

mo

b C(β) m sinh moaβb + D(β) m cosh moaβb ) sin2 mo(1x3 + b=2)b dx13

°(Q−

2 )

(β) m =

b/2 Z

−b/2

q

(β)

2 (°aβ=2; x13) sin mo(1x3 + b=2)

b dx13

=

b/2 Z

−b/2

°kβ

X∞

m=1

mo

b D(β) m sin2 mo(1x3 + b=2)b dx13 (141)

The above two equations can be expressed in matrix form as follows

∙ °Q−

2

Q+

2

 ̧(β)

m

=

2

6

4

0 °mokβ

2

°mokβ

2

sinh moaβ

b °mokβ2

cosh moaβ

b

3

7

5

∙ Cm

Dm

 ̧(β)

(142)

Case (ii)

Similar to case (i), evaluating the surface temperatures on the vertical boundaries of the Øth

strip i.e., Eq. (129) at (aβ=2; x13) and (°aβ=2; x13); respectively, yields

T(β)

(aβ=2; x13) = 100 + X∞

m=1

[C(β) m cosh moa

b + D(β) m sinh moaβb ] sin mo(1x3 + b=2)b (143)

or

[T(β)(aβ=2; x13) ° 100]m = X∞

m=1

[C(β) m cosh moaβ

b + D(β) m sinh moaβb ] sin mo(1x3 + b=2)b (144)

NASA/CR—2002-211910 42
and

T(β)

(°aβ=2; x13) = 100 + X∞

m=1

C(β) m sin mo(1x3 + b=2)

b (145)

or

[T(β)(°aβ=2; x13) ° 100]m = X∞

m=1

C(β) m sin mo(1x3 + b=2)

b (146)

The corresponding interfacial temperature harmonics at the left (°aβ=2; x13) and right (aβ=2; x13)

boundary for an arbitrary Øth strip can be obtained from Eqs. (144) and (146), respectively, using

orthogonality relations

(T +

2 )

(β) m =

b/2 Z

−b/2

h

T(β)

(aβ=2; x13) ° 100isin mo(1x3 + b=2)

b dx13

=

b/2 Z

−b/2

X∞

m=1

[C(β) m cosh moaβ

b + D(β) m sinh moaβb ] sin2 mo(1x3 + b=2)b dx13

and

(T −

2 )

(β) m =

b/2 Z

−b/2

h

T(β)

(°aβ=2; x13) ° 100isin mo(1x3 + b=2)

b dx13

=

b/2 Z

−b/2

X∞

m=1

C(β) m sin2 mo(1x3 + b=2)

b dx13

Upon performing the above integrations, these two equations simplify to

(T +

2 )

(β) m = b

2

(C(β) m cosh moaβ

b + D(β) m sinh moaβb ) (147)

(T −

2 )

(β) m = b

2

C(β) m (148)

Writing the above two equations in matrix form, we obtain

∙ T −

2

T +

2

 ̧(β)

m

=

2

6

4

b

2 0

b

2

cosh moaβ

b

b

2

sinh moaβ

b

3

7

5

∙ Cm

Dm

 ̧(β)

(149)

which has the same form as Eq. (138) obtained for case (i) with the first kind of boundary

conditions. Further, since the expression for the heat flux, Eq. (131), for both boundary conditions

does not change, the relationship between the coefficients C(β) m and D(β) m and the interfacial heat

flux harmonics, i.e., Eq. (142), is also the same. However, the specific values of the coefficients C(β) m ,

D(β) m and the corresponding interfacial temperatures (T −

2 )

(β) m ;(T +

2 )

(β) m , heat fluxes (Q−

2 )

(β) m ;(Q+

2 )

(β) m

are not the same due to the different expressions for the solution of Laplaceís equation with different

kinds of boundary conditions.

NASA/CR—2002-211910 43
Nevertheless, in both cases, upon substituting the expression for the coefficients C(β) m and D(β) m

given by Eq. (139) into Eq. (142), carrying out the matrix multiplication and simplifying, we obtain

the same relationship between the interfacial heat flux and interfacial temperature harmonics, i.e.,

the local conductivity matrix,

∙ °Q−

2

Q+

2

 ̧(β)

m

=

∙ ∙11 ∙12

∙21 ∙22  ̧(β)

m

∙ T −

2

T +

2

 ̧(β)

m

(150)

where

(∙11)

(β) m =

mokβ cosh moaβ

b

b sinh moaβ

b

(∙12)

(β) m = ° mokβ

b sinh moaβ

b

(∙21)

(β) m = °mokβ

b sinh moaβb +

kβmo cosh2 moaβ

b

b sinh moaβ

b

(∙22)

(β) m = °

mokβ cosh moaβ

b

b sinh moaβ

b

4.2.2 Global conductivity matrix assembly

Imposition of the continuity of interfacial heat flux and interfacial temperature along common

interfaces between adjacent strips, together with the external boundary conditions, gives rise to a

system of equations in the unknown interfacial temperature harmonics. The continuity of interfacial

heat flux is obtained by requiring the sum of the heat flux acting on the interface separating Ø and

Ø + 1 strip to be zero,

(Q+

2 )

(β) m + (Q−

2 )

(β+1) m = 0 Ø = 1; ::::::Nβ ° 1 (151)

whereas the continuity of interfacial temperature is given by

(T +

2 )

(β) m = (T −

2 )

(β+1) m = T β

2 Ø = 1; ::::Nβ ° 1 (152)

The system equations is constructed by using Eq. (151) at each interface along the x2 direction, starting with the left boundary where the boundary condition is prescribed, together with

application of the temperature continuity condition defined by Eq. (152), and proceeding to the

right boundary. Following this procedure, we obtain the following system of equations for the mth

harmonic of the interfacial temperatures

left boundary: ∙1

11Tlef t

2 + ∙1

12T12 = °Qlef t

2

1st interface: ∙1

21Tlef t

2 + (∙1

22 + ∙211)T12 + ∙212T22 = 0

:

Øth interface: ∙β−1

21 T β−22 + (∙β−122 + ∙β11)T β−12 + ∙β12T β2 = 0

Ø + 1th interface: ∙β

21T β−12 + (∙β22 + ∙β+111 )T β2 + ∙β+112 T β+12 = 0

:

Nβ ° 1th interface: ∙

Nβ−1

21 T Nβ−22 + (∙

Nβ

11 + ∙

Nβ−1

22 )T Nβ−12 + ∙

Nβ

12 Tright

2 = 0

right boundary: ∙

Nβ

21 T Nβ−1

2 + ∙

Nβ

22 Tright

2 = Qright2

NASA/CR—2002-211910 44
where Nβ is the number of strips, and the harmonic number m has been suppressed. The above

equations can be expressed in matrix form for the mth harmonic as follows

2

6

6

6

6

6

6

6

6

6

6

6

4

∙1

11 ∙112 : 0 0 : 0 0

∙1

21 ∙122 + ∙211 : 0 0 : 0 0

0 0 : 0 0 : 0 0

0 0 : ∙β−1

22 + ∙β11 ∙β12 : 0 0

0 0 : ∙β

21 ∙β22 + ∙β+111 : 0 0

0 0 : 0 0 : 0 0

0 0 : 0 0 : ∙

Nβ−1

22 + ∙

Nβ

11 ∙

Nβ

12

0 0 : 0 0 : ∙

Nβ

21 ∙

Nβ

22

3

7

7

7

7

7

7

7

7

7

7

7

5

m

2

6

6

6

6

6

6

6

6

6

6

6

4

Tlef t

2

T1

2

:

Tk−1

2

Tk

2

:

T Nβ−1

2

Tright

2

3

7

7

7

7

7

7

7

7

7

7

7

5

m

=

2

6

6

6

6

6

6

6

6

6

6

6

4

°Qlef t

2

0

:

:

:

:

0

Qright

2

3

7

7

7

7

7

7

7

7

7

7

7

5

m

The solution of the above system of equations for each interfacial temperature harmonic makes

possible the determination of the coefficients C(β) m and D(β) m using Eq. (139). These, in turn,

determine the temperature field in each strip of the plate using either Eq. (128) or (129).

5 Mesh Sensitivity, Validation and Application Studies

The user-friendly formulation of HOTFGM, developed in Section 2, dispenses with the two-level

volume discretization employed in the original HOTFGM formulation wherein generic cells are further discretized into subcells. In the user-friendly formulation the graded microstructure is built

up directly using subcells which play the role of fundamental subvolumes. This simplification in

the volume discretization, which also simplified the derivation of the interfacial heat flux continuity

conditions, provided the basis for the efficient reformulation outlined in Section 3. Therefore, the

present section begins with the validation of the user-friendly formulation using the analytical solutions developed in Section 4. This is accomplished in two steps. First, we examine the convergence

of the temperature field to the analytical solution as a function of the mesh discretization for the

problem of a homogeneous plate subjected to the thermal boundary conditions discussed in Section 4. Once the convergence characteristics of the solution based on the user-friendly formulation

are established, we then consider the problem of a piece-wise uniformly graded plate for which

the analytical solution has also been derived in Section 4. This problem is subsequently re-solved

using the reformulation of HOTFGM with variable thermal conductivity in order to demonstrate

the reduction in the graded microstructureís volume dicretization obtained by assuming linearly

varying material properties. In the last example, a graded thermal barrier coating subjected to a

through-thickness thermal gradient is considered, and the differences in the solutions based on the

actual and homogenized microstructures are discussed.

5.1 User-Friendly Formulation

5.1.1 Mesh sensitivity study: homogeneous plate

Let us consider a square plate, whose sides are unit length (i.e., a = b = 1), Figure 12, subjected

to the following two sets of boundary conditions:

(1) : Tlef t = 100◦C; Tright = 200◦C; Ttop = 100◦C; Tbottom = 100◦C

(2) : Tlef t = 100◦C; Tright = 200◦C; Ttop = 100◦C; Tbottom = 200◦C

NASA/CR—2002-211910 45
Figure 12. Problem definition and mesh discretization for a homogeneous plate.

NASA/CR—2002-211910 46
The first set produces a temperature field that is symmetric about a horizontal line halfway across

the plateís height, whereas the second set produces a temperature field that is symmetric about

a line connecting the upper left and lower right corners of the plate. Such temperature fields are

sufficiently complicated to test the higher-order theoryís predictive capability. The plateís thermal

conductivity k has a constant value, i.e., k = 25 W=m °◦ C. Since the plate is homogeneous,

however, the temperature field does not depend on the conductivity under steady-state boundary

conditions. This provided an additional check on the user-friendly formulationís correctness as

verified during preliminary numerical tests.

The convergence behavior of the user-friendly formulation is investigated as a function of volume

discretization by dividing the plate into 4 £ 4; 12 £ 12; 32 £ 32 uniformly-spaced subcells, Figure

12. The corresponding contour plots of the temperature distributions for the employed volume

discretizations are compared with the analytical solution in Figure 13 for the first set of applied

thermal boundary conditions. As is clearly seen in the contour plots, the temperature distributions obtained from the user-friendly formulation are symmetric with respect to the horizontal line

halfway across the plateís height irrespective of the volume discretization, and gradually approach

the analytical solution with increasing number of subcells. In the case of the 4 £ 4 subcell mesh,

we obtain only an approximate temperature field which, nevertheless, exhibits the basic characteristics of the analytical solution. Increasing the volume discretization to 12 £ 12 subcells generates

a temperature field which appears quite close to the analytical solution except at the upper and

lower right corners where the boundary temperature distributions change suddenly, producing discontinuities at these corners. Increasing the mesh to 32 £ 32 subcells eliminates these local corner

disturbances, producing smooth curve patterns in the temperature field which are visually the same

as those obtained from the analytical solution.

A more quantitive picture of the user-friendly formulationís convergence behavior with mesh

refinement is obtained by plotting the temperature distributions along specific cross sections. Three

horizontal cross sections and three vertical cross sections have been chosen as illustrated by the

arrows in Figure 12. The horizontal cross sections are situated at the elevations x3 = 0:0; 0:25; 0:5;

and the vertical cross sections are situated at x2 = 0:5; 0:75; 1:0. The cross sections include the lower

horizontal boundary defined by x3 = 0:0 and the right vertical boundary defined by x2 = 1:0. The

interior cross sections have been chosen such that they run along subcell interfaces. For example,

the cross-section situated at x3 = 0:25 runs along the horizontal interfaces that separate the sets of

subcells (Ø1) and (Ø2) in the case of the 4£4 subcell mesh. This same cross section runs along the

horizontal interfaces that separate the sets of subcells (Ø3) and (Ø4) for the 12 £ 12 subcell mesh.

As we know from the user-friendly formulation, the temperature distributions along the different interior cross sections can be obtained by using microvariables associated with the subcell on

either side of the common interface. Furthermore, recall that the thermal and heat flux continuity

conditions were applied in a surface-average sense. Therefore, the temperature distributions along

a given cross section separating two sets of subcells, obtained using two different sets of microvariables associated with either set of subcells, will generally be different. However, these distributions

should approach the actual temperature distributions with increasing mesh refinement. Furthermore, the surface-averaged interfacial temperatures, which exhibit piece-wise uniform distributions,

should approach the actual temperature distributions in a step-wise manner with increasing mesh

refinement.

In order to demonstrate the above behavior, three types of temperature distributions along

the chosen subcell interfaces have been generated for the interior cross sections. In the following

horizontal cross-section plots, we define the common interfacial surface-averaged temperature as T; 1

the interfacial temperature obtained using the microvariables of the subcells below the interfaces

as T +; and the interfacial temperature obtained sing the microvariables of the subcells above the

NASA/CR—2002-211910 47
Figure 13. Temperature fields for the three mesh discretizations and first set of boundary

conditions. Comparison with the analytical solution.

NASA/CR—2002-211910 48
interfaces as T −: Similar notation is used for the vertical cross-section plots. In these cases, T +

indicates the interfacial temperature obtained using the microvariables of the subcells lying to

the left of the subcell interfaces and T − denotes the interfacial temperature obtained using the

microvariables of the subcells lying to the right of the subcell interfaces. In the case of the boundary

subcells, only two curves are given since there is only one set of subcells forming the boundary.

As before, the designation T1 is used to denote the surface-averaged temperature distribution along

both the horizontal and vertical boundaries. The second curve is obtained using the properties of

the boundary subcells: for the lower horizontal boundary it is denoted by T − and for the right

vertical boundary it is denoted by T + for consistency with the established notation.

The temperature distributions along the interior horizontal cross sections at the two elevations

x3 = 0:25 and x3 = 0:5 for the first set of boundary conditions are shown in Figure 14. For the

4 £ 4 subcell mesh, the temperature distributions calculated on the different sides of the subcell

interfaces in the x3 = 0:25 cross section shown in the left column of Figure 14 are quite different,

with discontinuities or jumps evident as the subcell vertical interfaces are traversed along the x2

direction. The differences, however, gradually disappear with increasing mesh refinement, including

the above-mentioned jumps. In the case of the 12£12 subcell mesh the differences are already quite

small and barely visually discernible, while in the case of the 32 £ 32 subcell mesh the differences

cannot be identified visually. These distributions are graphically indistinguishable from the actual

distributions obtained from the analytical solution (not shown). The surface-averaged temperature

distributions also gradually converge to the actual temperature distributions in a step-wise manner.

The corresponding distributions in the x3 = 0:5 cross section are shown in the right column

of Figure 14. The same observations with regard to the solutionís convergence with increasing

mesh refinement hold in this cross section as in the preceding one, with one important difference

however. Interestingly, we only see two curves in this cross section, namely the surface-averaged

temperature distribution T1 and the two T + and T − distributions superposed on each other. This

is due to the symmetric boundary conditions employed in this example, which in turn produce a

temperature field that is symmetric with respect to the x3 = 0:5 cross section. This, in turn, ensures

that the subcell microvariables on either side of the line of symmetry are the same, producing the

same interfacial temperature distributions. This behavior is clearly captured by the user-friendly

formulation, providing further evidence of its correctness.

Figure 15 shows the temperature distributions along the vertical cross sections x2 = 0:5 and

x2 = 0:75 as a function of the volume discretization for the first set of boundary conditions. In

this case, the temperature distributions are not symmetric about any of the two cross sections,

and consequently three distinct sets of interfacial temperature distributions, T; T 1 + and T − are

present. Symmetry exists only about the x3 = 0:5 cross section, and this is reflected in each of the

individual distributions for each cross section (since the temperature distributions are given as a

function of the x3 coordinate) and for each volume discretization. As in the case of the horizontal

cross sections, the difference between the T + and T − distributions vanish with increasing mesh

refinement (including the jumps at the subcell horizontal interfaces), while the T1 distribution

approaches the actual distribution in a step-wise manner.

The temperature distributions at the two external boundaries x3 = 0 and x2 = 1 are given in

Figure 16 for the three volume discretizations and the first set of boundary conditions. From these

figures, we conclude that the T − temperature distributions for the lower horizontal boundary, shown

in the left column of Figure 16, and the T + temperature distributions for the right vertical boundary,

shown in the right column of Figure 16, gradually approach the applied boundary conditions.

However, even with 32£32 subcell mesh, they exhibit some fluctuations at the lower and upper right

corners of the plate. That is because the applied temperature distributions possess discontinuities

at these corners, requiring further mesh discretization in the vicinity of these points.

NASA/CR—2002-211910 49
Figure 14. Convergence of temperature distributions along the cross-sections x3 = 0:25; 0:5.

NASA/CR—2002-211910 50
Figure 15. Convergences of the temperature distributions along the cross-sections x2 = 0:5; 0:75.

NASA/CR—2002-211910 51
Figure 16. Convergence of the temperature distributions along the edges x3 = 0 and x2 = 1.

NASA/CR—2002-211910 52
The corresponding temperature contour plots and cross section distributions for the second set

of boundary conditions are given in Figures 17 through 20. The full field contour distributions

shown in Figure 17 illustrate temperature field symmetry with respect to the cross section passing

through the upper left and lower right corners of the plate for all three volume discretizations.

As in the preceding case, the results obtained from the user-friendly formulation of the higherorder theory converge to the actual temperature distribution obtained from the analytical solution

with increasing mesh refinement. The full field temperature distribution obtained from the 32 £ 32

subcell mesh is visually indistinguishable from the distribution generated by the analytical solution.

Similar observations with regard to the convergence behavior with increasing mesh refinement hold

for the cross section temperature distributions along interfaces separating the interior subcells,

presented in Figures 18 and 19, and along external boundaries presented in Figure 20. In this case,

the temperature distributions are not symmetric about the chosen cross sections as expected, and

therefore three distinct sets of interfacial distributions are evident in the interior, which converge

with increasing mesh refinement at the same rate as before.

In summary, as observed from the above figures, whatever interfaces are considered, for either

set of the boundary conditions (Figures 13-15 or Figures 17-20), the differences between the temperature distributions based on the microvariables belonging to adjacent subcells decrease with

increasing number of subcells, while the surface-averaged temperature distributions approach the

actual distributions in a step-wise manner, with the step increments becoming increasingly smaller.

In other words, the greater the mesh refinement, the better the point-wise convergence or accuracy

to the actual distribution.

5.1.2 Validation: discretely graded plate

In this section, the results from the user-friendly formulation for a plate with discretely graded

thermal conductivity in a piecewise uniform manner are compared with the corresponding analytical

solution outlined in Section 4 for the Nβ-strip plate. The effect of thermal conductivity variation

on the temperature field is also investigated.

We consider the same boundary value problem discussed in the preceding section but limit the

investigation to the first type of boundary conditions, i.e., the plate is subjected to a temperature of

200◦C on the right boundary and a temperature of 100◦C on the remaining boundaries, Figure 21.

The plate is divided into vertical strips and the thermal conductivity k (assuming isotropic thermal

properties) is assigned a different constant value in each strip, producing a piecewise uniformly

graded plate. The overall variation of the thermal conductivity is given by exponential functions

in the x2 coordinate. Two different thermal conductivity variations are employed. One is an

exponentially increasing function of x2

k = 0:5e5x2

while the other is an exponentially decreasing function of x2

k = 60e−5x2

These continuous functions are employed to produce piecewise uniform variations as follows.

The plate is divided into 25 vertical strips of equal width and the above functions are used to

calculate the thermal conductivity at the center of each strip, which is then assumed to remain

constant in that strip, as shown in the two lower graphs of Figure 21. The analytical solutions

for this problem are then obtained based on the Nβ-strip solution outlined in Section 4 using

a 200 harmonic representation of the temperature field in each of the 25 strips. The corresponding user-friendly formulation solution is obtained by dividing the plate into 25£25 subcells of equal

NASA/CR—2002-211910 53
Figure 17. Temperature fields for the three mesh discretizations and second set of boundary

conditions. Comparison with the analytical solution.

NASA/CR—2002-211910 54
Figure 18. Convergence of the temperature distributions along the cross-sections x2 = 0:5; 0:75.

NASA/CR—2002-211910 55
Figure 19. Convergence of the temperature distribution along the cross-sections x3 = 0:25; 0:5.

NASA/CR—2002-211910 56
Figure 20. Convergence of the temperature distributions along the edges x3 = 0 and x2 = 1.

NASA/CR—2002-211910 57
Figure 21. Problem definition and mesh discretization for a layered plate with piece-wise

uniformly graded thermal conductivity.

NASA/CR—2002-211910 58
Figure 22. Temperature field comparison between HOTFGM-2D and analytical solutions for the

layered plate with exponentially increasing (top) and decreasing (bottom) piece-wise uniformly

graded thermal conductivity.

NASA/CR—2002-211910 59
dimensions, shown in the upper right schematic of Figure 21, and assigning the same values of

thermal conductivity in each vertical column of subcells as in the case of the analytical solution.

The above boundary conditions and manner of grading will produce a temperature field which is

symmetric about the horizontal cross section halfway across the plateís height.

The comparison of the full field temperature distributions obtained from the analytical and

higher-order theory solutions is given in Figure 22 for the exponentially increasing (top) and decreasing (bottom) thermal conductivity variations. In both cases, excellent agreement is evident

throughout almost the entire plate with the exception of isolated regions in the vicinity of the

upper and lower right corners. These regions exhibit large temperature gradients caused by the

temperature discontinuities at the corners themselves. Local mesh refinement is therefore needed

to produce better correlation in these regions in the case of HOTFGM.

As is also observed in Figure 22, when the thermal conductivity increases in the x2 direction

according to the function k = 0:5e5x2 , the temperature field is shifted towards the left vertical

boundary with respect to the temperature field in the homogeneous plate subjected to the same

boundary conditions shown in Figure 13. Alternatively, when the thermal conductivity decreases

in the x2 direction according to the function k = 60e−5x2 ; the temperature field is shifted towards

the right vertical boundary. This can be explained by the Fourierís heat conduction law,

q

(βγ)

i = °k(βγ)ij @jT(βγ)

where the temperature field is seen to be modulated by the spatial variation in thermal conductivity.

5.2 Efficient Reformulation: Variable Thermal Conductivity

In this section, the efficient reformulation of the higher-order theory described in Section 3 is employed to generate the temperature field in a continuously graded plate whose thermal conductivity

variations are described by the exponential functions of the preceding section. To mimic the exponentially varying thermal conductivity in the x2 direction, three volumetric discretizations were

employed. In all cases the plate was discretized into 25 subcells along the x3 direction, while along

the x2 direction the plate was discretized into 10 and 12 uniformly spaced, and 12 nonuniformly

spaced subcells, resulting in two uniformly spaced 10 £ 25 and 12 £ 25 microstructures and one

nonuniformly spaced 12 £ 25 microstructure. The thermal conductivity k(βγ)

2 in each vertical column of (Ø°) subcells was varied linearly with the local subcell coordinate x1

(β)

2 such that continuity

was preserved from one column of subcells to another along the global x2 coordinate. The values

of the thermal conductivity at the subcell boundaries were calculated using the exponential functions employed in the preceding section, and subsequently employed in the construction of linear

approximation of the thermal conductivity k(βγ)

2 in each subcell. Figure 23 illustrates the piecewise linearly varying k(βγ)

2 in the manner that approximates the exponentially increasing thermal

conductivity for the three volume discretizations. The thermal conductivity in the x3 direction,

k(βγ)

3 , was kept constant in each column of subcells, but varied from one column of subcells to the

next by taking the average value of the thermal conductivities at the subcell boundaries calculated

from the exponential functions. The general expressions for variable thermal conductivities in the

x2 and x3 directions, k(βγ)

2 and k(βγ)3 , respectively, thus reduce to the following expressions,

k(βγ)

2 = k(βγ)20 + k(βγ)22 x1

(β)

2

k(βγ)

3 = k(βγ)30

NASA/CR—2002-211910 60
Figure 23. Thermal conductivity approximation by piece-wise linear variation.

The full field temperature distributions obtained for the 10£25 subcell mesh with the exponentially increasing (top) and decreasing (bottom) thermal conductivities, generated using the efficient

reformulation with locally linearly varying conductivities, are presented in the right column of Figure 24. The corresponding predictions generated using the user-friendly formulation with 25 £ 25

subcell mesh of the preceding section are included in the left column of Figure 24 for comparison.

The agreement is very good everywhere except in the upper and lower right corner regions where

the temperature jumps occur. Increasing the mesh size to uniformly spaced 12 £ 25 subcells improves the agreement in these corner regions, as seen in the left column of Figure 25. Retaining

the same mesh size but preferentially increasing the number of subcells in the corner regions where

large thermal gradients occur produces essentially the same temperature distributions in the entire

linearly graded plate as those in the discretely graded plate, as seen in the right column of Figure

25. This indicates that a smaller number of subcells with local linearly varying thermal conductivity is required to capture the essential character of the temperature field in a discretely graded

plate with piecewise uniform thermal conductivity.

NASA/CR—2002-211910 61
Figure 24. Temperature field comparison between user-friendly formulation and efficient

reformulation of HOTFGM-2D.

NASA/CR—2002-211910 62
Figure 25. Temperature field comparison between evenly and nonevenly distributed thermal

conductivity predicted by the efficient reformulation of HOTFGM-2D with linearly varying

thermal conductivity.

NASA/CR—2002-211910 63
5.3 A Thermal Barrier Coating Application

In this section, the temperature field in a graded thermal barrier coating subjected to a throughthickness temperature gradient is investigated using the efficient reformulation of the higher-order

theory. The zirconia coating is graded with metallic inclusions such that the volume fraction of the

inclusion phase increases with increasing distance from the hot surface of the coating exposed to the

temperature of 200◦C, Figure 26. The cold surface where the metallic inclusion content is greatest

is exposed to 100◦C. The coatingís graded microstructure, which is periodic along the vertical

direction, is simulated by the zero heat flux boundary conditions on the top and bottom surfaces

of the coating. The thermal conductivities of the zirconia coating and the metallic inclusions were

taken to be 0:5 and 100 W=m °◦ C, respectively.

Figure 26. Actual microstructure and volume fraction distribution of the inclusion phase.

The thermal analysis is conducted using the actual graded microstructure seen in Figure 26

which has been discretized into 10 £ 100 subcell mesh, and two additional discretely graded and

linearly graded microstructures with local homogenization. The locally homogenized microstructures were obtained from the actual microstructure by first determining the local metallic inclusion

content in each of the nine regions into which the actual microstructure was divided in the manner

shown in Figure 26. Subsequently, each region with its local inclusion content was modeled using

HOTFGM as a strip containing 16 uniformly spaced metallic inclusions subjected to a known macroscopic temperature gradient along its length and zero heat flux on the top and bottom surfaces,

Figure 27.

NASA/CR—2002-211910 64
Figure 27. Homogenized thermal conductivity vs fiber volume fraction.

The effective thermal conductivity for such strip was calculated from the ratio of the average heat

flux obtained from the solution of the specified boundary-value problem and the negative of the

imposed average or macroscopic temperature gradient. The choice of 16 uniformly spaced metallic

inclusions was dictated by the asymptotic behavior of the average thermal conductivity as a function of the inclusion number. The results of the calculations are included in Figure 27 in the form

of a graph of the macroscopic thermal conductivity as a function of the metallic inclusion volume

fraction which corresponds to each of the regions shown in Figure 26. The resulting discretely and

linearly graded in a piecewise manner thermal conductivities along the actual coatingís length are

given in Figure 28. The piecewise linearly graded microstructure was obtained from the piecewise

discretely graded microstructure by using the center points of adjacent subcells to define new subcells with linearly varying conductivities. Due to the homogenization of the microstructure and the

applied boundary conditions which produce a one-dimensional temperature field, the homogenized

microstructures were represented by 1 £ 9 subcell meshes.

The full field temperature distributions in the actual and homogenized graded coatings are

given in Figure 29. The temperature distribution in the actual coating is clearly two-dimensional,

with the disturbances created by the inclusion presence more pronounced near the upper and lower

boundaries than in the interior. In contrast, the temperature fields in the homogenized coatings

are clearly one-dimensional, with a slight shift present between the discretely and linearly graded

homogenized microstructures. A more quantitative picture of the differences in the temperature distributions among the three microstructures is obtained by comparing the cross section plots. These

are presented in the top portion of Figure 30 with three cross section temperature distributions

NASA/CR—2002-211910 65
Figure 28. Homogenized thermal conductivity as a function of position.

in the actual coating and one in each of the two homogenized coatings. The staircase temperature pattern observed in the actual coating is due to the pronounced differences in the thermal

conductivities of the metallic inclusion and ceramic matrix phases. In contrast, the temperature

distributions in the homogenized coatings are smoothly and continuously increasing functions of

the coordinate along the coatingís length, with the differences between the two distributions visible

near the hot surface. The bottom portion of Figure 30 compares the temperature distribution in

the actual coating averaged across the coatingís width with the distributions in the homogenized

coatings. In both cases, the temperature distributions in the actual coating, averaged or unaveraged

across the coatings width, are above the temperature distribution in the discretely graded coating.

This is not true in the case of the linearly graded coating wherein the temperature distribution is

actually higher in some regions near the hot surface than the temperature distribution in the actual

coating.

Despite the fact that the temperature distributions in the two homogenized coatings are almost

identical with the exception of the region near the hot surface, the average heat flux calculated

in the homogenized coating with piecewise linear variation of the thermal conductivity is closer to

the average heat flux in the actual coating than the averag flux in the discretely graded coating.

These values can be used to calculate the macroscopic thermal conductivity of the coating across

its thickness as shown below,

q1

actual

2 = °9:787 =) 1

kactual

2 = q1actual

2

¢T =L = 0:9787

q1

linear

2 = °9:768 =) 1

klinear

2 = q1linear

2

¢T =L = 0:9768

q1

discrete

2 = °9:841 =) 1

kdiscrete

2 = q1discrete

2

¢T =L = 0:9841

NASA/CR—2002-211910 66
Figure 29. Temperature distributions in TBC with actual and homogenized microstructures.

NASA/CR—2002-211910 67
Figure 30. Cross-section temperature distributions comparison.

NASA/CR—2002-211910 68
where ¢T = 100◦C is the temperature difference between the top and bottom surfaces of the

coating and L = 100 is the coatingís thickness. Consequently, the effective thermal conductivity

of the entire coating based on the linearly graded microstructure is closer to the effective thermal

conductivity of the actual coating than the discretely homogenized coating. Thus it appears that

the linear thermal conductivity variation is more appropriate for estimating the effective thermal

conductivities of actual graded microstructures than the discrete thermal conductivity variation

when the coatingís microstructure is relatively coarse.

6 Conclusions and Future Work

The goal of this research was to develop a generic cell-free, subcell-based formulation of the higherorder theory for functionally graded materials developed by Aboudi et al. (1999), which would

then provide a basis for a more efficient reformulation using the local/global conductivity matrix

approach. In order to extend the higher-order theoryís range of applicability, variable thermal

conductivity capability at the subcell level was also incorporated into the efficient reformulation.

This investigation was undertaken because the original formulation of HOTFGM is computationally

intensive. This, in turn, limits the size of problems that can be analyzed due to the large number

of equations required to mimic realistic microstructures of functionally graded materials.

In order to validate both the user-friendly formulation and the reformulation of HOTFGM with

variable thermal conductivity, analytical solutions to the steady-state heat equation were developed for homogeneous and layered plates. The homogeneous plate solution was used to study the

convergence behavior of the user-friendly formulation as a function of mesh discretization. The

layered-plate solution was employed to investigate temperature distributions in discretely graded

plates. Comparison between the results obtained from the user-friendly formulation and analytical

solutions validates the simplified set of equations of the higher-order theory based on the more

straightforward volume discretization which does not differentiate between generic cells and subcells. In particular, identical results were obtained for a homogeneous plate subjected to thermal

boundary conditions that produced internal temperature fields with different symmetries. The solutionís convergence was demonstrated to be quite rapid. A piecewise uniformly graded plate with

exponentially varying thermal conductivity was also analyzed and the results obtained from the

higher-order theory coincided with the analytical solution for a layered plate.

The user-friendly formulation is easy to understand and implement. More importantly, it provides a basis for the development of an efficient reformulation based on the local/global conductivity

matrix approach, which was also successfully accomplished, taking the locally variable thermal conductivity at the subcell level into account. The efficient reformulation has significantly decreased

the number of equations. The reduction of the number of equations by as much as 60% for some

problems has enhanced the theoryís capability, enabling solution of more practical and complicated

problems. Comparison of temperature fields in piecewise uniformly graded and linearly graded

cases, both of which approximate the same spatially varying thermal conductivity at the global

scale, revealed that the same temperature distribution can be obtained with a smaller volume discretization when the thermal conductivity is allowed to vary at the subcell level. Therefore, the

efficient reformulation of the higher-order theory with variable thermal conductivity made possible

the solution of a wide range of problems requiring large number of subcells. Thus, after transition from the user-friendly formulation to reformulation, the temperature distribution prediction

for arbitrary microstructures under many circumstances can be accurately accomplished in a more

efficient manner.

NASA/CR—2002-211910 69
However, this research was limited to steady-state heat conduction of graded materials in two

dimensions. In order to further extend the theoryís range of applicability to practical problems,

further research is required. The next step is to extend the presented efficient reformulation to

enable the analysis of mechanical problems. As the first step, the two-dimensional elastic reformulation of the higher-order theory with linearly varying material properties within each subcell

will be developed and then extended to three dimensions. This will complement the efficient thermomechanical reformulation of the higher-order theory with constant subcell material properties

recently completed by Bansal (2002). Further, in order to make full use of material properties and

optimize structural design, it is necessary to incorporate viscoelastic, viscoplastic or plastic effects

of material response into the theoretical framework. Therefore, viscoelastic, viscoplastic and plastic

capabilities should be included in the reformulated higher-order theory. Fracture analysis is also an

important area of FGMs. Therefore, the capability to analyze fracture mechanics problems should

also be included in the reformulated higher-order theory.

References

[1] Aboudi, J. (1991). Mechanics of Composite Materials: A Unified Micromechanical Approach,

Elsevier, Amsterdam.

[2] Aboudi, J., Pindera, M-J., Arnold, S. M. (1993). Thermoelastic Response of Large-Diameter

Fiber Metal Matrix Composites to Thermal Gradients. NASA TM 106344, NASA-Lewis Research Center, Cleveland, OH.

[3] Aboudi, J., Pindera, M-J., Arnold, S. M. (1994). Elastic Response of Metal Matrix Composites with Tailored Microstructures to Thermal Gradients. International Journal of Solids

Structures, Vol. 31, No. 10, pp 1393-1428.

[4] Aboudi, J., Pindera, M-J., and Arnold, S. M. (1996). Thermoelastic Theory for the Response

of Materials Functionally Graded in Two Directions, International Journal of Solids and Structures, Vol. 33, No. 7, pp. 931-966.

[5] Aboudi, J., Pindera, M-J., and Arnold, S. M. (1999). Higher-Order Theory for Functionally

Graded Materials, Composites Part B: Engineering, Vol. 30, No. 8, pp. 777-832.

[6] Amada, S. (1995). Hierarchical Functionally Graded Structures of Bamboo, Barley, and Corn,

Material Research Society Bulletin, January 1995, Volume XX, No.1. pp. 35-36.

[7] Bansal, Y. (2002). Efficient Reformulation of the Higher-Order Theory for Functionally Graded

Materials, M. S. Thesis, School of Engineering and Applied Science, University of Virginia,

Charlottesville, Virginia.

[8] Davison, T., Wadley, H. N. G, and Pindera, M-J. (1994). Elastic Bahavior of a Layered Cylinder

Subjected to Diametral Loading, Composites Engineering, Vol. 4, No. 10, pp. 995-1009.

[9] Delfosse, D., Cherradi, N., and IIschner, B. (1997). Numerical and Experimental Determination

of Residual Stresses in Graded Materials, Composites Part B: Engineering, Vol. 28, No. 1-2,

pp. 127-141.

NASA/CR—2002-211910 70
[10] Drake, J. T., Williamson, R. L., and Rabin, B. H. (1993). Finite Element Analysis of Thermal

Residual Stresses at Graded Ceramic-Metal Interfaces, Journal of Applied Physics, Vol. 74,

No. 2, pp. 1321-1326.

[11] Delale, F. and Erdogan, F. (1983). The Crack Problem for a Non-Homogeneous Plane, Journal

of Applied Mechanics, Vol. 50, pp. 609-614.

[12] Dollmeier, K., Thiele, W., and IIschner, B. (1995), Pilot Study Related to Continuous Fabrication of FGM Strip by P/M, in FGM 94, Proceedings of the Third International Symposium

on Structural and Functionally Graded Materials, IIschner, B. and Cherradi, N. (Eds.), pp.

27-32, Presses Polytecniques et Universitares Romandes.

[13] Erdogan, F. (1985). The Crack Problem for Bonded Non-Homogeneous Materials under AntiPlane Shear Loading, Journal of Applied Mechanics, Vol. 52, pp. 823-828.

[14] Erdogan, F., Kaya, A. C., and Joseph, P. F. (1991). The Crack Problem in Bonded NonHomogeneous Materials, Journal of Applied Mechanics, Vol. 58, No. 2, pp. 410-418.

[15] Erdogan, F. (1995). Fracture Mechanics of Functionally Graded Materials, Composites Engineering, Vol. 5, No. 7, pp. 753-770.

[16] Erdogan, F. and Wu, B. H. (1996). Crack Problems in FGM Layers under Thermal Stresses,

Journal of Thermal Stresses, Vol. 19, No. 3, pp. 237-265.

[17] Fukumoto, M., Yamasaki, T. and Okane, I. (1994). Influence of Coating Configuration on

Cyclic Thermal Shock Fracture Behaviour of Plasma Sprayed Coatings, in FGM 94, Proceedings of the Third International Symposium on Structural and Functionally Graded Materials,

IIschner, B. and Cherradi, N. (Eds.), pp. 425-429, Presses Polytecniques et Universitares Romandes.

[18] Horgan, C. O. and Chan, A. M. (1999). The Pressurized Hollow Cylinder or Disk Problem for

Functionally Graded Isotropic Linearly Elastic Materials, Journal of Elasticity, Vol. 55, No. 1,

pp. 43-59.

[19] Ichikawa K. (Ed.), (2001). Functionally Graded Materials in the 21st Century (A Workshop

on Trends and Forecasts), Kluwer Academic Publishers.

[20] Ilschner, B. and Cherradi, N., (1994). Proceedings of the 3rd International Symposium on

Structural and Functional Gradient Materials, Presses polytechniques et Universitaires Romandes.

[21] Ilschner, B. (1997). Lessons Learnt in 7 years of FGM Research at Lausanne in FGM 96,

Proceedings of the Fourth International Symposium on Functionally Graded Materials, Shiota,

I. and Miyamoto, Y. (Eds.), pp. 15-20, Amsterdam, Elsevier.

[22] Jian, C. Y., Hashida, T., Takahashi, H., and Saito, M. (1995). Thermal Shock and Fatigue Resistance Evaluation of Functionally Graded Coating for Gas Turbine Blades by Laser Heating

Method, Composites Engineering, Vol. 5, No. 7, pp. 879-889.

NASA/CR—2002-211910 71
[23] Jin, Z. H. and Paulino, G. H. (2001). Transient Thermal Stress Analysis of an Edge Crack in

a Functionally Graded Material, International Journal of Fracture, Vol. 107, No. 1, pp. 73-98.

[24] Kaysser, W. A. and IIschner, B. (1995). FGM Research Activities in Europe, Material Research

Society Bulletin, January 1995, Volume XX, No.1. pp. 22-26.

[25] Kawasaki, A. and Watanabe, R. (1997). Effect of Gradient Microstructure on Thermal Shock

Crack Extension in Metal/Ceramic Functionally Graded Materials, in FGM 96, Proceedings of the Fourth International Symposium on Functionally Graded Materials, Shiota, I and

Miyamoto, Y. (Eds.), pp. 143-154, Amsterdam, Elsevier.

[26] Kokini, K. and Choules, B. B. (1995). Surface Thermal Fracture of Functionally Graded Ceramic Coatings-Effect of Architecture and Materials, Composites Engineering, Vol. 5, No. 7,

pp. 865-877.

[27] Kokini, K., Takeuchi, Y. R., and Choules, B. B. (1996). Surface Thermal Cracking of Thermal

Barrier Coating Owing to Stress Relaxation: Zirconia vs Mullite, Surface Coating Technology,

Vol. 82, pp. 72-82.

[28] Lee, Y-D. and Erdogan, F. (1998). Interface Cracking of FGM Coatings Under Steady-State

Heat Flow, Engineering Fracture Mechanics, Vol. 59, No. 3, pp. 361-380.

[29] Needleman, A. and Suresh S. (1996). Mechanics and Physics of Layered and Graded Materials,

Journal of the Mechanics and Physics of Solids, Vol. 44, No. 5, pp. 647-821.

[30] Nomura, S. and Sheahen D.M. (1997). Micromehanical Approach to the Thermomechanical

Analysis of FGMs, in FM 96, Proceedings of the Fourth International Symposium on Functionally Graded Materials, Shiota, I and Miyamoto, Y. (Eds.), pp. 137-141, Amsterdam, Elsevier.

[31] Noda, N. (1999). Thermal Stresses in Functionally Graded Materials, Journal of Thermal

Stresses, Vol. 24, No. 4-5, pp. 477-512.

[32] Paley, M. and Aboudi. J. (1992). Micromechanical Analysis of Composites by the Generalized

Cells Model. Mechanics of Materials, Vol. 14, pp. 39-127.

[33] Pindera, M-J. (1991). Local/Global Stiffness Matrix Formulation for Composite Materials and

Structures, Composites Engineering, Vol. 1, No. 2, pp. 69-83.

[34] Pindera, M-J. and Lane, M. S. (1993). Frictionless Contact of Layered Half Planes, Part I:

Analysis, Journal of Applied Mechanics, Vol. 30, No. 3, pp. 633-639.

[35] Pindera, M-J., Freed, A. D., and Arnold, S. M. (1993). Effects of Fiber and Interfacial Layer

Morphologies on the Thermoplastic Response of Metal Matrix Composites, International Journal of Solids and Structures, Vol. 30, No. 9, pp. 1213-1238.

[36] Pindera, M-J., Arnold, S. M., Aboudi, J., and Hui, D. (Eds.), (1994a). Use of Composites in

Functionally Graded Materials, Composites Engineering, Vol. 4, No. 1. pp. 1-145.

[37] Pindera, M-J., Williams, T. O. and Arnold, S. M. (1994b). Thermoplastic Response of Metalmatrix Composites with Homogenized and Functionally Graded Interfaces, Composites Engineering, Vol. 4, No.1, pp 129-145.

NASA/CR—2002-211910 72
[38] Pindera, M-J., Aboudi, J., Arnold, S. M., and Jones, W. F. (Eds.), (1995a). Use of Composites

in Muti-Phased and Functionally Graded Materials, Composites Engineering, Vol. 5, No. 7.

pp. 743-974.

[39] Pindera, M-J., Aboudi, J., and Arnold, S. M. (1995b). Limitations of the Uncoupled, RVEbased Micromechanical Approach in the Analysis of Functionally Graded Composites, Mechanics of Materials, Vol. 20, No. 1, pp. 77-94.

[40] Pindera, M-J., Aboudi, J., and Arnold, S. M. (1995c). Recent Advances in the Mechanics of

Functionally Graded Composites, in Aerospace Thermal Structures and Materials for a New

Era, Progress in Astronautics and Aeronautics, 168 . Washington, DC, Thornton, E. A. (Ed.),

American Institute of Aeronautics and Astronautics Inc. pp. 181-203.

[41] Pindera, M-J., Aboudi, J., Glaeser, A. M., and Arnold, S. M. (Eds.), (1997). Use of Composites

in Muti-Phased and Functionally Graded Materials, Composites Engineering, Vol. 28, No. 1/2,

pp. 1-175.

[42] Pindera, M.-J. and Dunn, P. (1997). Evaluation of the Higher-Order Theory for Functionally

Graded Materials via the Finite Element Method, Composites: Part B: Engineering, Vol 28,

No. 1/2, pp. 109-119.

[43] Pindera, M-J., Aboudi, J., and Arnold, S. M. (1998). Higher-Order Theory for the Analysis

of Functionally Graded Materials, in Advanced Multilayered and Fiber-Reinforced Composites,

Haddad, Y. M. (Ed.), pp. 111-132, Kluwer Academic Publishers, The Netherlands.

[44] Pindera, M-J. and Aboudi, J. (2000). HOTFGM-2D: A Coupled Higher-Order Theory for

Cylindrical Structural Components with Bi-Directionally Graded Microstructures, NASA Contractor Report 2000-210350, NASA-Glenn Research Center, Cleveland, OH.

[45] Salzar, R. S. and Barton, F. W. (1994). Residual Stress Optimization in Metal-Matrix Composites using Discretely Graded Interfaces, Composites Engineering, Vol. 4, No. 1, pp. 115ó128.

[46] Salzar, R. S., Pindera, M-J., and Barton, F. W. (1996). Elastic/Plastic Analysis of Layered

Metal Matrix Composite Cylinders. Part I: Theory, Journal of Pressure Vessel Technology,

Vol. 118, No. 1, pp. 13-20.

[47] Shiota, I. and Miyamoto, Y., (1997). Proceedings of the 4th International Symposium on

Functionally Graded Materials, Functionally Graded Materials 1996, Elsevier.

[48] Suresh, S., Giannakopoulos, A. E., and Olsson, M. (1994). Elastoplastic Analysis of Thermal

Cycling-Layered Materials with Sharp Interfaces, Journal of the Mechanics and Physics of

Solids, Vol. 42, No. 6, pp. 979-1018.

[49] Suresh, S. and Mortensen A., (1998). Fundamentals of Functionally Graded Materials, 2 nd

edition, Cambridge University Press.

[50] Sutradhar, A., Paulino, G. H., and Gray, L. J. (2002). Transient Heat Conduction in Homogeneous and Non-Homogeneous Materials by the Laplace Transform Galerkin Boundary Element

Method, Engineering Analysis with Boundary Elements, Vol. 26, No. 2, pp. 119-132.

[51] Urquhart, E. E. and Pindera, M-J. (1994). Incipient Separation Between a Frictionless Flat

Punch and an Anisotropic Multilayered Half Plane, International Journal of Solids and Structures, Vol. 31, No. 18, pp 2445-2461.

NASA/CR—2002-211910 73
[52] Williamson, R. L., Rabin, B. H., and Drake, J. T. (1993). Finite Element Analysis of Thermal Residual Stresses at Graded Ceramic-Metal Interfaces. Part 1. Model Description and

Geometrical Effects, Journal of Applied Physics, Vol. 74, No. 2, pp. 1310-1320.

[53] Williams, T. O. and Pindera, M-J. (1997). An Analytical Model for the Inelastic Axial Shear

Response of Unidirectional Metal Matrix Composites, International Journal of Plasticity, Vol.

13, No. 3, pp. 261-289.

[54] Yamanouchi, M., Koizumi, M., Hirai, T. and Shiota, I. (Eds.), (1990). FGM 90, Proceedings

of the First International Symposium on Functionally Graded Materials, October 8-9, 1990,

Sendai, Japan.

NASA/CR—2002-211910 74
This publication is available from the NASA Center for AeroSpace Information, 301–621–0390.

REPORT DOCUMENTATION PAGE

2. REPORT DATE

19. SECURITY CLASSIFICATION

 OF ABSTRACT

18. SECURITY CLASSIFICATION

 OF THIS PAGE

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this

collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson

Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102

Form Approved

OMB No. 0704-0188

12b. DISTRIBUTION CODE

8. PERFORMING ORGANIZATION

 REPORT NUMBER

5. FUNDING NUMBERS

3. REPORT TYPE AND DATES COVERED

4. TITLE AND SUBTITLE

6. AUTHOR(S)

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT

13. ABSTRACT (Maximum 200 words)

14. SUBJECT TERMS

17. SECURITY CLASSIFICATION

 OF REPORT

16. PRICE CODE

15. NUMBER OF PAGES

20. LIMITATION OF ABSTRACT

Unclassified Unclassified

Final Contractor Report

Unclassified

1. AGENCY USE ONLY (Leave blank)

10. SPONSORING/MONITORING

 AGENCY REPORT NUMBER

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546–0001

Available electronically at http://gltrs.grc.nasa.gov

November 2002

NASA CR—2002-211910

E–13600

WU–708–73–35–00

NAG3–2524

80

Efficient Reformulation of HOTFGM: Heat Conduction With Variable

Thermal Conductivity

Yi Zhong and Marek-Jerzy Pindera

Micromechanics; Elastic; Thermal; Deformation; Composites; Analysis

Unclassified -Unlimited

Subject Categories: 24 and 39 Distribution: Nonstandard

University of Virginia

0 Mccormick Road

Charlottesville, Virginia 22904–2000

Project Manager, Steven M. Arnold, Structures and Acoustics Division, NASA Glenn Research Center, organization code

5920, 216–433–3334.

Functionally graded materials (FGMs) have become one of the major research topics in the mechanics of materials community during

the past fifteen years. FGMs are heterogeneous materials, characterized by spatially variable microstructures, and thus spatially variable

macroscopic properties, introduced to enhance material or structural performance. The spatially variable material properties make

FGMs challenging to analyze. Two distinct and fundamentally different computational strategies, called uncoupled macromechanical

and coupled micro-macromechanical approaches by some investigators, have been pursued. The uncoupled macromechanical approaches ignore the effect of microstructural gradation by employing specific spatial variations of material properties, which are either

assumed or obtained by local homogenization, thereby resulting in erroneous results under certain circumstances. In contrast, the

coupled approaches explicitly account for the micro-macrostructural interaction, albeit at a significantly higher computational cost. The

higher-order theory for functionally graded materials (HOTFGM) is representative of the coupled approach. However, despite its

demonstrated utility in applications where micro-macrostructural coupling effects are important, the theory’s full potential is yet to be

realized because the original formulation of HOTFGM is computationally intensive. This, in turn, limits the size of problems that can

be solved due to the large number of equations required to mimic realistic material microstructures. Therefore, a basis for an efficient

reformulation of HOTFGM, referred to as user-friendly formulation, is developed herein, and subsequently employed in the construction of the efficient reformulation using the local/global conductivity matrix approach.