515.44.843.17: 44:.
u... rs.‘ 1 .

 

. . 5v.
9 .E 5.1.! r. 5:...
.I. .cfiux... a
2.; L . fir. ,e ...
unhawra. x 1.. . 5
. . w. La. 2.0:”,
V g .7 .5

V, .57 7310.5.

A _:ve..vl.,
7 v
t

1:,
«(v K.“

 

 

q. , . .l. :f, ‘ .01. . .3 -1»: ! W151»: ?

 

lllllllllllllllllllllllllllillIllllllllllllll

kw.
. ‘f 3 1293 02080
I’ (/1

This is to certify that the
thesis entitled

A NEW THERMAL THEORY AND ITS USE IN THERMO-
MECHANICAL ANALYSIS OF LAMINATED COMPOSITES
AND SANDWICH PANELS

presented by
Antonio Pantano

has been accepted towards fulfillment
of the requirements for

f(. S. degree in HFCHA’WCS

 

MKW

Major professor

Date (/7/30’/?7

0-7639 MS U is an Affirmatiw Action/Equal Opportunity Institution

 

PLACE IN RETURN BOX to remove this checkout from your record.

 

 

LIBRARY

Micm'Qan State
UH‘VGFSEty

f

 
  
    

To AVOID FINE-S return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE

DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

moo mummwspu

A NEW THERMAL THEORY AND ITS USE IN THERMO-
MECHANICAL ANALYSIS OF LAMINATED COMPOSITES AND
SANDWICH PANELS

By

Antonio Pantano

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Materials Science and Mechanics

1999

ABSTRACT

A NEW THERMAL THEORY AND ITS USE IN THERMO-
MECHANICAL ANALYSIS OF LAMINATED COMPOSITES AND
SANDWICH PANELS

By

Antonio Pantano

A new thermal lamination theory and its associated finite element model are
presented for analysis of heat transfer in laminated composite structures. The form of the
present theory closely resembles that of recent zig-zag sublarninate structural laminate
theories. The through-thickness distribution of temperature is assumed to vary linearly
within each ply, and continuity of transverse flux at ply interfaces is enforced
analytically. This theory lends itself well to development of efficient finite element
models. Non-linear capabilities are implemented in the model to solve heat conduction
problems in the presence of material non-linearity.

Combining the presented thermal theory with a zig-zag structural theory a 3D
finite element model for thermal stress analysis of composite laminates has been realized.
The novel features of the formulation allow accurate results in predicting the distribution
of temperatures, displacements and stresses in laminated plates wherein the plies have
dissimilar thermal and/or structural properties. Both linear and nonlinear numerical

results are presented.

Cepyright by
Antonio Pantano

1999

To my Mother and Father

iv

ACKNOWLEDGMENTS

I will always be grateful to my beloved parents for years of unconditional support.
I hope they will be proud of me, now and in the future.

I would like to reserve a special gratitude to my sweet Giulia. Without her love
this result would have been harder to achieve.

I definitely recognize the immense fortune I had to meet a person like my research
advisor Prof. Ronald Averill. His uncommon humanity equals only his exceptional
academic knowledge.

I also greatly appreciate Prof. Robert WM. Soutas-Little and Prof. Tomas J. Pence

for making efforts to serve my thesis committee.

Antonio Pantano
November 1999

East Lansing, Michigan

TABLE OF CONTENTS

LIST OF TABLES ......................................................................................... ix
LIST OF FIGURES ........................................................................................ x
CHAPTER 1 INTRODUCTION ................................................................. l
1.1 Preliminary Information. ........................................................ 1

1.2 Literature Review ................................................................... 2

1.2.1 Equivalent Single Layer Approach .......................................... 2

1.2.2 FSDT with Postprocessing Techniques ................................... 4

1.2.3 Discrete-Layer Theories .......................................................... 4

1.2.4 Linear Flux Theory .................................................................. 4

1.2.5 Zig-Zag Laminate Theory ........................................................ 5

1.3 Objective of the Present Study ............................................... 6

1.4 Organization of the Thesis ..................................................... 7
CHAPTER 2 THERMAL LAMINATION THEORY ................................ 9
2.1 Formulation of the Theory ..................................................... 9

2.1.1 Demonstration of Zig-Zag Effects ......................................... 13

2.2 Finite Element Model. .......................................................... 16

2.2.1 Possible Forms ....................................................................... 16

2.2.2 Four-Noded Layered Rectangular Thermal Element ............ 19

2.3 The Computer Program “Therm2D” ................................... 21

2.4 The Computer Program “Sinload” ....................................... 22

2.5 Numerical Results ................................................................ 23

2.5.1 Symmetric 8 Layer Laminate ................................................ 24

vi

CHAPTER 3
3. 1

3.2
3.3

CHAPTER 4
4.1

4.2

CHAPTER 5
5. l

5.2

CHAPTER 6

APPENDIX A

APPENDIX B

2.5.2 Unsymmetric 6 Layer Laminate ............................................ 29

2.5.3 Shuttle Thermal Protection System ....................................... 33
NONLINEAR THERMAL MODEL ................................... 41
Basic Relations ..................................................................... 41
The Computer Program “ThermNL” ................................... 42
Numerical Results ................................................................ 43
3.3.1 Mixed Boundary Conditions Analysis ................................... 46
3.3.1.1 Growing Thermal Conductivity coefficients .................. 46
3.3.1.2 Decreasing Thermal Conductivity coefficients ............... 49
3.3.2 Essential Boundary Conditions Analysis ............................... 53
3.3.2.1 Growing Thermal Conductivity coefficients. .................. 53
3.3.2.2 Decreasing Thermal Conductivity coefficients ............... 54
FIRST ORDER ZIG-ZAG PLATE THEORY .................... 58
Formulation of the Theory ................................................... 58
Finite Element Model ........................................................... 64
THERMAL STRESS MODEL ............................................ 66
Equivalent Nodal Force Vector Method .............................. 66
Numerical Results ................................................................ 68
5.2.1 Sandwich Panel ...................................................................... 69
5.2.1.1 Demonstration of Thermal Zig-Zag Effects .................... 70
5.2.1.2 Sandwich with Aspect Ratio a/h=20 ............................... 72
5.2.1.3 Sandwich with Aspect Ratio a/h=4 ................................. 81
5.2.2 Four-Layered Antisyrnmetric Cross-Ply Laminate ............... 90
CONCLUSIONS .................................................................. 98
FORTRAN CODE, THERMZD ............................ 101
FORTRAN CODE, SINLOAD ............................. 129

vii

APPENDIX C FORTRAN CODE, THERMNL ........................... 142

REFERENCES ........................................................................................... 1 82

viii

Table 1.
Table 2.

Table 3.

Table 4.

Table 5.

Table 6.

LIST OF TABLES

Variation of temperature versus ply thermal conductivity ratio. ................... 16
Material properties used in the numerical analyses. ...................................... 25

Temperature distribution along the thickness at x1 = 2.5 for the 8 layer
symmetric plate with aspect ratio of 10. (Temperatures are normalized by the
exact solution.) ............................................................................................... 26

Temperature distribution along the thickness at x1= 5 for the 8 layer
symmetric plate with aspect ratio of 10. (Temperatures are normalized by the
exact solution.) ............................................................................................... 26

Shuttle thermal protection system properties adopted from [12]. ................. 33

Comparison among the temperature distributions along the thickness at x1= 5
at the end of the nonlinear analysis ................................................................ 54

ix

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6
Figure 7

Figure 8

Figure 9

Figure 10

Figure 11

Figure 12

Figure 13

LIST OF FIGURES

Geometry and coordinate system used in the thermal sublarninate theory. .. 10

Through-thickness variation of the temperature of a three-ply laminate
subjected to T=1 at the top surface and T = O at the bottom. ........................ 15

Quadrilateral and triangular finite elements based on the zig-zag thermal
lamination theory. Each node contains two computational DOFs. ............... 17

(a) Eight-noded brick and six-noded wedge three-dimensional finite elements
based on the zigzag thermal lamination theory. Each node contains one
computational DOF - the temperature at that node. (b) An example of
discretizing the thickness of a laminate using multiple elements
(sublarninates). ............................................................................................... 17

Element topology and nodal degrees of freedom for the three-dimensional

form of the thermal model. ............................................................................ 18
Four-noded two-dimensional sublarninate finite element model .................. 18
BC’s for the symmetric 8-layer laminate ....................................................... 24
Through-thickness distribution of temperature at x1 = 2.5 for an 8-layer
symmetric laminate with aspect ratio of 10. .................................................. 27
Through-thickness distribution of temperature at x1 = 5 for an 8-layer
symmetric laminate with aspect ratio of 10. .................................................. 28
Through-thickness distribution of temperature at x1 = 2 for an 8-layer
symmetric laminate with aspect ratio of 4. .................................................... 30
Through-thickness distribution of temperature at x1 = 5 for a 6-layer
unsyrnmetric laminate with aspect ratio of 10. .............................................. 31
Through-thickness distribution of temperature at x1 = 2 for a 6-1ayer
unsyrnmetric laminate with aspect ratio of 4. ................................................ 32
Shuttle thermal protection system. ................................................................ 34

Figure 14

Figure 15

Figure 16

Figure 17

Figure 18

Figure 19

Figure 20

Figure 21
Figure 22
Figure 23
Figure 24

Figure 25

Figure 26

Figure 27

Figure 28

Figure 29

Through-thickness distribution of temperature at x1 = 5 for the Shuttle
thermal protection system with aspect ratio of 10. ........................................ 35

Through-thickness distribution of temperature at x, = 5 for the Shuttle
thermal protection system with aspect ratio of 10. Zoom view showing the

last 4 layers. ................................................................................................... 36
Through-thickness distribution of temperature at x] = 2 for the Shuttle
thermal protection system with aspect ratio of 4. .......................................... 37

Through-thickness distribution of temperature at X] = 2 for the Shuttle
thermal protection system with aspect ratio of 4. Zoom view showing the last
4 layers. .......................................................................................................... 38

Through-thickness distribution of temperature at x1 = 5 for the Shuttle
thermal protection system with aspect ratio of 10. Test case with top and
sides of the plate insulated. ............................................................................ 39

Through-thiclmess distribution of temperature at x1 = 5 for the shuttle thermal
protection system with aspect ratio of 10. Test case with top and sides of the
plate insulated. Zoom view showing the last 4 layers. .................................. 40

Points of the sublaminate element where the temperature needs to be

evaluated. ....................................................................................................... 42
Thermal conductivities variation with T (°C) for M and 3.2 positive ............. 44
Thermal conductivities variation with T (°C) for M and A2 negative. ........... 45
Model for mixed B.C.’s analysis ................................................................... 46
Final through-thickness distribution of temperature at x1 = 2.5 for values of
the conductivity coefficients k increasing with T. ......................................... 47
Final through-thickness distribution of temperature at X] = 5 for values of the
conductivity coefficients k increasing with T ................................................ 48
Solution at every increments at x1 = 5 and x3 = l for values of the

conductivity coefficient k increasing with T. ................................................ 50
Solution at every increments at x1 = 5 and X3 = 0.5 for values of the
conductivity coefficient k increasing with T. ................................................ 51

Solution at every increments at X] = 5 and X3 = l for values of the
conductivity coefficient k decreasing with T ................................................. 52

Model for essential B.C.’s analysis. .............................................................. 53

xi

Figure 30

Figure 31

Figure 32

Figure 33
Figure 34

Figure 35

Figure 36

Figure 37

Figure 38

Figure 39

Figure 40

Figure 41

Figure 42

Figure 43

Figure 44

Figure 45

Solution at every increments at x1 = 5 and X3 = 0.5 for values of the
conductivity coefficient k increasing with T. ................................................ 55
Final through-thiclmess distribution of temperature at X] = 5 for values of the
conductivity coefficient k increasing with T. ................................................ 56
Solution at every increments at x. = 5 and x3 = 0.5 for values of the
conductivity coefficient k decreasing with T ................................................. 57
Element topology and nodal degrees of freedom of the structural model. 65
Through-thickness distribution of temperature at x. = X2 = a/2 in the

sandwich plate with aspect ratio 20. .............................................................. 71

Through-thickness variation of the displacement 111 at x1 = O and X2 = a/2 in

the sandwich plate with aspect ratio 20. ........................................................ 73
Through-thickness variation of the displacement u2 at x1 = a/2 and X2 = O in
the sandwich plate with aspect ratio 20. ........................................................ 74

Through-thickness variation of the stress 0'11 at x1 = a/2 and X2 = a/2 in the
sandwich plate with aspect ratio 20. .............................................................. 75

Through-thickness variation of the stress 022 at x1 = a/2 and X2 = a/2 in the
sandwich plate with aspect ratio 20. .............................................................. 76

Through-thickness variation of the stress 033 at x] = a/2 and X2 = a/2 in the
sandwich plate with aspect ratio 20. .............................................................. 77

Through-thickness variation of the stress 023 at X] = a/2 and X2 = O in the
sandwich plate with aspect ratio 20. .............................................................. 78

Through-thickness variation of the stress 013 at x1 = O and X2 = a/2 in the
sandwich plate with aspect ratio 20. .............................................................. 79

Through-thickness variation of the stress 0'12 at x] = O and X2 = O in the
sandwich plate with aspect ratio 20. .............................................................. 8O

Through-thickness variation of the displacement ul at x, = 0 and X2 = a/2 in
the sandwich plate with aspect ratio 4. .......................................................... 82

Through-thickness variation of the displacement u2 at x1 = a/2 and X2 = O in
the sandwich plate with aspect ratio 4. .......................................................... 83

Through-thickness variation of the stress on at x] = a/2 and X2 = a/2 in the
sandwich plate with aspect ratio 4. ................................................................ 84

xii

Figure 46

Figure 47

Figure 48

Figure 49

Figure 50

Figure 51

Figure 52

Figure 53

Figure 54

Figure 55

Figure 56

Figure 57

Through-thickness variation of the stress 022 at X1 = a/2 and X2 = a/2 in the
sandwich plate with aspect ratio 4. ................................................................ 85

Through-thickness variation of the stress 033 at X] = a/2 and X2 = a/2 in the
sandwich plate with aspect ratio 4. ................................................................ 86

Through-thickness variation of the stress 023 at X] = a/2 and X2 = O in the
sandwich plate with aspect ratio 4. ................................................................ 87

Through-thickness variation of the stress 0'13 at X] = O and X2 = a/2 in the
sandwich plate with aspect ratio 4. ................................................................ 88

Through-thickness variation of the stress 0'12 at x1 = 0 and X2 = O in the
sandwich plate with aspect ratio 4. ................................................................ 89

Through-thickness variation of the displacement u. at X] = O and X2 = a/2 in
the 4-layer antisymmetric cross-ply laminate. ............................................... 91

Through-thickness variation of the displacement u2 at x; = a/2 and X2 = O in
the 4-layer antisymmetric cross—ply laminate. ............................................... 92

Through-thickness variation of the stress on at X] = a/2 and x2 = a/2 in the 4-
layer antisymmetric cross-ply laminate. ........................................................ 93

Through-thickness variation of the stress 022 at x1 = a/2 and X2 = a/2 in the 4-
layer antisymmetric cross-ply laminate. ........................................................ 94

Through-thickness variation of the stress 023 at X] = a/2 and x2 = 0 in the 4-
layer antisymmetric cross-ply laminate. ........................................................ 95

Through-thickness variation of the stress 013 at x1 = O and X2 = a/2 in the 4-
layer antisymmetric cross-ply laminate. ........................................................ 96

Through-thickness variation of the stress 012 at X] = 0 and X2 = 0 in the 4-
layer antisymmetric cross-ply laminate. ........................................................ 97

xiii

CHAPTER 1 INTRODUCTION

1.1 Preliminary Information

The technology of composite materials is quickly developing to satisfy the
demand for materials with advanced properties. A typical laminated composite sandwich
structure allows the possibility to achieve new multi-functional such as high bending
stiffness, low weight, variable heat transfer capacities, sound insulation and extended
operational life. The sandwich is one of the composite structures that better expresses
these qualities. It consists of one or more high strength, stiff face layers separated by a
thick low-density flexible core. Whereas the facings provide the primary inplane load
capacity, the core carries the majority of the transverse load.

The increasing use of these materials requires a more accurate analysis of stresses
induced by thermal as well as mechanical loads. In particular, thermal stresses are
achieving a growing consideration in structural design since they come into play in
numerous situations, including the manufacturing process. Several components of the
aerospace industry, such as engine parts, heat shields, nozzles and nose cones are
subjected to repeated and prolonged exposure to severe heat environments. The rising

production of electronic devices, such as circuit boards, requires knowledge of thermal

stress behavior of composite materials to improve reliability. Many other examples from
the automotive and marine industry could be reported.

The accuracy in predicting thermal stresses in composite materials mainly
depends upon the correctness of two basic steps: the evaluation of temperature
distributions (heat transfer model) and the computation of the resulting stresses (plate
theories for structural analysis). To follow industry needs for innovation, precise
calculations should be maintained also in case of laminates composed of layers with

highly dissimilar material properties.

1.2 Literature Review

The development of heat transfer models for laminated structures has followed
closely that of laminate theories for structural analysis. The assumed through-thickness
distribution of temperature in thermal theories resembles the through-thickness
distribution of displacements in structural laminate theories. It follows that if the
structural model has a linear through-thickness distribution of displacements, then the
temperature distribution is computed by a thermal model based on a linear through-

thickness variation.

1.2.1 Equivalent Single Layer Approach

The most common approach is the equivalent single layer approach, wherein the
material properties of all the layers are "smeared", and the laminate is modeled as an
equivalent single anisotropic layer. The most popular of these theories is the first order

shear deformation theory (F SDT) [1, 2]. This theory yields good predictions of overall

laminate behavior and inplane stresses provided the plate is thin and the material
properties of adjacent layers do not differ significantly. However, F SDT does not account
for warpage of the cross section, which may be significant in laminated composites. In
order to reduce the inaccuracies of the F SDT, higher order shear deformation theories
(HSDT) were proposed [3-5]. These theories permit nonlinear variations of
displacements, strains and stresses through the thickness, and thus warpage of the cross
section. However, even though they often can predict well the gross behavior of thin and
some moderately thick laminates, all equivalent single layer theories have a common
shortcoming. They are unable to account for the sometimes-severe discontinuities in
transverse shear strains that occur at interfaces between two adjacent layers with
drastically different stiffness properties. In these cases, the local deformations and
stresses, and even the overall laminate response, are not well predicted.

The heat transfer models derived from the equivalent single layer approach share
the limits of the displacement theory. They assume a polynomial (usually linear)
variation of temperature through the thickness of the laminate and replace the discrete-
layer thermal pr0perties with an equivalent set of homogeneous anisotropic thermal
properties using a Thermal Lamination Theory (TLT) (see, for example [6-11]). For
greater accuracy, higher order thermal lamination theories were proposed [4,6,12].
However, when (i) the span-to-thickness ratio of the laminate is small, (ii) the thickness-
direction components of the thermal conductivity of adjacent plies differs significantly,
and (iii) the time gradients of temperature are large, both the first and high-order TLT

produce approximate results.

ac:

of:

1.1

dis
are
cor
err
sir:

ten

inc

1.1

ant

1.2.2 FSDT with Postprocessing Techniques

In recent years, since finite element models based on the first order shear
deformation theory are not computationally expensive, their use in combination with
postprocessing techniques has received increased attention. These approaches allow more
accurate evaluations of transverse shear stresses in thermally loaded laminates by the use
of: 1) three-dimensional equilibrium equations [13,14], 2) predictor-corrector approaches

and 3) simplifying assumptions.

1.2.3 Discrete-Layer Theories

In an effort to overcome the deficiencies of the equivalent single layer approach,
discrete-layer (or layerwise) theories have been proposed (e.g., [IS-18]). These theories
are based on a unique displacement field for each layer and enforce interlarninar
continuity of displacements and sometimes of transverse stresses, as well. They predict
excellent global and local distributions of in-plane and out-of-plane displacements and
stresses. The related thermal model assumes a piecewise, or layerwise, variation of the
temperature, generating accurate temperature distributions [19]. The major shortcoming
in these theories is their large computational expense, for as the number of layers

increases, the number of degrees of freedom increases proportionally.

1.2.4 Linear Flux Theory

This approach employs the assumption that the fluxes vary linearly as the
dependent variables over the element [20]. Finite element equations derived for thermal
and structural analysis yield finite element in integral form that are different from those

appearing in conventional formulations. The advantage of this method consists in a

 

reduc

first c

1.2.5

“asc
lanfin
finch
numt
stress

very .

resea
impri
dISpl
Great

bend

dtfle

gene
Orde;
trPins
00nd

ElEm

reduced computational time, but the accuracy of the result is comparable to that of the

first order shear deformation theories.

1.2.5 Zig-Zag Laminate Theory

A new class of laminate theory, called here the first order zig-zag theory (FZZT),
was developed by DiSciuva in the mid 1980's [21,22]. The in-plane displacements in a
laminate are assumed to be piecewise (layerwise) linear and continuous through the
thickness, yet the total number of degrees of freedom is only five (not dependent on the
number of layers). This is accomplished by analytically satisfying the transverse shear
stress continuity conditions at each interface in the laminate. This theory was shown to be
very accurate for many cases, especially symmetric laminates.

On the basis of the concept introduced in [21,22], DiSciuva as well as other
researchers made significant improvements to the F ZZT [23-30]. The primary
improvement was achieved by superimposing a piecewise linear variation of in-plane
displacements on a continuous cubic function of the transverse coordinate [24-29],
creating a displacement field that can better account for the warping that occurs during
bending of asymmetric laminates.

However, these theories have the unfortunate lirrritation that the transverse
deflection degree of freedom w0 is required to be Cl continuous. Averill proposed a

generalized form of the first-order zig-zag theory [23] and a generalized form of the high-
order zig-zag theory [24] for beams to alleviate the requirement of Cl continuity on the
transverse deflection. A new variable was introduced along with an associated constraint
condition, which was enforced by employing an interdependent (anisoparametric)

element interpolation scheme and the penalty method.

5

 

Yip and Averill developed a model for laminated beams [25] and plates [26] that
combined the benefits of the discrete-layerwise and high-order zig-zag theories. This
model allows the laminate to be subdivided into a number of sublarninates, with each
sublarninate containing several, even many, physical layers. A zig-zag kinematic
assumption is used within a sublarninate. Each finite element represents one sublarninate,
and, if cast in the form of a four-noded quadrilateral (for beams) or an eight-noded brick
(for plates), refinement of a model by through-thickness discretization can be achieved
without the use of any special constraints. When only one sublarninate is used through
the entire thickness of the laminate, nodal degrees of freedom are three translations and
two rotations. Nevertheless, when multiple elements (sublarninates) are needed to model
a laminate, additional shear stress degrees of freedom are present, making the element
unsuitable for implementation into many of the current commercial finite element codes.

Cho and Averill [31], presented a refined plate theory based on sublarninate linear
zig-zag kinematics and a new 3-D finite element based on that theory. The new element
contains eight nodes of which each has only five engineering degrees of fieedom -- three
translations and two rotations. The element has shown to be accurate, simple to use, and

compatible with the requirements of commercial finite element codes.

1.3 Objective of the Present Study

The first task of this research project is to introduce a new thermal lamination
theory for evaluating the temperature distribution in composite plates and sandwich
panels. This theory applies to heat transfer problems the zig-zag sublarninate concept

described above. Accordingly, a layerwise through-thickness distribution of temperature

6

is assumed, and the interior DOFs are eliminated by analytically satisfying the condition
of transverse flux continuity at each ply interface. The finite element model is developed
in a sublarninate version that allows the laminate to be subdivided into a number of
sublarninates, with each sublarninate containing a number of physical layers.

The second objective of the work is to combine the presented thermal theory with
the structural theory developed in [31], in order to realize a finite element formulation for
accurate thermal stress analysis of composite laminates. With respect to others theories,
the presented formulation is expected to show superior results in predicting the
distribution of temperatures, displacements and stresses in thin or thick laminated plates

wherein the plies have dissimilar thermal and/or structural properties.

1.4 Organization of the Thesis

In Chapter Two, the basic concepts of the new zig-zag sublaminate thermal theory
are presented and the possible forms of the finite element models based on the theory are
investigated. A four-noded two-dimensional sublarninate finite element model is
developed and implemented in a Fortran code. Several numerical results are reported to
demonstrate the effectiveness of the theory.

In Chapter Three, to test the validity of the theory in the presence of material non-
linearity, the thermal conductivity coefficients of the individual layers are assumed to be
a quadratic ftmction of the temperature. This has required non-linear capabilities to be
implemented in the model and in a new F ortran code adopting the incremental Newton-
Raphson technique. Numerical outcomes are used to show the high degree of adherence

of the model to the physical phenomenon.

This
(lot

impl

The]

code

pro;

In Chapter Four, a description of the zig-zag structural plate theory is provided.
This task is accomplished neglecting some details related to the mathematical
development and focusing the attention on the necessary changes required by the
implementation of the thermal stress capabilities.

In Chapter Five, the method of the equivalent nodal force vector is explained.
Then, the solutions of numerous thermal stress analyses, obtained by use of a Fortran
code developed fi'om the model, are compared with the exact ones to show that the
proposed procedure is efficient and accurate.

Final conclusions are given in Chapter Six.

CHAPTER 2 THERMAL LAMINATION
THEORY

2.1 Formulation of the Theory

It is assumed that the laminate is composed of N T perfectly bonded layers. The
thickness and material thermal properties may vary arbitrarily from layer to layer. The
laminate can be modeled as M sublarninates, with each sublaminate containing NM

layers, where m is the sublarninate number. Mathematically, this is represented as:

NTzsz (1)

In order to facilitate the development of the theory, all expressions in the
following derivation pertain to the m-th sublaminate. The sublarninate number
designation m is omitted for brevity. The reference surface of the m-th sublarninate is
assumed to be the bottom surface of the sublarninate, as shown in Figure 1.

The balance equation for thermal conduction within a single ortlrotropic

homogeneous layer can be written as:

6T _
PCE=V¢1+Q (2)

where:

6T
‘1.- = “k.- '5; (3)

is the flux in the i-th direction (i=l,2,3), p is the density, 0 is the specific heat, k, is the
thermal conductivity in the i-th direction, and Q is the internal heat generated per unit

volume.

When the material principle direction is oriented at an angle 6 to the coordinate
axis x], the thermal conductivity coefficients kl and k2 are calculated as follows. If kL is
the thermal conductivity of the lamina in the direction of the fibers and k, is the
transversal one, then:

kl =kL -cos20+k, .sin26?

(4)

k2 = kL -sin20+kT-c0526

 

 

 

x,

Figure 1 Geometry and coordinate system used in the thermal sublarninate
theory.

The set of essential boundary conditions is of the form:

10

T = T on S 1 (5)
and the set of natural boundary conditions is of the form:

kl—a—Zfil+k2£fi2+k3§£fi3—cj+h(T—Tw)=0 on S, (6)
6 0x 6x

1 2 3
where h is the coefficient of convection and Tan is the temperature of the external fluid.
The conditions of temperature and normal flux continuity at interfaces between

adjacent plies are:

= T(k+l)

Isaak

at x3 = x3, k (k) 61*“)
3

TU‘)

 

 

13=x3k

aTum) (7)

 

_ (lt+l)
_ k3

 

 

3 3

Xr=xrr xfi’sr

for k = 1 to N-l, where N is the number of plies in the sublarninate. k3”) is the thermal

conductivity in the x3 -direction of the (k)'h ply.
The simplest zig-zag temperature distribution that can be assumed is a piecewise

linear one, where for the k-th layer:

k-l
T“) (x,,x2,x3,t) = T0 (x,,x2,t)+ x37”, (x,,x2,t)+ 2 (x3 — x3,.)T,.(x,,x2,t) (8)
i=1

In Eq. (8), 7}, is the temperature and T1 is the thickness gradient of the temperature at the

bottom of the sublarninate. The first set of interface conditions in (7) is identically

satisfied by the assumed temperature distribution in (8).

A

Following the example of the zig-zag structural theories, the DOFs T. can be

eliminated by enforcing continuity of the normal flux at ply interfaces, the second
condition in (7).
This can be performed analytically, resulting in the following expression:

11

T}=a.l] (9)

where:

 

kgi) i-l
a, =[k§‘”’ —1]- 1521a,. (10)
F

Now, the theory contains only the DOFs 7}, , 7; for each sublarninate. If the entire
thickness of a structure is modeled using a single sublarninate, then the present theory
resembles other commonly used thermal lamination theories, except that here a piecewise
linear (or zig-zag) distribution of temperature is assumed. Note that if the transverse
component of thermal conductivity is the same from ply to ply, then the parameters a,
are all zero, according to (10). In this case, the present theory reduces to one in which the
assumed through-thickness variation of temperature is linear through the entire
sublarninate.

When multiple sublarninates are used, it is necessary to enforce continuity of
temperature across the sublanrinate interfaces. This can be facilitated by replacing 7;, , I;
with the new DOFs Tb , I; , the temperature at the bottom and top surfaces, respectively,
of the sublarninate:

Tb =__T(l) , I; =T(N) (11)

x3=0 x3=h
where h is the thickness of the sublarninate.

The temperature distribution within a sublarninate is now written as:
T""=T,(1—<r>,)—T,cr>, (12)

where:

12

k-l
x3 + Z[x3 —x3,.]a,
" if]. (13)

- h + 21]}: -x3,.]a,.

 

(DI: (x3)

is a zig-zag through-thickness shape function that satisfies all interface conditions within
the sublarninate. Note that continuity of flux across sublarninate interfaces is not
enforced.

The heat transfer model described above is the simplest possible form of zig-zag
thermal lamination theory, and is apparently the first of its kind. It is not difficult to
envision higher-order forms of the model similar to the cubic zig-zag models used in

structural analysis [24,25,28].

2.1.1 Demonstration of Zig—Zag Effects

In this section a numerical example is presented in order to demonstrate the types

of temperature distributions that exist in laminates. Consider a plate having a unit
thickness h in the x3 direction and subjected to a uniform distribution of temperature

T =1 at the top surface, and T = 0 at the bottom surface. As a consequence of these
boundary conditions the value of the inplane dimension L has no effect on the results.
Further, the temperature distribution described in Eq. (12) represents the exact analytical
solution for this case. A three-ply laminate is examined with the top and bottom plies

consisting of the same material and the middle ply composed of a different one. The
material thermal conductivities in the x3 direction k3 are varied to demonstrate their

effect on laminate response. This variation is performed in such a way that the overall

laminate thermal conductivity, denoted 1?, , remains constant. In our case:

13

when

respe

ratio

them

 

abort
thick]
them
V'ariai
Iempi
repic

Conti:

norm

first

1amir

h 3
Efqaa=2@a no
0 isl

where k; and h, are the thermal conductivity in the x3 direction and the thickness,
respectively, of the i-th ply. The plies are of equal thickness, so h, = 1/3 for i = 1,2,3. The
ratio of thermal conductivity properties in the material of the first and third layer, k3I , and
the material of the second layer, k: , is defined as:

E
fl-E- on

With ,6 and K3 specified, the values of k; and k: are easily determined from the

above equations. In the present case, 1?, is taken to be 3.5, and ,6 is varied. The through-
thickness variation of the temperature is shown in Figure 2 for different values of the ply
thermal conductivity ratio ,6 . It can be that seen when ,6 =1, the through-thickness
variation of the temperature is linear. However, when ,6 is not equal to one, the
temperature field exhibits a distinct piece-wise continuous shape that cannot be
represented by even a high-order smooth polynomial. In these cases, a piece-wise

continuous firnction is required to capture the essential aspects of the temperature field.
In Table l, for the four values of ,6 considered in Figure 2, the temperatures,

normalized by the result that occurs when ,6 =1, are reported at the top surface of the
first and second layer. It is interesting to note that the temperature distribution in the

laminate depends very strongly upon the value of ,6.

l4

 

 

 

 

 

 

 

 

 

 

 

 

1.0
x3 / r"
/ ' l
0.9 4 I " I
i I ' '
/ .' |
0.8 . 5:1 / .. .
/ . l
/ ,' -
0.7 T fi_2 ,’ .' l
I / , , J
------ =10 I ,.-';'.”’
0.6+ /’ ,.-v,'.»'
/ . ‘ ’ '
/ / . 2‘ 3 T '
0 5 ,6—100 ‘5) .v
. i ’0’?-
I. " ‘ /
I ' .I' : . I
0.4- _;,. z’
a ' ’ ' ' ’
av o ‘ ' ’
a ’ 1 , ‘ ' . / I
0.3 . .r .' /
. /
l .. ,
I ' /
0.2 -' : /
. l I
l . /
I .0 l
0.1 4‘ ' /
I '0 /
I I / T
0'0 r r T l r T r 1 r
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Figure 2 Through-thickness variation of the temperature of a three-ply laminate

subjected to T =1 at the top surface and T = 0 at the bottom.

15

Table 1. Variation of temperature versus ply thermal conductivity ratio.

 

 

 

 

'6 l_(_fl_)_ at x3 =1/3 —T—(fl)— at x3 =2/3
T (B = 1) T (13 =1)
1 1 .0000000 1 .0000000
2 0.7500000 1 .1250000
10 0.2500000 1.3750000
100 0.02941 18 1 .4852941

 

 

 

 

 

 

2.2 Finite Element Model

2.2.1 Possible Forms

A finite element model based on the above theory could take many forms. In the
case where the entire laminate thickness is modeled using a single sublarninate, the

simplest elements would take the form of a three-noded triangle or four-noded
quadrilateral, with two DOFs per node, 7}, , 7; or T, , T, , as shown in Figure 3.

Alternatively, a three-dimensional finite element model of the forms shown in
Figure 4a can be developed with one DOF per node. This approach has the additional
advantage that multiple elements can be used through the thickness of a laminate, if
needed, to further increase the accuracy of the approximation (see Figure 4b). The model
can then be considered a multi-sublaminate model with one degree of freedom per node,
see Figure 5.

For the purpose of demonstrating and assessing the new thermal lamination
theory, a four-noded two-dimensional sublarninate finite element model has been

developed, as shown in Figure 6.

l6

 

/ / Q

Figure 3 Quadrilateral and triangular finite elements based on the zig-zag
thermal lamination theory. Each node contains two computational
DOFs.

4-/ 4...;
(a) ’

1.”
i

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

""+“‘\

 

(b)

Figure 4 (a) Eight-noded brick and six-noded wedge three-dimensional finite
elements based on the zig-zag thermal lamination theory. Each node
contains one computational DOF - the temperature at that node. (b) An
example of discretizing the thickness of a laminate using multiple
elements (sublarninates).

 

Figure 5 Element topology and nodal degrees of freedom for the three-
dimensional form of the thermal model.

x3

 

 

 

 

2 x1
—-—I*

 

Figure 6 Four-noded two-dimensional sublaminatefinite element model

18

2.2.2 Four-Noded Layered Rectangular Thermal Element
If the temperature does not vary in the x2 -direction, the balance equation for

thermal conduction takes the form:

"€11: 6T) a [k.aT]=pc-a§—Q (16)

ax ‘5; ’5,— 5}:

Ignoring transient effects and internal heat generation, the following equation governs our

problem:
__a_ .12: __6_ .32: -0 (,7,
6x, 6x, 6x3 fix,
The corresponding weak statement is:
aar 6T 55T 6T
k —— +k —— dxdxz 6Td 18
J'{ 1[6x1 6x1] 3( 6x3 6163]} 1 3 J14, ]s ( )
0' r‘
where:
6T 6T
qn =kr‘a—x—"r‘1'k3'é‘na (19)

l 3

is the flux normal to boundary 1“.
The following approximation scheme is employed within the four-noded

rectangular element:
(k) 4 k
T =2T1N1(xvx3) (20)
J'=1

where:

19

Ni =[1_(pk(x3)].\111(x1)
Ni =[l—d>,(x3):|-‘I’2(x,)
N: =<Dk (x3)-‘I’2(x,)
N: =<Dk(x3)r‘l’,(x,)

(21)

‘P, and ‘1’, are the Lagrangian linear interpolation functions adopted to approximate T

long xl :

Wick;- (l-é‘)

 

1 (22)
as): (1+5)
And (D, is defined in equation (13) to be:
k-l
x3 + 2h, - x3,]a,.
<1). (x3) = r). (23)
h + Z[h — x3i]a,.
i=1
The resulting finite element model is:
4
ZKin =9: (24)
j=l
where:
” 6N." aN" aN.‘ aN’f
K; = Z I k: ——'-—1— +k," _.___2 dx,dx, (25)
k I 6x, 6x, 6x3 6x3
" A.
and
Q: = I [q. N." ]ds (26)
re

The integration along the x1 direction is performed using Gauss-Legendre
quadrature. The integration along x3 direction utilizes Newton-Cotes quadrature using 3

20

base points per layer.

The final expression for calculating the coefficients K ,j is:

1:
K5 = 5.“fo la: [xfm <k>l-¢i[xf°"°" will??? 651*
kt

+—6— - JK,‘ [40, [2.1 (01' [2.]+ o;[x,T°P(k)]-<1>,{[x,T°P(k)]]-[2

 

06 06
area]

.1014: lair...)-‘Pi(¢.-.)+W.(é.=2)-‘I’j(5p=2)]'i 6x 6x

where:

, (l-Cbk) if i=l,2
(I), = , .
(I), If 1:3,4

<D£ =[(1—<1>,) if j = 1,2

 

(1),‘ 1f j = 3,4
(28)
. (k)-x1°"<k>+xf°"°" (k)
"' 2

J = value of the J acobian for the change
in coordinate system from x, to f

2.3 The Computer Program “T hermZD ”

In order to apply the developed formulation for solving thermal conduction
problems a computer program has been realized. The computer language chosen is
Fortran77, generally recognized fiom the scientific community as one of the most
appropriate for writing finite element code. The code is reported in Appendix A. Several

comments have been inserted in the code to facilitate the reader’s understanding. In

21

particular, a detailed heading section explains the function of the main variables.

M. The input data file is very simple to build, since some time wasting
operations like writing the Boolean matrix are automatically accomplished by the
program. The material properties and the height must be specified for each layer
composing the laminate.

M. The output file in its first part summarizes the characteristics of the
analysis: element data, mesh data, Boolean matrix, nodal boundary conditions, number
and property of each layers element by element. It provides a complete output of the
analysis avoiding the necessity of printing also the input data file. Moreover, it makes
possible to accurately verify of the correctness of the input file. Subsequently, the thermal
solution at each node in the mesh is written. The last section of the file provides the
temperatures element by element at each layer interface along the two vertical sides.
Since the temperatures inside a single layer vary linearly, from the given information the

solution in every point of the domain can be calculated.

2.4 The Computer Program “Sinload”

In order to verify the effectiveness of the model the solutions coming from the
code had to be compared to the exact solution. Then, another F ortran code, called
“Synload”, has been realized. It is based on the publication [32], where the problem of a
laminated plate subjected to a sinusoidal thermal load is solved exactly. The code is
reported in Appendix B. The computer language is Fortran77, like for the other codes.

Structure. In [32] are provided the solution equations for the differential equation

governing each layer. The program solves the system of equations deterrrrined by

22

imposing the boundary and interface conditions. The results are the unknown coefficients
of the solution equations. Several comments have been inserted in the code to facilitate
the reader’s understanding.

input. The input data file is simpler than in the previous case, since no
information regarding the meshing of the domain is required. What needs to be inserted

are: the material properties and the height of each layer, the boundary conditions and the

horizontal xl coordinates where the solution is desired.
Output. The output file initially reports the characteristics of the analysis. Then, it

supplies the temperatures at each layer interface for the required xl coordinates.

2.5 Numerical Results

In this section the model is used to solve linear heat conduction problems to
demonstrate its applicability for thermal analysis of layered composites. For convenience,
all values of temperature are normalized by a reference temperature so as to non-
dimensionalize the results. Further, all dimensions and coordinates are normalized by the
total thickness of the laminate, h.

We have studied three types of laminates: a symmetric 8-layer laminate, an
unsyrnmetric 6-layer laminate and a laminate similar to the Shuttle thermal protection
system. In each case results are presented for two different aspect ratios of the composite
plate: 4 and 10.

All computational experiments, except the last one, consider a heat conduction
problem for which an analytical solution can be computed by use of the program

“Sinload”. It is the steady-state analysis of a laminated plate subjected at the top and
23

bottom to a sinusoidal distribution of normalized temperature of the type:

 

. Ir-Jcl
T...(x.) -1:-sm[ L j
(29)

 

T.....(x.) = T. -sin[”'L"‘ )

On the two lateral sides, at x, = O and x, = L , a uniform temperature T = 0 is

imposed. The value of the normalized temperatures T, and 7}, have been taken equal to l

and 0.1 respectively. For the case of the symmetric 8-layer laminate the model with the

applied boundary conditions is shown in the figure 7.

_ M.tcf..'2
ran-1 25-1: Mitcfllll

IOI
T-0.0 T=0.3826 T=0.7l7l T=0.9238 T ’

l T=0.9238 T=0.7l7l T=0.3826 T=0.0

T = 0 =—_=——=== T g 0

0,0 T=0.0 T=0.03826 T=0.07l7l T=0.09238 T=0.l T=0.09238 T=0.07l7l T=0.03826 T-0.0

 

Figure 7 B C ’s for the symmetric 8-layer laminate.

For the Shuttle thermal system an additional case is considered where three
surfaces are insulated: the top and the two sides.

It should be noted that when a uniform distribution of temperature is imposed on
the top and on the bottom of the laminated plate, the current model yields the exact
solution, with just one element through the thickness. Hence, in the examples presented,

only the more complex load schemes are considered.

2.5.1 Symmetric 8 Layer Laminate

Aspect ratio l/h = 10. In this test we have considered a steady-state heat

 

24

conduction problem in an 8-layer symmetric laminate composed of two materials whose

thermal properties are shown in Table 2.

 

 

 

 

 

Table 2. Material properties used in the numerical analyses.
Material ItI — Conductivity in the k,-— Conductivity in the
xl direction, x3 direction,
"I ' m .
Material 1 10.0 5.0
Material 2 1.0 0.5

 

 

 

 

 

The lamination scheme is (1/2/1/2/2/1/2/1). For the comparison between the exact

solution and our finite element model, the temperature is reported along the thickness at
x, = 2.5 and x, = 5. Eight different uniform meshes are used: 8X8, 8X4, 8x2, 8X1, 4X8,
4X4, 4X2 and 4x1, where the two numbers indicate respectively the subdivisions in x1 and

x3 directions. For example, since the laminate contains 8 layers of equal thickness, it

follows that in a mesh with 4 divisions along the thickness each element includes 2
layers.

Normalized temperature distributions are shown in tabular form in Tables 3 and 4
and in graphic form in Figures 8 and 9.

A close agreement between the exact analytical solution and the new theory is
obtained with the 4 by 1 mesh, and convergence toward the exact solution is achieved as

the mesh is refined. The linear solution is provided for the sake of comparison.

25

Table 3.

Temperature distribution along the thickness at x; = 2.5 for the 8 layer

symmetric plate with aspect ratio of 10. (Temperatures are normalized
by the exact solution.)

 

Exact

8x8

8X4

8x2

8x1

4X8

4x4

4x2

4x1

 

0.000

0.070711

1 .00000

1 .00000

1 .00000

1 .00000

1 .00000

1 .00000

1 .00000

1.00000

 

0.125

0. 083926

0.99981

1.00135

1.00566

1.01486

0.99925

1.00086

1.00533

1.01487

 

0.250

0.21 7500

0.9993 1

0.99917

1.01751

1.05659

0.99717

0.99710

1.01610

1 .05 660

 

0.375

0.231230

0.99926

1 .00074

1.01630

1.05640

0.99716

0.99875

1.01485

1.05641

 

0.500

0.3 72420

0.99944

0.99936

0.99863

1.04426

0.99776

0.99772

0.99697

1.04428

 

0.625

0.514760

0.99953

1.00101

1.01847

1.03648

0.99815

0.99975

1.01782

1.03649

 

fiso

0.5298 70

0.99957

0.99947

1.01818

1.03422

0.99825

0.99823

1.01761

1.03423

 

0.875

0. 689930

0.99994

1.00149

1.00280

1.00392

0.99980

1.00142

1.00277

1.00393

 

1.000

 

 

0. 707107

 

1 .00000

 

1 .00000

 

0.99999

 

0.99999

 

1.00000

 

1.00000

 

1.00000

 

1.00000

 

Table 4.

Temperature distribution along the thickness at x 2=

5 for the 8 layer

symmetric plate with aspect ratio of 10. (Temperatures are normalized
by the exact solution.)

 

Exact

8x8

8x4

8x2

8x1

4x8

4x4

4x2

4x1

 

0.000

0. 10000

1 .00000

1 .00000

1 .00000

1 .00000

1 .00000

1 .00000

1 .00000

1.00000

 

0.125

0.11869

0.99983

1.00135

1.00566

1.01487

0.99924

1.00084

1.00532

1.01487

 

0.250

0.30759

0.99932

0.99922

1.01754

1.05660

0.99717

0.99714

1.01610

1.05660

 

0.375

0.32 701

0.99927

1 .00076

1.01632

1.05641

0.99716

0.99872

1.01484

1.05641

 

0.500

0.52668

0.99943

0.99937

0.99865

1.04428

0.99776

0.99772

0.99697

1.04428

 

0.625

0. 72 798

0.99953

1.00104

1.01849

1.03649

0.99815

0.99974

1.01782

1.03649

 

0.750

0. 74934

0.99957

0.99951

1.01821

1.03424

0.99827

0.99824

1.01762

1 .03424

 

0.875

0. 9 7571

0.99995

1.00151

1.00281

1.00393

0.99980

1.00141

1.00277

1.00393

 

1.000

 

 

1. 00000

 

1 .00000

 

1.00000

 

1 .00000

 

1 .00000

 

1.00000

 

1 .00000

 

1.00000

 

1.00000

 

26

 

 

 

 

 

 

     

 

 

 

 

 

 

 

 

 

1.0

x 3
0.9 .

Exact

0.8 _

_ . _ . 8X4
0'7 " ...... 8X2 0

/

0. 6 _ _ _ _ 8x1 . .J

_ . . _ Linear Theory
0.5 .
0.4 -
0.3 -
0.2 -
0.1 -

T
0.0 , T . , , , Y
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
Figure 8 Through-thickness distribution of temperature at x) = 2.5 for an 8-layer

symmetric laminate with aspect ratio of 1 0.

27

 

 

 

 

    

 

 

 

 

 

 

1.0
x 3
0.9 _
Exact

0.8 -

_ . _ . 8X4
0.7 - ...... 8X2

/

_ _ .. 8X1 . 'I
0.6 -

_ . . _. Linear Theory
0.5 1
0.4 -1
0.3 -
0.2 .
0.1 .

T
0.0 , T . , , T , , ,
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

 

 

 

Figure 9 Through-thickness distribution of temperature at x; = 5 for an 8-layer
symmetric laminate with aspect ratio of 1 0.

28

Aspect ratio l/h = 4. A smaller aspect ratio produces a larger gradient of

 

temperature inside the plate in the x, direction. So this case tests the limits of our
theoretical approach. Materials, lamination scheme and the imposed sinusoidal
distribution of temperature are the same as in the previous analysis; the length-to-
thickness aspect ratio is taken to be 4.

Results are shown in Figure 10 only at the center of the plate, since the level of
agreement does not change significantly at different values of the x, coordinate. In this
case the solution obtained with an 8 by 4 mesh follows almost exactly the analytical
solution. The solution based on the 8 by 2 keeps a good degree of accuracy, instead the 8

by l meshes is able to give only a reasonable approximation.

2.5.2 Unsymmetric 6 Layer Laminate

Here we consider an unsyrnmetric 6-layer laminate composed of the same two
materials as in the previous analysis. The applied sinusoidal distribution of temperature is
also the same, but the lamination scheme is (1/2/2/1/1/1).

For aspect ratios of 10 and 4, respectively, results are reported in Figures 11 and
12. An analysis of the presented graph confirms that the present approach works well for
unsymmetrical laminates. Some degree of approximation is registered for the 8X1 mesh

in the plate having an aspect ratio of 4.

29

 

 

 
  
  
   
   

 

 

 

 

 

 

 

 

1.0
x3
0.9 -
Exact

0.8 .. _ . - . 8x4

...... 8x2 .4’
0.7 - , . l

-.’ I

_. _ _ 8X1 . 4; I
0‘6 ‘ _.._LinearTheory 'I'lz’
0.5 -
0.4 .
0.3 -
0.2 .
0.1 _

T
0‘0 1 r r 1 T r r r r
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

 

 

 

Figure 10 Through-thickness distribution of temperature at x; = 2 for an 8-layer
symmetric laminate with aspect ratio of 4.

30

 

 

 

 

 

 

 

 

 

 

 

 

 

1.0
x 3
0.9 . O/
0.8 . Exact ‘ '
/
_ . _ . 8X3 '
0.7 .
...... 8x2 '/
I
0.6- ___8x1 '.’ ,
.’ /
.. - . ._ Linear Theory ,’ I
0.5 .
0.4 1
0.3 .
0.2 -
0.1 .
0.0 , , . , , , , .
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Figure 11 T hrough-thickness distribution of temperature at x 2 = 5 for a 6-layer

unsyrnmetric laminate with aspect ratio of] 0.

31

 

 

 

 

 

 

 

 

 

 

 

 

1.0
x 3

0.9 .
0.8 . Exact ‘

_ . _ . 8X3
0.7 4

...... 8X2 .
0.6 . —. _ _ 8X1 1

l

_ . . .. Linear Theory 1
0.5 . J
0.4 _
0.3 -
0.2 .
0.1 -
0’0 r r f r r a r r r

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

 

 

 

Figure 12 Through-thickness distribution of temperature at x 2 = 2 for a 6-layer
unsyrnmetric laminate with aspect ratio of 4.

32

2.5.3 Shuttle Thermal Protection System

The Space Shuttle thermal protection system is composed of five different
materials, whose thermo physical properties are assumed to be constant and are given in
Table 5 (adopted from [12]). This case tests the accuracy of the present model when there
is a large difference in the values of thermal conductivity in adjacent plies. The thickness
of the RSI layer is much greater than all the other layers, so we will focus on the results
for the last four layers (RTV / FELT / RTV / ALUMINIUM) in order to get a better view
of the quality of the solution. An image of the model, for the aspect ratio L/h = 4, is given

in Figure 13.

Table 5. Shuttle thermal protection system properties adopted from [12].

 

 

 

 

 

 

 

Material Conductivity, ;W—K Density, -I,%_ 0 Specific heat, Kg-K
O42 coating 0.05 1665.92 500
RSI 0.043 114.17 500
RTV 0.32 1409.62 1200
Felt 0.0286 96.1 1 1200
Aluminum 125.0 2851.28 950

 

 

 

 

 

 

The characteristics of this analysis are similar to those of the symmetric and
unsyrnmetric laminates previously studied, but the applied sinusoidal distribution of
temperature is inverted. In this way the higher temperature is applied at the bottom to the
042-coating layer, that is the one having the stronger resistance to the severe environment
produced by the friction with the atmosphere.

Results in Figures 14 and 15 confirm that for an aspect ratio of 10 the current

33

approach provides adequate accuracy with only one element through the thickness. For
an aspect ratio of 4, only a small decrease in accuracy is observed in Figures 16 and 17.

In the final test case, the top and the sides of the plate are insulated and a
sinusoidal distribution of temperature is applied on the bottom surface. For this type of
boundary conditions, the temperature on the top is unspecified. Because an exact
analytical solution is not available for the present problem, a reference solution is
obtained using the finite element program Ansys 5.3 [38]. The mesh adopted for the
analysis with Ansys is made of 16 elements along xI and 8 along x,. As seen in Figures
18 and 19, even one element through the thickness is able to accurately capture the

temperature distribution in this case.

Alumlnlum
Shuttle thermal protection system Rrv

 

042 coating

Figure 13 Shuttle thermal protection system.

34

 

1.0

 

0.9 .

0.8 .

0.7 -

0.6 .

 

 

0.5 .1

 

Exact

 

0.4. _._.8X3

0.3- ___8x1

0.2 1 _ . . _ Linear Theory

 

 

 

0.1 .

 

 

0'0 r T r r T r r l
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

 

 

 

 

Figure 14 Through-thickness distribution of temperature at x, = 5 for the Shuttle
thermal protection system with aspect ratio of 1 0.

35

 

 

 

 

 

 

 

    

 

 

 

 

 

 

 

 

1.00 ‘
x 3 \.
0.99 - \“
\\
.\
0.98 -
0.97 -
0.96 -
0.95 .
Exact
0.94 - —-—-3x3
_ _ _ 8x1
0.93 - _
.. . - _ Linear Theory \
. l T
0‘92 l r l r r r r :\
0.09 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18
Figure 15 T hrough-thickness distribution of temperature at x; = 5 for the Shuttle

thermal protection system with aspect ratio of 10. Zoom view showing
the last 4 layers.

36

 

1.0

 

0.9. x3
0.8 1
0.7 -
0.6 _
0.5 .

0.4 .

 

Exact

 

0.3.
_._.8X3

0.2. ———8X1

 

 

0.1 . — - . — Linear Theory
L J

 

 

T
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.0

 

 

 

 

Figure 16 Through-thickness distribution of temperature at x; = 2 for the Shuttle
thermal protection system with aspect ratio of 4.

37

 

 

 

 

 

 

 

 

 

 

  

 

 

 

 

 

1.00 ‘
x 3 \“
\
0.99 q ‘\
\ \
\ \ ‘
0.98 t
0.97 _
0.96 A
0.95 - t
\
Exact \
0.94 a
_ . _ . 8x3
\\
0.93 .4 ———8X1 “\\
\‘ \
_.._.LineaITheory ‘\ T
0'92 T f r I I I l
0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18
Figure 1 7 Through-thickness distribution of temperature at x; = 2 for the Shuttle

thermal protection system with aspect ratio of 4. Zoom view showing the
last 4 layers.

38

 

1.0

 

0.9 a

0.8 t

0.7 4

0.6 t

0.5 1

0.4 _

 

0-3 ~ [ Ansys 16x8

 

0.24 -'-'8X3

0.1 4

 

 

 

 

,3 T
0.0 , T -— —

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

 

 

 

 

 

Figure 18 Through-thickness distribution of temperature at x; = 5 for the Shuttle
thermal protection system with aspect ratio of 10. Test case with top and
sides of the plate insulated.

39

 

 

 

  

 

 

 

 

 

 

 

 

1.00
x 3 I
0.99 _ I
0.98 . K
0.97 .
0.96 -
0.95 .
Ansys 16x8

0.94 _

.. . .. . 8x3
0.93 - ------ 8Xz

_ _ _ 8x1 ‘
0'92 r I J I I I I I T

0.14 0.15 0.16 0.17 0.18 0.19 0.20 0.21 0.22

 

 

 

Figure 19 T hrough-thickness distribution of temperature at x; = 5 for the shuttle
thermal protection system with aspect ratio of 10. Test case with top and
sides of the plate insulated. Zoom view showing the last 4 layers.

40

CHAPTER 3 NONLINEAR THERMAL
MODEL

To test the validity of the theory in the presence of material non-linearity, the
thermal conductivity coefficients of the individual layers were assumed to be a quadratic
function of the temperature. This has required non-linear capabilities to be implemented
in the model. Among the solution techniques more adopted for solving non-linear

problems, the incremental Newton-Raphson has been chosen.

3.1 Basic Relations

The mathematical formulation of the presented thermal theory obliges to evaluate
the temperatures in all layers composing the element at each iteration.
The stiffness matrix, Eq. (23) is a function of the thermal conductivity coefficients

of the materials composing the laminate, because they appear directly in the expression.

But this is not the only dependence, since the N f terms are also function of (1),:
N? =Nf(<1>»‘1’n‘1’2) (30)

And (D, is directly related to k," , as it has been previously defined in Eq. (10).

In Figure 20 the four-noded two-dimensional sublarninate finite element model

and the points where we need to evaluate the temperature are shown.

41

 

 

 

 

 

Figure 20 Points of the sublaminate element where the temperature needs to be
evaluated

According to [7], a quadratic dependence of the coefficients of thermal

conductivity on the temperature has been assumed. The adopted function is:
k=k0+11-T+A,-T2 I (31)
where the material parameters k0 , ,1, and 2,, define the behavior of the thermal
conductivity coefficient. Since the materials composing a laminate are usually
orthotropic, we will have different values of this parameters in the two directions x, and

x,.

3.2 The Computer Program “Th ermNL ”

The program has been realized in F ortran77.
Structure. The non-linear part of the structure is a typical example of a Processor
subroutine for applying the Incremental Newton-Raphson method. However, since in the

mixed B.C.’s we increment the flux while in presence of essential B.C.’s we scale the

42

temperatures, two processor subroutines have been written for managing in a different
way the two cases.

Input. The input file requires more information than the linear one, among these
we have: number of load steps, maximum number of iterations, convergence tolerance
and all the parameters required to evaluate the thermal conductivity in each material for
the two directions.

Output. The output file in its first part summarize the characteristics of the
analysis: element data, non-linear analysis data, mesh data, Boolean matrix, nodal
boundary conditions, number and property of each layer element by element. Then, the
thermal solution at each increment for a chosen node is written. After the last increment
the temperature at each node is presented, and for all the layers included in an element is

produced the solution at the four comers and at the center.

3.3 Numerical Results

The same symmetric 8-layer laminate studied in the first chapter has been
analyzed. It has the following lamination scheme (1/2/1/2/2/1/2/1), is composed by two
different materials and has an aspect ratio 1:10. Two types of analysis have been made:
one in which only essential boundary conditions are specified, and another where the
boundary conditions are mixed. For both problems, two possible sub-cases have been
studied: in the first one the value of A and A, are positive, while in the second they are
negative. This has been done since the response of the model will be different for thermal
conductivity coefficients that grow with the temperature as opposed to decreasing ones.

In the first case, the dependence on the temperature for the two materials, as shown in

43

Figure 21 in a range going from O to 100 °C, has been assumed:

k1] =10+O.OOl-T+O.Ol-T2

. _ 2 (32)
k,3-5+0.001-T+0.01-T

k: =1+0.0001-T+O.001.T2

2 _ 2 (33)
k, —O.5+0.0001-T+0.001-T

where the upper number marks the material and the conductivity coefficients are given in

 

 

 

 

 

unitsof
m-°C
120 "
/.
/ I
//
//
/ I'
/ I'
/ 1'
//'
//'
kxI.I(T) /"
__ /
kx3.1(T) ’f
kx1.2(T) //1 I
"""" /
kx3.2(T) /
/ .
/ .1
/ .I
/ .’
/ .’
/ .I
/ /
/ ,x
/ r
/ 1'
// ."
we!” 1'"
o -------------- ¥ """""""" + t : a‘

Figure 21 Thermal conductivities variation with T (°C) for 111 and A; positive

44

For decreasing thermal conductivity coefficients the dependence on the

temperature for the two materials has been reduced.

It; =10—0.001-T--0.0004-T2
kl, =5—O.001~T-0.0004-T2

(34)

k: =1-—0.0001-T—0.00004-T2
k: = 0.5 - 0.0001 - T — 0.00004 - T2

(35)

The assumed parameters let us reach very small values of thermal conductivity as

the temperature approaches 100 °C, avoiding negative values that are not admissible (see

 

 

 

 

 

 

Figure22).
Io‘__”““'“-~...
\\\
\\
__ \
\
\
\
\
\
.. \
\
\
\
\
u \
\

\\
kxI.I(T)4 \\
kx3.1(T)
kx1.2(T) -~-~ ~ ‘ ~
[013.2(an \ "\

CI!- -\ \-
\.\‘
\s
.l. \ ‘
\.\‘
\§

_‘_ ------------------------------------------------------------------------- ‘
n l 1 . l

I ' I I 1
0 T 100

Figure 22 Thermal conductivities variation with T ("C) for 11 and It; negative.

45

The previous choice of the parameters comes out from several analyses that have
been nm in order to find values producing a high level of non-linearity. Being careful at
the same time to avoid problems that can converge only by the use of special techniques

like the arc-length method. These tests, however, worked on the values of the
parameters 11 and A, , because the constant term k0 has not been changed for comparison
reasons with the linear results.

Because an exact analytical solution is not available for the present problems, a
reference solution is obtained using the finite element program Ansys 5.3 [3 8] with an

8x8 mesh.

3.3.1 Mixed Boundary Conditions Analysis

As it is shown in the Figure 23, a uniform heat flux has been specified on the top
surface of a laminated plate while a uniform temperature is assigned to the bottom
surface. Lateral surfaces are insulated.

El‘-."tn‘.‘$.-ES res: Mltfl'llll

 

 

 

H (Fl “400 65
ea II! or 10.1
Insullted =——= lllflllfled
—_ ‘
__—— “—
’ Temperature-0

Figure 23 Model for mixed B. C. ’s analysis

3. 3. 1.] Growing Thermal Conductivity coefficients

Validation of final results. We will compare the final solutions presented from

 

AnsysS.3, with those produced by our finite element model. The specified heat flux is

400 3;.
m

46

By the use of Figures 24 and 25, the temperature is reported along the thickness at
x, =2.5 and x1 =5. Four different meshes are used in ThermNL: 8x4, 8x2 and 8x1,
where the two numbers indicate respectively the subdivisions in 1:1 and x3 directions. A
close agreement between the solutions from AnsysS.3 and the program is obtained with

all the meshes; only a small error can be notice at x1 = 2.5 for the 8x1 mesh.

 

L0

 

(19.

(18-

 

(l7-

(16.

05 -

 

 

 

(14.

03 -

(12.

OJ _

 

T

 

 

01)

 

 

 

0 10 20 30 40 50 60 70 80

 

Figure 24 Final through-thickness distribution of temperature at x; = 2.5 for
values of the conductivity coefficients k increasing with T.

47

 

L0

 

(L9.

(l8-

 

(l7.

(16.

05 .

 

 

 

(14-

03 -

(12-

OJ .

 

T

 

 

 

0'0 j I I I I I T
0 10 20 30 4O 50 60 70 80

 

 

 

Figure 25 Final through-thickness distribution of temperature at x; = 5 for values
of the conductivity coefficients k increasing with T.

48

Validation of the solution through the increments. In this validation test, the
solutions given by Ansys and ThermNL have been recorded at the end of each increment,
or load step. In particular, in Figure 26 are reported the solutions at the central node of
the top surface obtained using 100 increments. The linear solution allows us to see how
soon the non-linear ones diverge from the temperatures that we obtain if the thermal
conductivity coefficients are constant. Since the results from the 8x2 and 8x4 meshes
perfectly agree with the Ansys, only those from the 8x2 and 8x1 meshes are reported in
order to get a clearer graph.

We could have predicted the observed behavior, knowing that for a specified heat
flux, the temperature in all points of the mesh decrease as soon as the values of the
thermal conductivity coefficients increase. In our analysis, as the flux keeps grong with
the increments the values of the temperature in the mesh increase. As a consequence, the
thermal conductivity coefficients grow making the temperature at the end of the
following increments smaller than those in the linear solution.

The adherence among the solutions is excellent for all the meshes. The same test
has been done at the central point of the mesh. The quality of the solution is unchanged,

as Figure 27 demonstrates.

3. 3. 1.2 Decreasing Thermal Conductivity coefficients

This validation analysis is similar to the previous one, with the only change that
now we will apply the negative values of the thermal conductivity coefficients, defined in
(34) and (35). In order to maintain the solution in a comparable range of temperature like
the other analysis, the specified heat flux is reduced from 400 to 65.

In Figure 28 is reported the solutions at the central node of the top surface. In this

49

case the linear solution goes under the others, because the thermal conductivity
coefficients decrease while the temperatures in the mesh increase.

Again, the adherence among the solutions is excellent for all the meshes.

 

 

 

__ Ansys 8x8
. 8x2

........... 8x1

 

 

 

 

Iter.

 

 

 

 

 

l 11 21 31 41 51 61 71 81 91

 

Figure 26 Solution at every increments at x; = 5 and x3 = 1 for values of the
conductivity coefficient k increasing with T.

50

 

 

 

 

 

 

 

 

 

 

60

50 -l

40 -
.____ Ansys 8x8

30 d 9 8X2
........... 8x1

20 _ .
..___..Lmear

10 J

Iter.
0 I I T I I I if I I
1 11 21 31 41 51 61 71 81 91

 

 

 

Figure 27 Solution at every increments at x; = 5 and x3 = 0.5 for values of the
conductivity coefficient k increasing with T.

51

 

7O

 

60 - f
.____ Ansys 8x8 .55;

 

  
 
 
    

50 _ 0 8X2

........... 8x1

40-

 

Linear

 

 

 

30-

20-

10.

Iter.

 

 

O T Y I I T T T I l
1 11 21 31 41 51 61 71 81 91

 

 

 

 

Figure 28 Solution at every increments at x; = 5 and x3 = 1 for values of the
conductivity coefficient k decreasing with T.

52

3.3.2 Essential Boundary Conditions Analysis
In this problem the same laminated plate has been subjected at the top and bottom

to a sinusoidal distribution of temperature of the type used in the first chapter:

 

T...<x.> =T.-sin(”;“)
(36>

 

Tbottom(’xl) = 7b . Sin(fl-.Lxl )

On the two lateral sides, at x, = 0 and xI = L , a uniform temperature T = 0 is

imposed. The value of the constants T, and I}, have been taken equal to 100 and 10

respectively. In this analysis at each load step, the program will increment the values of
the temperature. They will start from zero and reach their specified value at the last
increment. The model with the applied boundary conditions is shown in the Figure 29.

_ M‘tef‘.l 2
IS“: “er-Mum M‘ttflfl I

10 l
T=0.0 T=38.26 T=7 l .71 T=92.38 T= 100 T=92.38 T-7l .7l T-3 8.26 T-0.0 ’

;.‘ . -

 

 

 

 

T=° = “0
0,0 T=0.0 T-3.826 T=7.l7l T=9.238 T=10 T=9.238 T-7.l7l T-3.826 T-0.0

Figure 29 Model for essential B.C. ’s analysis.

3.3.2.1 Growing Thermal Conductivity coefficients

As in the previous two tests, we compare the solution at the end of each
increment. Since the temperatures are specified on all the surfaces, here we can check the

solution only in the internal points. The central node of the mesh has been chosen; its
coordinates are x, = 5 and x3 = 0.5. In Figure 30 are reported the solutions at this node
for an analysis of 100 increments. The adherence among the solutions is excellent for the

53

N‘J

 

 

IlC

SC

61

“'1

8k

re;

8x8, 8x4 and 8x2 meshes. A small error is registered for the 8x1 mesh. The results fiom
the 8x8 and 8x4 meshes are not reported in order to get a clearer graph. The final
temperatures obtained along the thickness at x, = 5 are shown in Table 6 and in Figure
31. A close agreement between the two solutions is obtained even with the 8 by 1 mesh,

and convergence toward the solution from Ansys is obtained as the mesh is refined.

Table 6. Comparison among the temperature distributions along the thickness at
x1= 5 at the end of the nonlinear analysis.

x3 Ansys 8x8 8x8 8x4 8x2 8x1

 

3. 3.2.2 Decreasing Thermal Conductivity coeflicients

This validation analysis resembles the previous one, except that the values of the
nonlinear thermal conductivity coefficients are negative. In Figure 32 is reported the
solutions in the central node of the mesh. Again the adherence among the solutions is
excellent for the 8x4 (not shown in the graph) and 8x2 meshes, but an error of about
4.7% is registered for the 8x1 mesh. There is a reason why the solutions are less accurate

with respect to those in which the flux is specified. Here the gradient of the temperature

along the xl direction is greater and the thermal theory has more difficulty in

representing it with few elements.

54

 

 

 

 

 

 

 

 

 

 

 

 

 

80
7o . 4y
———“‘ Ansys 8X8 2.1"“.
60 . 2
8x2 *3"
----------- 8x1 2" /
50 a ‘0‘ ///
Linear 2' /
,.¢" /,/V
40 r 4.41” / //
. "0;. I,//
30 -
20 - /
‘3'. / /
+' //’
It //
10 ~< *eggt”l
"JVt/(I
o ,u/ Incr.
l l l 21 31 41 51 61 71 81 91
Figure 30 Solution at every increments at x; = 5 and x3 = 0.5 for values of the

conductivity coefl‘icient 1: increasing with T.

55

 

1.0

 

0.9 -

 

0.8
Ansys 8x8

 

0.7. ___8x4 ,

0.6 .

 

 

0.5 _

 

0.4 _

0.3 .

0.2 .

0.1 .

 

 

 

0.0

 

I I I I I I I I I

0 10 20 30 40 50 60 70 80 90 100

 

 

 

Figure 31 Final through-thickness distribution of temperature at x; = 5 for values
of the conductivity coefficient k increasing with T.

56

 

 

 

 

 

 

 

 

 

 

 

 

50 - T

__ Ansys 8x8
4o .

9 8X2 x/x

........... 8x1
30 . Lrnear .
20 4 ‘//‘o‘/ .
10 .

Iter.
O I I I I I I I T I I
1 11 21 31 41 51 61 71 81 91 101

 

Figure 32 Solution at every increments at x; = 5 and x3 = 0.5 for values of the
conductivity coefficient k decreasing with T.

57

CHAPTER 4 FIRST ORDER ZIG—ZAG PLATE
THEORY

4.1 Formulation of the Theory

In this section a brief description of the first order zig-zag plate theory will be
provided, neglecting details related to the mathematical development and focusing the
attention on the implementation of the thermal stress capabilities. For a complete
explanation of the theory and the numerical validation studies, the reader is referred to
[31].

As in the thermal theory the laminate will be modeled as M sublarninates, with

each of them consisting of N," perfectly bonded layers. Again the global X3 axis is taken
perpendicular to the inplane X1, X2 coordinate axes and has its origin at the bottom of

the laminate, while x3 is a local sublarninate coordinate in the direction of X3 with its

origin at the bottom of the sublarninate (see Fig. 1). In the following derivation all
expressions are related to the m-th sublarninate, but all indices indicating the sublaminate
number are omitted for brevity.

The constitutive relations for a three dimensional stress state in the k-th layer of

the m-th sublarninate in the presence of a thermal load can be written as:

58

luff)‘ "cffl Cg) Cf? 0 0 C3)“ [81(1) _ a1“) . ATW
k k k k k k k

Uiz) Cf.) Cl.) Cl.) 0 o cg; 8g) _a§ ), ATI)

Gig) C1?) C2) Cg) O O Cg? 8%) “(19) .ATU)

1= < 1 (37)
cg) o o 0 C3? Cl? 0 7:?
of? 0 0 0 Ci? Cl? 0 r1?

 

 

 

 

 

 

05;), _C.‘? Cl? C? 0 0 Ct? .r.‘?-a.‘?-AT"".

where: AT”) = 7:2“. — T”)

initial
For most practical laminated composites, the 13 coefficients of the material
stiffness matrix are obtained from 9 independent material constants and a transformation
law (see [39]).
In the case of the thermal expansion coefficients the change between the material

and the laminate coordinate system follows the usual transformation rule:

a, = aL cos2 15l+arsin2 0
a2 = a, sin2 6 + a, cos2 19
. (38)
a!2 = 2(aL -a,)srn6cos6

a3 = “1'
where a, and a, are, respectively, the thermal expansion coefficient in the directions

parallel and perpendicular to the fibers, and I9 is the angle between fiber direction and

laminate x1 axis.

The sublarninate displacement fields are initially assumed in the following form:

1-1
u:”(x,,x2,x3,t) = “b +x3 "yr +2053 ’x31)'§t

i=1
k-l

u§*)(xl,x2,x3,t)=vb+x3o‘I’2+Z(x3-x3,.)-77,. (39)

i=1

.;~(.,.,,.,,.)=w,.(._:hs.)...(%)

59

where h is the thickness of the sublarninate, ab and V, are the axial displacements in x1
and x2 directions, respectively, at x3 = 0 , ‘1’, and ‘1’, are rotations of the normal at

x3 = 0 , and W, and w, are the transverse deflections of the bottom and top surfaces,

respectively, of the m-th sublarninate. Thus, it is assumed that uf") and ué") vary in a

piecewise linear fashion through the thickness of a sublarninate and ug“ varies linearly

through the thickness. The parameters 5,. and 77,. in Eq. (39) are eliminated by enforcing
the continuity of transverse shear stress at each interface [21,22]. The shear stress
continuity conditions at the k-th interface can be expressed as:

05? =C§§*". 0.? =05?" (40)

Assuming that infinitesimal strain theory holds, the form of the assumed

displacement fields in Eq. (39) is such that the inplane displacements uf") and ugk’

(It)

contribute only constant terms to the transverse shear stresses, 0'1 3 k)

and <75, . For

consistency, then, the transverse displacement u?) should also contribute only constant

terms to these stress components. In order to achieve this consistency, terms in the
transverse shear strain expressions that involve the transverse displacement variables are

evaluated at the midplane of the sublarninate, thereby taking a thickness averaged value

of these terms. Substituting Eqs. (37, 39) into Eq. (40) and solving for 5, and 77,, it can

be seen that 4‘, and I], depend on the ratios of shear properties between adjacent layers

and the shear deformation in each sublarninate. If either of these quantities is small, then

discrete-layer effects will also be small.

The rotational variables ‘I‘, and ‘1’, in the displacement fields are now

60

eliminated by introducing the variables u, and v, , the inplane translations at the top

surface of the sublarninate [25,26], in order to expedite the development of versatile finite
element models. Thus, rather than describing the inplane displacement field by a
translation and rotation at one point, it can more conveniently be described here by the
translation at two points.

The displacement fields must be further manipulated in order to develop a C0

fimte element. Because the derrvatrves of transverse deflections, b , ' , ” and —’-
6xl 6x, 6x2 6x2

 

appear in the displacement field, their second derivatives will be present in the strain
energy functional. C1 continuity of wb and w, , is thus required. It is desirable to alleviate

such a requirement by introducing new rotational degrees of freedom as follows:

%=—62b1 %=-62v gig-=61» 291:6" (41)
x1 6x, 6x2 6x2

These relationships will be constrained in the strain energy via a combination of
the interdependent interpolation and the penalty method. The displacement fields can

now be written using indicial notation as follows:
uf") =5“, (b513, ug“ =17“, 1113;}, ug“ = w, (2,, a =1...4; p =1,2 42)
where the index ,6 is used to denote the top and bottom of the sublarninate with 1:
bottom and 2: top, and (13:3,)ng and 9,, are shape fimctions in the thickness direction
[31]. Unless noted otherwise, summation on repeated indices is implied. The variables

5a,, are:

17119:“13’ 172p=vfi, Y‘s/9:61.61 1745:1921, (43)
The transverse normal strain arising from the assumed displacement field (42) is

61

constant through the thickness of a sublarninate. This assumption may generate

substantial errors in composites that have a soft core or layer, especially when thermal

loads are applied. Consider a laminate composed of 3 layers, one of which is very

compliant. When it is under loading the uniform transverse normal strain distribution

through the thickness is predicted to be uniform, but the actual distribution is not. Thus, it

is desirable to improve the distribution of the transverse normal strain through the

thickness. The improvement can be achieved by instead assuming a constant transverse

normal stress, 5",, through the thickness of a sublarninate. 633 can be determined using

Reissner's mixed variational principle [40]. From Eq. (37), the transverse normal strain is

 

obtained as:
_ k 1 k k I k k 1:
53(3)- = k (5'53 ) —C1(3)£11)( ) ”($952) Ci6)71(2)+
C( )
33
1 k k k k 1:
CT? —(C,‘, ’a, + C§,’a, + C§,’a,, + C§3)a3)AT( ’

Where, the part of the strain due to the thermal load is denoted by:

_n.(k)_

k k k k k
533 Cis'i')a1 d3)a2+C§6)a12 +C9513)A AT( )

d:)(

633 is then determined analytically from the following relation:

0 = f: fun (531;) h‘ffiifljdx,

kal ’31“)
The constant transverse normal stress 6"33 is found to be:
—TL k
:fl(uafia Wfir¢g2 {11:20 Qp)_C§;).£33()
The newly defined transverse normal strain can be written as:

k —Mk —71.k
53(3): 533” 533”

62

(44)

(45)

(46)

(47)

(48)

where:

—M(k)

5,, = j}(u'afl,wfl,¢ff;,‘l’ffg,flfl) (49)

The firnctions fl and f; are defined in [31].
The equations of motion, the essential and natural boundary conditions, as well as
the displacement-based finite element model can be developed from Hamilton's principle.
Again, attention will be limited to a single sublarninate. The energy functional is

modified here to include the imposition of the constraints in Eq. (41) via the penalty

method. Hamilton's principle for the m-th sublarninate can be defined as follows:
5][U-K—W+P]dt=o (50)
‘0

U is the internal strain energy:

” 1
_ (k) (k) (k) (k) (k) (k) (k) (k) (k) (k) (k) (k)
U—Z [5'37” 8” +022 622 +033 .933 +612 712 +033 7,3 +023 723 :1de (51)
k=1
k

where V, is the volume of the k—th layer.

K is the kinetic energy:

K = N [lp<’°[(a§*>)’ 402;“)2 +(a§*>)’]dV, (52)

k=l V1 2

where p“) is the density of the k-th layer and a superposed dot indicates differentiation

with respect to time. Wis the work of external forces:

W = [wag d!) + [1,... d1“ (53)
Q l"

BI= concentrated forces; ’1 = 03.1.77 1.; a =1,2; i =1,2

P is the penalty constraint:

63

p =§ “2%., 193:4] {1.4, +%:4]2“do (54)

a fl=l
where 7 is the penalty parameter. As 7 is assigned successively larger values, the
constraints of Eq. (41) are satisfied more exactly in the least squares sense. In the internal
strain energy expression, Eq. (51), the transverse normal stress term requires special

attention. The constant transverse normal stress leads to an asymmetric stiffness matrix.

Thus, instead of the constant transverse normal stress 6%) , 0g" from the constitutive

equations should be employed. Moreover, the term 3;? to be substituted in (51) will

consist only of the part 53?”). The thermal component of the total strain 532“") will play

its role in the nodal force vector. Making the necessary substitutions into Eq. (50), taking
the variation, and integrating by parts, the equations of motion are obtained (see [31] for

details).

4.2 Finite Element Model

The novel features of the theory suggest developing the finite element model in
the form of an 8-noded brick element with five degrees of freedom per node, Figure 33.

This approach has been used with great success [21,22], and has the very
important advantage of allowing discretization in both inplane and through-thickness
directions without the need for any special multi-point constraints. Elements are
assembled in the standard way. Therefore, the element can be used to predict the global
response of a laminate using only one (or a small number) of elements through the

thickness, and local effects such as interlaminar stresses can be captured by refining the

64

mesh in the thickness direction near the interface(s) of interest. The capability of
discretization in the through-thickness direction also enables one to simulate the

separation of laminated composites due to delamination.

ut’vt’wt’glt’ 62!

ub’vb’wb’glb’ 621»

 

Figure 33 Element topology and nodal degrees of freedom of the structural model.

Inplane displacement and rotational degrees of freedom are approximated by the
bilinear Lagrange interpolation fimctions. For transverse deflection degrees of freedom,
an interdependent interpolation concept similar to that developed by Tessler and Hughes
[41] is utilized (see [26,31] for details).

The interdependent interpolation scheme alleviates the shear locking problem but
does not eliminate it totally. In order to prevent the shear locking phenomenon, the
transverse shear strain fields have been made field consistent using the approach

presented in [42] (see also [26,31]).

65

CHAPTER 5 THERMAL STRESS MODEL

5.1 Equivalent Nodal Force Vector Method

The nodal temperature vector is obtained from a thermal analysis adopting the

zig-zag thermal lamination theory. Then, the effect of the thermal strain is imposed

through an equivalent nodal force vector {P}. In order to evaluate {P} it is necessary to

compute at the elemental level {Pg} by the following relation: g

{41% I[B""]’[D‘“]{e’“’}dn (55>

k=l y.

where [B‘“] and [D‘”] are, respectively, the nodal strain-displacement relation matrix

Tm} is the thermal strain vector

and the elasticity matrix governing the layer k. {3

evaluated at the layer level by:

I k ‘
£1T1( ) [611(k) - ATWW

ag") AT“)

T“) (I) (k)

8 .

?(k)i=<a3 AT 1 (56)
723 O

712(k) O

 

 

 

 

k 7112(1) lag) , AT“);
The integration is accomplished numerically by Gauss quadrature, using one point

66

along x, and x2 directions and two along x3 in each layer. The zig-zag thermal

lamination theory is once again used to find the temperature at the integration points for
each layer in terms of the nodal temperatures (see Eqs. (12,13)).

It has been shown in the previous section, that the term 53?“) , defined in (45),
should be subtracted from the total strain 53‘," to obtain the corrected transverse strain
and to maintain consistency with the constant transverse normal stress 63,. This is
accomplished by modifying the term 3313”) in (56):

42.2“) = a?” ~AT—Esl“? (57)

The equivalent nodal force vector {P} will be added to the external forces applied

to get the total force vector. Then the system is solved and the displacement vector {u} is

found. The strain vector is calculated at the integration points for each layer by use of the

relation:
{a‘*’}=[B‘*’]{u“”} <ss>

Finally, the stress vector {W} for a single layer is defined by:
{0?}=[D‘*’]({6‘“}-{6"“}) <59>
The possibility of having inconsistent thermal stress evaluation, see [34,35], is

avoided by calculating the stress at the center of the element with respect to the x1 and x2

dimensions. Moreover, at this point stresses are super-convergent [43].

67

5.2 Numerical Results

A computer program, realized in F ortran77, has been used to apply the presented
thermal stress model to the analysis of laminate plates. The original structural code of
Cho [31] was modified for this purpose. It is not reported in an appendix due to its length.

To assess the effectiveness of the presented model in evaluating thermal stresses
and displacements, two different simply supported laminated plates were studied: a
sandwich and a four-layered antisymmetric cross-ply laminate ([0/90/0/90]).

In both cases the plates are subjected to a temperature rise at the external faces

that has a double sinusoidal variation in the x1 - x2 plane:

 

 

 

 

 

. N . \
‘Tiop(x19x2) : 7; 'Sin[ 7: x1 osin[ It x2 ,
L, 2 L2 )
\ \ (60)
Ibottom(x19x2)=];'5in[ fl .sin[ fl..x2
Ll ) L2 )

For the sandwich the value of the temperatures T, and 1}, have been taken equal
to -l and 1, respectively, yielding a temperature gradient through the thickness. The
antisymmetric laminate is loaded with a uniform temperature through the thickness

produced by the conditions 1; = Tb =1.

The plates are square L1 = L2 = a and are constrained by the following boundary
conditions:

u2=0, u3=0, 61:0 at x,=0,Ll

61
ul=0,u3=0,92=0 at x2=0,L2 ( )

These boundary conditions allow an exact solution to be obtained. In the case of

the sandwich panel the exact solution is provided by Cho and Park [36]. The exact

68

 

solutio

[13,14

lamina

follow

Where

coeffic

Of the
“here
directl
folio“

1ayers

are C0

the ex

 

 

5.2,]

first '

 

6X16

solution for the antisymmetric laminate has been calculated by Rolfes, Noor, and Sparr
[13,14].
In all graphs the x3 dimension is normalized by the total thickness of the

laminate, h. Moreover, displacements and stresses will be given according to the

following non-dimensional forms:

5,.=-——_‘Z’—, ".= " (62)

 

-6
where a, = fl— and E, =1066Pa are the reference thermal expansion and elastic

coefficients, and T =|Tb|=|1;|.

Taking advantage of the symmetry of load and boundary conditions, one quarter
of the plate has been analyzed with three different meshes: 8x8x2, 4x4x4, and 2x2x2.
Where the three numbers indicate respectively the subdivisions in x1, x2 and 1:3
directions. For example, given that the antisymmetric laminate contains 4 layers, it
follows that in a mesh with 2 divisions along the thickness each element includes 2

layers.
Since stresses are evaluated at the centroid of each element, the results presented

are computed at the center of the element that is nearest the point of the plate for which

the exact solution is given.

5.2.1 Sandwich Panel
The steady—state thermo elastic bending of a sandwich square plate is analyzed
first. The thickness of the soft core is equal to 0.6h , while that of each of the two stiff

external layers is 0.2h , where h is the thickness of the plate. The physical properties of

69

the two materials are the following:

Stifl' external layers (Unidirectional graphite-epoxy composite)
EL = ZOOGPa, E, = SGPa, G” = SGPa, GT, = 2.2 GPa, v” = 0.25, v" = 0.35,
10" 10'6 W W (63)

=—2 —, 0: =50 —, k =50—. k =0-5—
“L xx ’ XK ‘ K-m ’ K-m

Sofl core
E, = lGPa, E, = ZGPa, G” = 0.8GPa, 0,, = 3.7 GPa, V” = 025,1!" = 0.35,
10*S W (64)

a, :07. =30X-K—, k, =kL =50K—;

Where subscripts L and T stand for the directions parallel and perpendicular to the
fibers for the external layers, whereas I and T mark inplane and transverse directions for
the soft core. Results are presented for two different aspect ratios of the sandwich plate: 4

and 20.

5.2.1.1 Demonstration of Thermal Zig-Zag Effects

In subchapter 2.1.1 a numerical example was presented in order to demonstrate
the types of temperature distribution that exist in a laminate. It has been proved that the
often assumed linear variation of temperature through the thickness of the laminate can
produce considerable errors. To further underline this statement, the exact temperature

distribution (from [32]) through the thickness at the center of the studied sandwich plate
(% ,g), is compared with the solution produced by the zig-zag thermal theory and with

that obtained by a linear variation in Figure 34. A remarkable difference in the quality of
the results between the two theories can be observed. The errors in temperature

distribution produced by the assumption of a linear variation will have a negative effect

70

 

 

 

 

 

on tht

graph

assun

0.9
{0.8 _
0.7 .
0.6
0.5
0.4
1 0.3
0.2

0.1

 

0.0

 

F 1°qu

 

 

on the accuracy of stresses and displacements. To emphasize this problem, in all the
graphs presenting the solutions for the sandwich plate, the results obtained with an

assumed linear variation of temperature have been included.

 

1.0

 

 

Exact

 

0.9 .

0.3 - . TLT 8x8x1

 

 

0.7 . \ .
\ _ _ _ Linear

 

 

 

0.6 .

0.5 -

0.4 _

0.3 -

 

0.2 -

0.1-

 

 

 

0.0

 

 

-l

 

 

 

Figure 34 Through-thickness distribution of temperature at x, = x2 = o/2 in the
sandwich plate with aspect ratio 20.

71

5.2.1.2 Sandwich with Aspect Ratio a/h=20

Figures 35 shows the predicted through the thickness variation of the inplane

displacement ul at x1 = 0 and x2 = g , whereas Figure 36 illustrates the distribution of u2
at x1 =% and x2 = 0. The normal stresses 0'“ , an and 0,, in the center of the plate

(xl =3- and x2 =%) are plotted in Figures 37, 38 and 39, respectively. The transverse

shear stresses are represented in Figures 40 and 41 (023 at x1 =% and x2 = 0, a", at

x1 = 0 and x2 = 521-). The inplane shear stress 0'12 is evaluated at the corner of the plate

(x1: 0 and x2 = 0) and is shown in Figure 42.
A close agreement between the exact analytical solution and the new theory is

obtained with all meshes, and convergence toward the exact solution is obtained as the

mesh is refined. The only exception is represented by 033 , where the assumed constant
through-thickness variation is not sufficient to capture the true variation of 0'33. Because

0'33 is about 5 order of magnitude smaller than 0'11 and an in this problem, the current

predictions are considered to be satisfactory.

72

 

 

 

 

 

 

 

 

 

 

 

 

1.0 - ,

x 3 /'
0.9 . /'

Exact l ,o
I /
0.3 . _ _ _ 128el 8x8x2 /'
...... 64el 4x4x4
0.7 .
_ . _ . 8e12x2x2

0.6 -

l - . . _ Linear '

{ Temperature I
0.5 .
0.4 3
0.3 .
0.2 l /'

0"
0.1 . /'
I.
x u 1
0‘0 4 I I . I I .
-4 -3 -2 -l 0 l 2 3 4

 

 

 

Figure 35 T hrough-thickness variation of the displacement u 1 at x; = 0 and x2 =
W in the sandwich plate with aspect ratio 20.

73

 

 

 

 

 

  
  
   

 

 

 

 

 

 

 

 

 

 

1.0 ,’
x 3 /'
0.9 . /
Exact ."
0.8 - '
_ _ .. 128el 8x8x2
'/
0-7 — ...... 64el 4x4x4 ,'
8 12 2 2 ’O
0.6. "m 6 xx /'
_ . . .. Linear '
0.5 . Temperature
0.4 _
0.3 s
0.2 s
.’
0.1 - /'
/' u 2
0.0 I . . . . T
-6 4 -2 o 2 4 6
Figure 36 T hrough-thickness variation of the displacement u; at x; = a/2 and x2 =

0 in the sandwich plate with aspect ratio 20.

74

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.0
x 3
0.9 4
0.8 -
0.7 -
Exact l
0.6 - .
_ _ _ 128el 8x8x2 {I
0'5 . ...... 64el 4x4x4
i
04‘ —-_.8e12x2x2 I:
_ . . _ Linear '
0.3 _ Temperature
0.2 4
0.1 -
0'11
0'0 I I I I I I I I I
-500 -400 -300 -200 -100 0 100 200 300 400 500

 

 

Figure 37 T hrough-thickness variation of the stress on at x; = a/2 and x2 = a/2 in
the sandwich plate with aspect ratio 20.

75

 

 

1.0

 

 

 

 

 

 

 

 

 

 

 

 

 

x 3
0.9 _
0.8 s
V—' Exact .'
0.7 - I
_ _ _ 128el 8x8x2 .'
l
0.6 _ .
...... 64el 4x4x4 {l
0,5 _ _ . _ . 8e12x2x2 {
_. . . _. Linear '
0.4 - 1 Temperature 1
l
0.3 . .-
l
0.2 . l
0.1 -
0'22
0'0 I I I . fl. I I I I
-500 -400 -300 -200 -100 0 100 200 300 400 500

 

 

Figure 38

76

Through-thickness variation of the stress on at x; = a/2 and x2 = a/2 in
the sandwich plate with aspect ratio 20.

 

 

 

 

 

 

 

 

 

 

 

1.0
x 3 I I l
0.9 - t I I
I I i
I I .
0.8. .; l l |
l
: i : l
0.7 - : , I I
E i : I
0.6 . Exact ; I I
I I
. |
_ _ .. 128el 8x8x2 '
0.5 ,_ .._ . -__-,_ .1
...... 64el4x4x4 I I i
0.4 s l I | :
_-_.8e12x2x2 i I: .
. | l ;
0.3. _.._Linear ' II I
Temperature I : .
0.2- | l | ..
' I
l I i z
o I | '
0.1 .
! I I '
I l
0,0 . . 1 . 1 1 g . 633
-0.0075 -0.0050 -0.0025 0.0000 0.0025 0.0050

 

 

Figure 39 Through-thickness variation of the stress 0'33 at x; = a/2 and x2 = a/2 in
the sandwich plate with aspect ratio 20.

77

 

 

 
   

 

 

 

 

 

 

     

 

 

 

1.0 . . a
I
0.9 . |
I I
: 1 1
0.8 . ! .......... I. . . . , I
I I I I
' |
o 7 ' I i :
I I | 1 I '
Exact | I I
0.6 4 I I z :
_ _ _. 128el 8x8x2 I I I I
. 1 I :
0.5 - ...... 64el 4x4x4 ! I |
I I I I
_ . _ .8e12x2x2 . I E I
0.4 _ I I '
. I 1 : '
_ . . - Llnear I - I
Temperature I I I 1
0'3 . I 1 I I
I I I
I 1 : :
0.2 . 1 , ......... 1 ..... : 1
I . I I
I
I I I
0.1 .
I I I
I I
0'0 r T T I II 1 T i 623
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0

 

 

 

Figure 40 Through-thickness variation of the stress 023 at x; = a/2 and x2 = 0 in the
sandwich plate with aspect ratio 20.

78

 

 

 

 

 

 

 

 

 

 

 

1.0 , | . .
x 3 I I II
0.9 . | I I
I l '
I I |
0.8 - | , .I ........... L
I I I I
I | 1
0.7 . - I
I I |
. I Exact
I I I
I I |
I . I .
0.5 I I I ! ...... 64e14x4x4
i I ' _ . _ . 8e12x2x2
I I
0.4 - - 1 ,
I I I _ . . _ Lmear
I I I Temperature
0'3 “ I 1 I
I I
: 1 I
0.2 . I .. .I ........... 1..
1 I 1 :
I l I .
0-1 - | I I I
I I
I I I : 0'13
0‘0 l T T T l T I I i T I T
-9 -8 -7 -6 -5 -4 -3 -2 -l 0

 

 

 

Figure 41 Through-thickness variation of the stress 013 at x; = 0 and x2 = a/2 in the
sandwich plate with aspect ratio 20.

79

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.0
x 3

0.9 I
0.8 I
0-7 - Exact

_ _ _ 128el 8x8x2 II
0.6 _ .

...... 64el 4x4x4 I
0.5 _

— . _ . 8e] 2x2x2

I
0.4 4 —--—Linear .
Temperature
0.3 _
0.2 I
0.1 _
0'12
0.0 I I I .
-7 l 2 3 4 5 6 7

 

 

 

 

Figure 42 Through-thickness variation of the stress on at x, = 0 and x2 = 0 in the
sandwich plate with aspect ratio 20.

80

5. 2. 1.3 Sandwich with Aspect Ratio a/h=4

An aspect ratio of 4 was chosen in the present example in order to highlight the
ability of the present model to obtain accurate results for very thick laminates.

All displacements and stresses are shown at the same positions in the plate as for

the previous case. Figures 43 and 44 show the plot of the displacements uI and 14,.
Normal stresses 0'“ , 0'22 and 033 in the center of the plate are shown in Figures 45-47.
The shear stresses 0'23 , a” and on are represented in Figures 48-50.

In this case the highly non-linear displacements, especially ul , are not captured

exactly by use of the adopted meshes. As a consequence, the accuracy of the transverse

shear stresses 0,, and 0'1 3 show a small decrease. Prediction of normal stresses maintains

a high level of precision, with the only exception being 0'33 for the same reasons

expressed above.

81

 

 

 

 

 

 

 

 

 

 

 

1.0
x 3

0.9 .
0.8 -
0.7 I
0.6 -
0-5 ~ Exact

_ _ _ 128el 8x8x2
0.4 .

...... 64el 4x4x4
0.3 -

_ . _ . 8e12x2x2
0-2 « _ . . .. Linear

Temperature
0.] -
ll
1
0'0 I I I I T
-1.5 -l.0 -O.5 0.0 0.5 1.0 1.5

 

 

 

Figure 43 T hrough-thickness variation of the displacement u 1 at x, = 0 and x2 =
a/2 in the sandwich plate with aspect ratio 4.

82

 

 

 

 

 

 

 

 

   
  
 
 
  

 

 

 

 

1.0 ,, /-
I/ ,
x 3 /'/ /
0.9 -
. /
0.8 - '
Exact

0.7 _

_ _ _ 128el 8x8x2

/

0.6 -. ...... 64el 4x4x4 '

_ . .. . 8e12x2x2
0.5 .

_ - - .. Linear

Temperature
0.4 .
0.3 -
0.2 .I /
I / '
0.1 - / '
. x
' /
’ ' / lo, u 2
0.0 l T fill , . , I . I
-10 -8 -6 -4 -2 O 2 4 6 8 10
Figure 44 Through-thickness variation of the displacement u 2 at x; = a/2 and x2 =

0 in the sandwich plate with aspect ratio 4.

83

 

 

 

1.0

0.9 .

0.8 -

 

0.7 _

0.6 -

0.5 .
_ . _ - 8612X2X2

0.4 -

 

 

_ _ _ 128el 8x8x2

_ - . _ Linear

Exact

.64e14x4x4

 

 

Temperature

 

0.3 _

0.2 -

0.1 _

 

 

 

 

 

 

 

0.0 r , . I

-400 -300 -200 —100 O 100

Figure 45 Through-thickness variation of the stress on at x; = a0 and x2 = a/Z in
the sandwich plate with aspect ratio 4.

 

 

84

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.0
x 3
0.9 .
0.8 .
I
0'7 “ Exact I
0 6 _ _ _ 128618X8X2 I
...... 64el4x4x4 I
0.5 A
_ . _ . 8e12x2x2 .
I
0.4 - _.._Linear I
Temperature I
I
0.3 . I
I
0.2 .
0.1 .
0'22
0'0 I ' [I ' I I I I I I I
-500 -400 -300 -200 -100 0 100 200 300 400 500
Figure 46 Through-thickness variation of the stress 0'22 at x; = a/2 and x2 = a/2 in

the sandwich plate with aspect ratio 4.

85

 

 

 

 

 

 

 

 

 

 

 

 

1.0
x 3
0.9 .
0.8 -
Exact
0.7 . _ _ _ 128e18x8x2
_ - _ . 64el 4x4x4
0.6 .
...... 8el 2x2x2
0.5 I _.._Linear
Temperature
0.4 I
0.3 .
0.2 .
0.1 .
0'33
0.0 , . . a
-0.100 -0.050 0.000 0.050 0.100

 

 

Figure 47 Through-thickness variation of the stress 0'33 at x; = a/2 and x2 = a/2 in
the sandwich plate with aspect ratio 4.

86

 

 

 

 

 

 

 

 

 

 

 

1.0 . . l .
x 3 I I I I
0.9 . I I l
I I I
| I l
l I
0.8 . I ........... , I
I I I
0-7 . Exact I I I
i . - '
_ _ _ 128el 8x8x2 - ' : '
0-6« l l : !
...... 64el4x4x4 i I 5 I
0-5 * _ . _.8e12x2x2 I I 5 I
i ' : I
_ . . _ Linear ' I I '
.4 . . I
0 Temperature I I : I
I I ' :
0.3 . I I
I I
0.2 - | .- ....... 1...: l
l . I I
| : I
I I I l
0.1 . I . I l
I - I I

0.0 I | I IT I T 023

0 10 20 30

 

 

 

Figure 48 Through-thickness variation of the stress 0'23 atx, = a/2 and x2 = 0 in the
sandwich plate with aspect ratio 4.

87

 

 

 

 

 

 

 

 

 

 

 

1.0 . I .
x3 : : ;I
0.9 q I I .I
l : II
! I .'
0.88 I l'l‘ .....--:!
E : ' I
0.7 I . . |
' : I I
: 2 I I Exact
0.6 a I : I
i : : I _ _ _ 128e18x8x2
. I I
. . I .
0.5 a | g ! ...... 64e14x4x4
: : I
I ; : I _._ . 8e12x2x2
0.4 - : : I .
E i I | _. ._L1near
i : ' I Temperature
0.3 t : : : I
l I I
0.2 4 u L:- _ ....... I
| I :I
| I :I
0.1 a I I :I
I I
I I 'l 013
0.0 l T l f u I
-40 -30 _20 -10 0

 

 

 

Figure 49 T hrough-thickness variation of the stress 013 at x, = 0 and x2 = a/2 in the
sandwich plate with aspect ratio 4.

88

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.0
x , '

3 /
0.9 A
0.8 a
07 - ExaCt

_ _ _ 128el 8x8x2
0.6 e

...... 64614X4X4
0'5 ~ _. . _ - 8612X2X2

. l

_ . . .. Lmear

0.4 - T
emperature
0.3 a
0.2 J
0.1 a . ‘
' /
' .’ 0'12
0.0 I r 1 T I
-40 -3O 10 20 3O 40
Figure 50 Through-thickness variation of the stress on at x, = 0 and x2 = 0 in the

sandwich plate with aspect ratio 4.

89

 

5.2.2 Four-Layered Antisymmetric Cross-Ply Laminate
In this case the 4 layers have equal thickness, 0.25h , where h is the thickness of
the plate. The physical properties of the material are those typical of high-modulus

fibrous composites:

EEL :15, .91.; = 0.5, 911 = 0.3378, v” = 0.3, v" = 0.48,
T

E, E, (65)
ca = 0.139x10'6, a, = 9x 10"5

Where subscripts L and T stand for the directions parallel and perpendicular to the
fibers. The fibers of the top layer are parallel to the x1 axis. Results are presented for an
aspect ratio of 10.

For this problem the exact solution is presented in [13,14] only for the transverse
shear stresses 0'23 and 0,3. However, even though the others stresses and the inplane
displacements cannot be compared to an exact solution, they are reported here. All
displacements and stresses are evaluated at the same positions in the plate as for the
previous problem and presented in Figures 51-57.

The average values of the transverse shear stresses are well captured by the two

finer meshes, whereas the coarse 2x2x2 mesh produces slightly underestimated

predictions.

9O

 

 

 

 

 

 

 

 

 

 

 

 

1.0 .
I. I
x 3 ’1'. .
0.9 . I." .I
I; I
It '
0.8 . I: .1
II: .I
0.7 a I: I
_ _ _ 128el 8x8x2 [If I
0.6 _ ['0' I
It '
...... 64el4x4x4 ,; f
0.5 _ I,‘ /
I . o
I.‘ I
_._.8612X2X2 ,0 '
0.4 a o I
I: .
I; I
I - '
0.3 . I .. .I
I; I
I , n
0.2 . /,' .I
I,’
I,’ .I
0.1 - I: I
I; ,‘ u
I' o 1
0.0 . 1: 1 .
—4 -3 -2 -1
Figure 51 Through-thickness variation of the displacement u; at x, = 0 and x2 =
a/2 in the 4-layer antisymmetric cross-ply laminate.

91

 

 

 

 

 

 

 

 

 

1.0 ‘ .
\ s ‘
x3 \‘\ \‘
0.9 . \'. \
\z .
\; \
\ ~ '
0.8 . \.‘ \‘
O 7 \\“- \
' ‘ _-._128e18x8x2 \z -
' \
\1 .
0.6 - \-‘ \
...... 64e14x4x4 \ -. ‘\
\ ‘ I
0.5 - ; \'-. \
i _ . _.8e12x2x2 \: '\
0.4 - l V. ‘\
\‘. ~
\; \
\‘ '
0.3 . \‘. \‘
\\\‘ \‘
0.2 - \'. \
r. ‘
\‘. \-
0.l - \.| \
\- '
- \
\ “ t u 2
0.0 , \. 7 \
-4 -3 -2 -1

 

 

 

Figure 52 Through-thickness variation of the displacement u; at x, = a/Z and x2 =
0 in the 4-layer antisymmetric cross-ply laminate.

92

 

1.0

 

! _ _ .. 128el 8x8x2
0.9 . I '
I
I

0.8 1 ------ 64el 4x4x4

 

0.7 . _ . _ . 8e12x2x2

 

 

 

| l
I

|

I; .

0.6 . I: I
|.

| I

0.5 . l

0.4. |

0.3 . I

 

1
l
|
l— .
l

0.2 .

0.1 -

 

 

P--——_
DOIOIOOO

0‘0 0'1]

 

T

-20 -15 -10 -5 O 5 10 15 20

 

 

 

Figure 53 Through-thickness variation of the stress 0', 1 at x; = a/2 and x2 = a/2 in
the 4-layer antisymmetric cross-ply laminate.

93

 

1.0

 

0.9 .

0.8 .

0.7 _

0.6 -

0.5

0.4 -

0.3 .

0.2 _

0.1

 

0.0

 

_ _ _ 128el 8x8x2

...... 64el 4x4x4

——-———
Inna
—.—.—.—.-.

p:

 

5, _._.8e12x2x2

 

 

 

_ _I J

. 7.7:: :1 .-. n

in

 

.' l 0'22

 

_.—.—.- . -
—

 

-20

 

-15 ~10 -5 O 5 10 15 20

 

 

Figure 54

Through-thickness variation of the stress 0'22 at x; = a/2 and x2 = a/2 in
the 4-layer antisymmetric cross-ply laminate.

94

 

 

 

 

 

 

 

 

 

 

1.0 1
I :l
x ' 1'
0.9 . 3 I ::
. ;|
I :I
0.8 _I I ::
l ..:|
I I
0.7 . . I
I l
i l
0'6 ‘ ' . Exact
I : I
! : I _ _ _ 128el 8x8x2
0.5 . _l - .',_ _l
I I : ...... 64el4x4x4
I |
0'4 . i , _. _. . 8e12x2x2
. : |
I : |
0.3 . I z :
I 2.1
| I
0.2 a . I
I I
I I
0.1 . I !
' I
0.0 T 1 , I . 1 , . . 023
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

 

 

 

Figure 55 T hrough-thickness variation of the stress 0'23 at x; = a/Z and x2 = 0 in the
4-layer antisymmetric cross-ply laminate.

95

 

 

 

 

 

 

 

 

 

 

1.0 ,
I I
x3 I
I
0.9 - I
I .
I I
0.8 _ I I
I I
I- . I
' ' I
0.7 _ | I .
I : I
l I i
Exact l . .
0.6 . .
I . |
_ _ _ 128el 8x8x2 : : I
0.5 - _ _I, .: ,.. I
...... 64c] 4x4x4 I I l
I . I
I 1
0.4 . _._.8e12x2x2 I - I
I E I
0.3 . : I l
I, . .2 |
I' I
0.2 . II .
I: I
I: I
0.1 . I:
l-
l: I
0.0 , , , . . L' . l .
-0.8 -O.7 -0.6 -0.5 -O.4 -0.3 -O.2 -0.1 0.13 0.0

 

 

 

Figure 56 Through-thickness variation of the stress 073 at x, = 0 and x2 = a/2 in the
4-layer antisymmetric cross-ply laminate.

96

 

1.0

 

0.9 I

 

0.8 A — _ _ 12861 8x8x2

0.7 1 ...... 6461 4x4x4

0.6 -I _ . _ . 8612x2x2

 

 

 

0.5 H

0.4 I

0.3 A

0.2 A

0.14

 

I
I
I
I
I
I
I
I
I
l
I
I
I
I
I
I
I
I
I
1

_—-—_=.—-—-—_———————-__— —-_—_-

 

0'12
0.0

-l.0 -0.5 0.0

 

 

 

 

Figure 57 ”rough-thickness variation of the stress 0'12 at x; = 0 and x2 = 0 in the
4-layer antisymmetric cross-ply laminate.

97

CHAPTER 6 CONCLUSIONS

An improved thermal lamination theory has been presented that is both accurate
and effective for predicting the distribution of temperature in laminated structures
wherein the plies have dissimilar thermal properties. The theory lends itself well to
several forms of finite element models, including two-dimensional shell-type models and
three-dimensional solid-type models. In the latter case, a sublarninate form of the theory
is utilized that allows the thickness of a laminate to be discretized, if necessary, to
increase solution accuracy. The four-noded two-dimensional sublarninate finite element
model used to test the thermal theory behaved exceptionally well even in the presence of
material non-linearity, wherein thermal conductivity coefficients of the individual layers
are defined as quadratic functions of the temperature.

Starting fi'om the proposed thermal lamination theory and a zig-zag sublarninate
plate theory [31] that has been previously proven to be very effective, accurate thermal
stress analysis of laminated structures can be performed. With respect to others theories,
the novel features of the formulation allow superior results in predicting the distribution
of temperatures, displacements and stresses in laminated plates wherein the plies have
dissimilar thermal and/or structural properties. The three-dimensional solid-type models
permit a sublarninate form of the theory in order to increase solution accuracy by refining
the mesh through the thickness when needed. However, the accuracy of the first-order

zig-zag sublarninate kinematics approximation minimizes the need for multiple

98

sublarninates for many problems of practical interest. The finite element model is
"regenerated" in the form of an eight-noded three-dimensional element with one DOF
(temperature) per node for the thermal analysis and five (three translations and two
rotations) for the structural one. Thus, it is suitable for implementation into commercial
finite element codes. Numerical results demonstrate that the current model is accurate,
efficient, and robust for analysis of a wide variety of thick or thin laminated plates and
sandwich panels.

The presented theory could be applied to other fields of engineering. For example,
models for moisture absorption would take the same form as above, with temperature
being replaced by moisture content. For the case of electric fields, Maxwell’s equations
governing electric fields are similar in form to those of heat conduction, and so a model
essentially equivalent to that given above could readily be developed in terms of the

electric scalar potential.

99

APPENDICES

100

APPENDIX A FORTRAN CODE, THERMZD

101

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C

THERMZD

C
C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

OOOOOOOOOOOOOOOOOOOO

OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO

>
>
>

TITLE
PRTSTF

I N P U T

DOM_X, DOM_Z, NELM_X, NELM_Z

IF (IMESH =

DO LOOP:

> NODE(I,J),

0) THEN

I = 1 TO NELMSH

J:

END OF LOOP

> X(I),
ELSE

DO LOOP:

> DX(I)

Z(I) I

1 TO NODELM

1 TO NODMSH

I = 1 TO NELM_X

END OF LOOP

DO LOOP: I = 1 TO NELM_Z
> DZ(I)
END OF LOOP
END IF
> NSGD NSGF
> ISGD(I,J), J = 1 TO 2, I = 1 TO NSGD
> VSGD(I), I = 1 TO NSGD
> ISGF(I,J), J = 1 TO 2, I = 1 TO NSGF
> VSGF(I), I = 1 TO NSGF
DO LOOP I=1, NELMSH
NLAYER(I)
DO LOOP J=1, NLAYER(I)
KX(I,J), KZ(I,J), ZlLAY(I,J), ZZLAY(I,J)

END OF LOOP

END OF LOOP

PARAMETERS:

IUNIT -
OUNIT

BNDMAX -

EDFMAX -
NBCMAX -
NDFMAX -
NEMMAX -
NEQMAX -

NGPMAX -
NNMMAX -
NPEMAX -
NLAMAX -

INTEGER,
INTEGER,

INTEGER,
GSTIF IS
GSTIF IS
INTEGER,
INTEGER,
INTEGER,
INTEGER,
INTEGER,

NUMBER OF D.O.F.

INTEGER,
INTEGER,
INTEGER,
INTEGER,

CONTROL VARIABLES:

V A R I A B L E S

UNIT NUMBER OF INPUT FILE.
UNIT NUMBER OF OUTPUT FILE.

COLUMN DIMENSION OF GLOBAL STIFFNESS:
BANDED IF IMASS = 0,

SQUARE IF IMASS = 1.

MAXIMUM NUMBER OF D.O.F. PER ELEMENT.
MAXIMUM NUMBER OF BOUNDARY CONDITIONS.
MAXIMUM NUMBER OF D.O.F. PER NODE.
MAXIMUM NUMBER OF ELEMENTS IN THE MESH.
MAXIMUM NUMBER OF EQUATIONS (OR MAXIMUM
IN THE MESH).

MAXIMUM NUMBER OF GAUSS POINTS.

MAXIMUM NUMBER OF NODES IN THE MESH.
MAXIMUM NUMBER OF NODES PER ELEMENT.
MAXIMUM NUMBER OF LAYERS PER ELEMENT.

102

OOOOOOOOOOOOOOOOOOO000000000000

C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

OOOOOOOOOOOOOOOOOOOOOOO

PRTGRD - INTEGER, FLAG FOR PRINTING GRADIENTS:
PRTGRD = 0 FOR NO PRINT,
PRTGRD = 1 FOR PRINTING.

PRTSTF - INTEGER, FLAG FOR PRINTING STIFFNESS MATRICES:
PRTSTF = 0 FOR NO PRINT,
PRTSTF = 1 FOR PRINTING STIFFNESS OF ELEMENT ONE,
PRTSTF = 2 FOR PRINTING GLOBAL STIFFNESS.

ELEMENT PROPERTIES:

NDFELM - INTEGER, NUMBER OF D.O.F. PER ELEMENT.
NODELM - INTEGER, NUMBER OF NODES PER ELEMENT.
NLAYER - INTEGER, NUMBER OF LAYERS PER ELEMENT.

MESH:

BNDWTH - INTEGER, HALF BANDWIDTH OF GLOBAL STIFFNESS MATRIX
FOR BENDING PROBLEMS.

EX - DOUBLE PRECISION(NPEMAX), EX(I) IS THE X-COORDINATE
OF THE I-TH NODE OF THE ELEMENT UNDER CONSIDERATION.

IMESH - INTEGER, FLAG TO INDICATE HOW MESH WILL BE INPUT:
IMESH = 0 IF THE MESH IS INPUT BY HAND,
IMESH = 1 IF THE MESH IS COMPUTED BY THE PROGRAM.

NDFMSH - INTEGER, NUMBER OF D.O.F. IN THE MESH.

NDFNOD - INTEGER, NUMBER OF D.O.F. PER NODE.

NELMSH - INTEGER, NUMBER OF ELEMENTS IN THE MESH.

NODE - INTEGER(NEMMAX,NPEMAX), CONNECTIVITY MATRIX:
NODE(I,J) IS THE J-TH NODE OF THE I-TH ELEMENT.

NODMSH - INTEGER, NUMBER OF NODES IN THE MESH.

X - DOUBLE PRECISION(NNMMAX), X(I) IS THE GLOBAL
COORDINATE OF THE I-TH NODE.

NELM_X - INTEGER, NUMBER OF ELEMENTS IN THE MESH IN X DIR.

NELM_Z - INTEGER, NUMBER OF ELEMENTS IN THE MESH IN Z DIR.

DOM_X - DOUBLE PRECISION, DOMAIN IN X DIRECTION

DOM_Z - DOUBLE PRECISION, DOMAIN IN Z DIRECTION

DX - DOUBLE PRECISION, DIMENSION OF ELEMENTS
IN X DIRECTION STARTING FROM THE LEFT

DZ - DOUBLE PRECISION, DIMENSION OF ELEMENTS
IN Z DIRECTION STARTING FROM THE BOTTOM

KX - DOUBLE PRECISION, TERMAL CONDUCTIVITY COEFFICIENT OF
THE LAYER ALONG X DIRECTION

KZ - DOUBLE PRECISION, TERMAL CONDUCTIVITY COEFFICIENT OF

THE LAYER ALONG Z DIRECTION

ZlLAY - DOUBLE PRECISION, Z COORDINATE OF THE BOTTOM OF THE
LAYER

Z2LAY - DOUBLE PRECISION, Z COORDINATE OF THE TOP OF THE
LAYER

BOUNDARY CONDITIONS:

ISGD - INTEGER(NBCMAX,2), ISGD(I,J) INDICATES THAT THE J-TH
DISPLACEMENT D.O.F. OF THE I-TH NODE IS SPECIFIED.

ISGF - INTEGER(NBCMAX,2), ISGF(I,J) INDICATES THAT THE J-TH
FORCE COMPONENT OF THE I-TH NODE IS SPECIFIED.

NSGD - INTEGER, NUMBER OF DISPLACEMENT BOUNDARY CONDITIONS.

NSGF - INTEGER, NUMBER OF NODAL FORCE BOUNDARY CONDITIONS.

VSGD - DOUBLE PRECISION(NBCMAX), VSGD(I) IS THE SPECIFIED
VALUE OF THE I-TH DISPLACEMENT BOUNDARY CONDITION.

VSGF - DOUBLE PRECISION(NBCMAX), VSGF(I) IS THE SPECIFIED

VALUE OF THE I-TH NODAL FORCE BOUNDARY CONDITION.

000000000000000OO00000000000000000000000OOOOOOOOOOOOOOOOOOOOOOO

FINITE ELEMENT MODEL:

103

OOOOOOOCOCOOOOOOOOOOOOOOOOOOOOOOO000(")0OOOOOOOOOOOOOOOOOOOOOOOOO

DSF

EF
ESTIF
EU
GAUSPT
GAUSWT
GDSF
GF
GSTIF
NGPSTF

NGPGRD

SF

DOUBLE PRECISION(NPEMAX), DSF(I) IS THE LOCAL (XI)
DERIVATIVE OF THE I-TH SHAPE FUNCTION.

DOUBLE PRECISION(EDFMAX), ELEMENT FORCE VECTOR.
DOUBLE PRECISION(EDFMAX,EDFMAX), ELEMENT STIFFNESS.
DOUBLE PRECISION(EDFMAX), ELEMENT EXTENSION VECTOR.
DOUBLE PRECISION(NGPMAX,NGPMAX) GAUSPT(I,J) IS THE
I-TH GAUSS POINT OF J-TH GAUSS RULE.

DOUBLE PRECISION(NGPMAX,NGPMAX) GAUSWT(I,J) IS THE
I-TH GAUSS WEIGHT OF J-TH GAUSS RULE.

DOUBLE PRECISION(NPEMAX), GDSF(I) IS THE GLOBAL (X)
DERIVATIVE OF THE I-TH SHAPE FUNCTION.

DOUBLE PRECISION(NEQMAX), GLOBAL FORCE VECTOR,
WHICH, ON RETURN, CONTAINS THE SOLUTION VECTOR;
ALSO USED AS A DUMMY VECTOR IN EIGENPROBLEMS.
DOUBLE PRECISION(NEQMAX,BNDMAX), GLOBAL STIFFNESS.
INTEGER, NUMBER OF GAUSS POINTS TO BE USED TO
INTEGRATE THE STIFFNESS MATRIX.

INTEGER, NUMBER OF GAUSS POINTS AT WHICH THE
GRADIENTS ARE TO BE EVALUATED.

INTEGER(NPEMAX), SF(I) IS THE I-TH SHAPE FUNCTION.

IMPLICIT LOGICAL (A-Z)

INTEGER IUNIT , OUNIT , NEMMAX, NPEMAX, NDFMAX, NGPMAX

NNMMAX, NBCMAX, NEQMAX, EDFMAX, BNDMAX,
NELM_X, NELM_Z, NLAMAX

PARAMETER (IUNIT = 10, OUNIT = 12, NEMMAX = 400,
NPEMAX = 4, NDFMAX = 1, NGPMAX = 5,
BNDMAX = 400, NLAMAX = 40,

EDFMAX = NDFMAX * NPEMAX,

NNMMAX = NEMMAX * (NPEMAX - 1) + 1,
NEQMAX = NNMMAX * NDFMAX,

NBCMAX = NEMMAX*4)

INTEGER PRTSTF, PRTGRD, nodelm, ndfelm, ngpstf, ngpgrd
nelmsh, ndfmsh, nodmsh, ndfnod, node , bndwth
nsgd , nsgf , isgd , isgf,

NLAYER
DOUBLE PRECISION X , vsgd, vsgf , ex , DOM_X, DOM_Z,
GAUSPT, GAUSWT, ESTIF , EF , GSTIF ,
GF , DSF , GDSF ,
Z , KX , KZ, ZlLAY, ZZLAY, EZ,
AK , FI_D , D_FI_N, D_FI, SFZ, GUL
DIMENSION NODE(NEMMAX,NPEMAX), ISGD(NBCMAX,2),

ISGF(NBCMAX,2), X(NNMMAX) , EX(NPEMAX) ,
VSGD(NBCMAX) , VSGF(NBCMAX) ,
EF(EDFMAX) , GF(NEQMAX) ,
DSF(NPEMAX) , GDSF(NPEMAX) ,

GAUSPT(NGPMAX,NGPMAX), GAUSWT(NGPMAX,NGPMAX),
ESTIF(EDFMAX,EDFMAX) ,

GSTIF(NEQMAX,BNDMAX),

Z(NNMMAX), NLAYER(NEMMAX),

KX(nemmax,NLAMAX), KZ(nemmax,NLAMAX),

104

OOOOOOOOOOOOOOOOOOOOOO

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

I

I

I

ZlLAY(nemmax,NLAMAX),

ZZLAY(nemmax,NLAMAX),

EZ(NPEMAX), AK(nemmax,NLAMAX), FI_D(nemmax),
D_FI_N(nemmax,NLAMAX), D_FI(nemmax,NLAMAX),
SF2(2,2),GUL(nemmax,2,NLAMAX)

open(unit=10,file='TERM2D.data',status='old')

open(unit=12,file='TERM2D.out',

INPUT OPTIONS.

CALL INPOPT (IUNIT , OUNIT ,

INPUT ELEMENT DATA.

call INPELM (OUNIT ,
NDFNOD)

call INPMSH (IUNIT ,
NDFNOD,
NELMSH,

OUNIT ,

NEMMAX,

INPUT MESH PARAMETERS AND FORM THE MESH.

NPEMAX,
DOM_X, DOM_Z, NELM_X, NELM_Z,
NDFMSH, NODMSH, NODE, BNDWTH, x, z )

status='unknown')

PRTSTF, PRTGRD)

NODELM, NDFELM, NGPSTF, NGPGRD,

NNMMAX, NODELM,

--- INPUT NODAL BOUNDARY CONDITIONS.
call INPNBC (IUNIT , OUNIT , NBCMAX, NNMMAX, NODMSH, NDFNOD,
NSGD , NSGF , ISGD , ISGF , X Z, VSGD,
VSGF )
--- INPUT ELEMENT LAYERS
call INPLAY (IUNIT , OUNIT , NELMSH, NLAYER, NLAMAX, NEMMAX,
KX , KZ , ZlLAY , ZZLAY)

CALL GQUAD (NGPMAX,

(OUNIT ,
BNDMAX,
NDFMSH,
NSGD ,
GAUSPT,
ESTIF ,

NEMMAX,
NBCMAX,
NELMSH,
ISGD ,
GAUSWT,
EF ,

call Procsr_T

GAUSPT, GAUSWT)

NPEMAX,
EDFMAX,
NODE ,
NSGF ,
DSF ,
GSTIF ,

105

NGPMAX,
NDFNOD,
NGPSTF,
ISGF ,
GDSF ,
GF ,

NNMMAX,
NDFELM,
PRTSTF,
X I

VSGD ,

INITIALIZE THE ARRAYS WHICH CONTAIN THE GAUSS POINTS AND WEIGHTS.

NEQMAX,
NODELM,
BNDWTH,
EX ,

VSGF ,

z , EZ , KX , KZ, ZlLAY, ZZLAY,
NLAMAX, NLAYER, AK , FI_D, D_FI_N, D_FI, SF2 )

C _____________________
C --- POSTPROCESSOR ---
C _____________________
call PRTGU (OUNIT , NEQMAX, NNMMAX, NDFNOD, NODMSH, GF,
X, Z )
call PR_L_T (OUNIT, NPEMAX,NLAMAX,NELMSH, NEMMAX,
NLAYER, NEQMAX, ZlLAY, ZZLAY,
AK, FI_D, NODE, GUL, GF)
STOP
END
C ---------------------------------------------------------------------- C
C C
C SUBROUTINES C
C C
C ---------------------------------------------------------------------- C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: ASMBLE C
C PURPOSE: ASSEMBLE THE GLOBAL MATRICES IN BANDED FORM. C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE ASMBLE (NEQMAX, BNDMAX, EDFMAX, NEMMAX, NPEMAX, NODELM,
NDFNOD, NODE , ESTIF , EF , GSTIF , GF ,
ELMNUM)

IMPLICIT LOGICAL (A-Z)
INTEGER NEQMAX, BNDMAX, EDFMAX, NEMMAX, NPEMAX, NODELM,
NDFNOD, NODE , ELMNUM, INOD , JNOD , IDOF ,
JDOF , ROW , COL , COLl , EROW , ECOL
DOUBLE PRECISION ESTIF , EF , GSTIF , GF
DIMENSION NODE(NEMMAX,NPEMAX) , ESTIF(EDFMAX,EDFMAX), EF(EDFMAX),
GSTIF(NEQMAX,BNDMAX), GF(NEQMAX)
DO 40 INOD = l, NODELM
ROW = NDFNOD * (NODE(ELMNUM,INOD) - 1)
DO 30 IDOF = l, NDFNOD
ROW = ROW + 1

EROW = (INOD - 1) * NDFNOD + IDOF
GF(ROW) = GF(ROW) + EF(EROW)

106

10
20
30
40

DO 20 JNOD = l, NODELM
COLl = NDFNOD * (NODE(ELMNUM, JNOD) - 1)

DO 10 JDOF = 1, NDFNOD
COL = COLl - ROW + JDOF + l
ECOL = (JNOD - l) * NDFNOD + JDOF
IF (COL .GT. 0) THEN
GSTIF(ROW,COL) = GSTIF(ROW,COL) + ESTIF(EROW,ECOL)
END IF
CONTINUE
CONTINUE
CONTINUE
CONTINUE
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C

C
SUBROUTINE: CR_FI C
PURPOSE: EVALUATE FUNCTIONS FI AND ITS DERIVATIVES C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

30

SUBROUTINE CR_FI(NLAMAX, NELMSH, NEMMAX,
NLAYER, KZ, ZlLAY, ZZLAY,
AK, FI_D, D_FI_N, D_FI)

IMPLICIT LOGICAL (A-Z)

integer NLAMAX, NELMSH,NEMMAX,
NLAYER, I, K, L , N

double precision KZ, ZlLAY, ZZLAY, AK, FI_D,
H, R, D_FI_N, D_FI, SUM_AK

DIMENSION NLAYER(NEMMAX), KZ(nemmax,NLAMAX),
ZlLAY(nemmax,NLAMAX), ZZLAY(nemmax,NLAMAX),
AK(nemmax,NLAMAX), FI_D(nemmax), D_FI_N(nemmax,NLAMAX),
D_FI(nemmax,NLAMAX)

DO 100 I = l, NELMSH

if (NLAYER(I) .gt. 1) then
SUM_AK = 0
AK(I,1) = ( KZ(I,1)/KZ(I,2) — 1 )
SUM_AK = AK(I,1)

DO 30 K = 2, (NLAYER(I)-l)
AK(I,K) = ( KZ(I,K)/KZ(I,K+1) - 1 )*( 1 + SUM_AK)
SUM_AK = SUM_AK + AK(I,K)
CONTINUE

END IF

H = ( ZZLAY(I,NLAYER(I)) - ZlLAY(I,1) )

107

4O

50

60

100

FI_D(I) H
DO 40 R
FI_D(I)
CONTINUE

1,
= FI_D(I) +

D_FI_N(I,1) = 1

DO 50 L = 2, NLAYER(I)
D_FI_N(I,L)
CONTINUE

DO 60 N
D_FI(I,N)
CONTINUE

1, NLAYER(I)
D_FI_N(I,N)

CONTINUE

RETURN
END

D_FI_N(I,L-l) +

(NLAYER(I)-l)
( H -

(ZZLAY(I,R)- ZlLAY(I,l)) )*AK(I,R)

AK‘IrL-l)

/ FI_D(I)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C

C
SUBROUTINE: procsr_T C
PURPOSE: DRIVER ROUTINE TO SOLVE THE equations. C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE procsr_T (OUNIT , NEMMAX, NPEMAX, NGPMAX, NNMMAX,

NEQMAX,
BNDMAX, NBCMAX, EDFMAX, ndfnod, ndfelm, nodelm,
ndfmsh, nelmsh, node , ngpstf, prtstf, bndwth,
nsgd , isgd , nsgf , isgf , X , EX,
GAUSPT, GAUSWT, dsf , GDSF , ESTIF,
EF , GSTIF , gf , vsgd , VSGF ,
Z , EZ, KX , KZ, ZlLAY, ZZLAY,
NLAMAX, NLAYER,AK, FI_D, D_FI_N, D_FI,
SF2)
IMPLICIT LOGICAL (A—Z)
integer OUNIT , NEMMAX, NPEMAX, NGPMAX, NNMMAX, NEQMAX,
BNDMAX, NBCMAX, EDFMAX, ndfnod, ndfelm, nodelm,
ndfmsh, nelmsh, node , ngpstf, prtstf, bndwth,
nsgd , isgd , nsgf , isgf ,
ELMNUM, ENOD , GNOD , IRES , I , J,
NLAMAX, NLAYER
double precision X , EX , GAUSPT, GAUSWT, dsf ,
GDSF , ESTIF , EF , GSTIF , gf , vsgd ,
VSGF ,
Z, EZ, KX , KZ, ZlLAY, ZZLAY,
AK, FI_D, D_FI_N, D_FI, SF2
DIMENSION ISGD(NBCMAX,2), ISGF(NBCMAX,2),
X(NNMMAX) , EX(NPEMAX) , EF(EDFMAX) ,
DSF(NPEMAX) , GDSF(NPEMAX) ,
VSGD(NBCMAX) , VSGF(NBCMAX) , GF(NEQMAX) ,

GAUSPT(NGPMAX,NGPMAX) ,

GAUSWT(NGPMAX,NGPMAX),

108

C -__
10
20
C ___
C _--
C -__
30
C ___

4O

ESTIF(EDFMAX,EDFMAX) ,
GSTIF(NEQMAX,BNDMAX) , NODE(NEMMAX,NPEMAX),

Z(NNMMAX) , EZ(NPEMAX), NLAYER(NEMMAX),
KX(nemmax,NLAMAX), KZ(nemmax,NLAMAX),
ZlLAY(nemmax,NLAMAX), ZZLAY(nemmax,NLAMAX),

AK(nemmax,NLAMAX), FI_D(nemmax), D_FI_N(nemmax,NLAMAX),
D_FI(nemmax,NLAMAX), SF2(2,2)

INITIALIZE THE GLOBAL STIFFNESS ARRAY AND FORCE VECTOR.

DO 20 I = 1, NDFMSH
GF(I) = 0.000
DO 10 J = 1, BNDWTH
GSTIF(I,J) = 0.000
CONTINUE
CONTINUE

CALL CR_FI(NLAMAX, NELMSH, NEMMAX,
NLAYER, Kz, ZlLAY, ZZLAY,
AK, FI_D, D_EI_N, D_FI)

BEGIN DO-LOOP ON NUMBER OF ELEMENTS TO CALCULATE THE ELEMENT
MATRICES AND TO ASSEMBLE THE GLOBAL MATRICES.

DO 50 ELMNUM = l, NELMSH
DETERMINE THE LOCAL (ELEMENT) DATA.

DO 30 ENOD = 1, NODELM
GNOD = NODE(ELMNUM,ENOD)
EX(ENOD) = X(GNOD)
EZ(ENOD) = Z(GNOD)

CONTINUE

call STIFF_T(EDFMAX, NPEMAX, NGPMAX, ndfelm,
NGPstf, ELMNUM,
DSF, GDSF, Ex, ESTIF, BF, GAUSPT, GAUSWT,
EZ, KX , KZ, ZlLAY, Z2LAY, NEMMAX,
NLAMAX, NLAYER,AK, FI_D, D_FI, SFZ)

PRINT STIFFNESS MATRIX AND FORCE VECTOR FOR ELEMENT ONE.

IF (PRTSTF .EQ. 1 .AND. ELMNUM .EQ. 1) THEN
WRITE(OUNIT,1000)
DO 40 I = 1, NDFELM
WRITE(OUNIT,1010) I
WRITE(OUNIT,1020) (ESTIF(I,J), J = 1, NDFELM)
CONTINUE
WRITE(OUNIT,1030)
WRITE(OUNIT,1020) (EF(I), I = l, NDFELM)
END IF

109

C ___

END IF

C ___

60

C -...-

C ___

1000
1010
1020
1030
1040
1050

ASSEMBLE THE GLOBAL MATRICES IN BANDED FORM.

CALL ASMBLE (NEQMAX, BNDMAX, EDFMAX, NEMMAX, NPEMAX, NODELM,
NDFNOD, NODE , ESTIF , EF , GSTIF , GF ,
ELMNUM)

END OF LOOP TO COMPUTE GLOBAL STIFFNESS AND FORCE MATRICES.
CONTINUE
IMPOSE SECONDARY BOUNDARY CONDITIONS.

IF (NSGF .GT. 0) THEN
CALL FBNDRY (NEQMAX, NBCMAX, NDFNOD, NSGF , ISGF , VSGF ,
GF )

PRINT GLOBAL STIFFNESS MATRIX AND FORCE VECTOR.

IF (PRTSTF .EQ. 2) THEN
WRITE(OUNIT,1040)
DO 60 I = 1, NDFMSH
WRITE(OUNIT,1010) I
WRITE(OUNIT,1020) (GSTIF(I,J), J = 1, BNDWTH)
CONTINUE
WRITE(OUNIT,1050)
WRITE(OUNIT,1020) (GF(I), I = l, NDFMSH)
END IF

IMPOSE PRIMARY BOUNDARY CONDITIONS.

CALL UBNDRY (NEQMAX, BNDMAX, NBCMAX, NDFNOD, NDFMSH, BNDWTH,
NSGD , ISGD , VSGD , GSTIF , GF )

SOLVE THE SYSTEM OF EQUATIONS.

IRES = 0

CALL SOLVE (NEQMAX, BNDMAX, NDFMSH, BNDWTH, GSTIF , GF ,
IRES )

FORMAT(' ',//3x,'STIEENEss MATRIX FOR ELEMENT ONE: ')

FORMAT(' ',/5X,'ROW: ',IZ,/)

FORMAT(' ',5X,6E12.4)

FORMAT(' ',//3X,'FORCE VECTOR FOR ELEMENT ONE: ',/)

FORMAT(' ',//3X,'GLOBAL STIFFNESS MATRIX: ',/)

FORMAT(' ',//3X,'GLOBAL FORCE VECTOR: ',/)

RETURN

END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C
C
C
C

C

SUBROUTINE: FBNDRY C
PURPOSE: IMPOSE THE PRESCRIBED FORCE B.C.‘S ON THE C
BANDED SYMMETRIC STIFFNESS MATRIX AND MODIFY C

THE FORCE VECTOR. C

C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE FBNDRY (NEQMAX, NBCMAX, NDFNOD, NSGF , ISGF , VSGF ,

110

10

GF )
IMPLICIT LOGICAL (A-Z)

INTEGER NEQMAX, NBCMAX, NDFNOD, NSGF , ISGF , IBC ,
IDOF

DOUBLE PRECISION VSGF , GF
DIMENSION ISGF(NBCMAX,2), VSGF(NBCMAX), GF(NEQMAX)

DO 10 IBC = 1, NSGF

IDOF = NDFNOD * (ISGF(IBC,1) - 1) + ISGF(IBC,2)
GF(IDOF) = GF(IDOF) + VSGF(IBC)

CONTINUE

RETURN

END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C
C
C

C ___
10
20
C ___
C ___
C --..

C
SUBROUTINE: GQUAD C
PURPOSE: ASSIGN VALUES TO THE ARRAYS "GAUSPT" AND "GAUSWT" C

USED IN GAUSSIAN QUADRATURE. C
C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE GQUAD (NGPMAX, GAUSPT, GAUSWT )
IMPLICIT LOGICAL (A-Z)

INTEGER NGPMAX, I , J
DOUBLE PRECISION GAUSPT, GAUSWT

DIMENSION GAUSPT(NGPMAX,NGPMAX), GAUSWT(NGPMAX,NGPMAX)

INITIALIZE ARRAYS.

DO 20 I = l,NGPMAX
DO 10 J = l,NGPMAX
GAUSPT(I,J) = 0.000
GAUSWT(I,J) 0.000
CONTINUE
CONTINUE

GAUSSIAN QUADRATURE OF ORDER 1.

0.000
2.000

GAUSPT(1,1)
GAUSWT(1,1)

GAUSSIAN QUADRATURE OF ORDER 2.

GAUSPT(1,2) = -1.000/SQRT(3.000)
GAUSPT(2,2) = -GAUSPT(1,2)
GAUSWT(1,2) = 1.000

GAUSWT(2,2) = GAUSWT(1,2)

GAUSSIAN QUADRATURE OF ORDER 3.

GAUSPT(1,3) = -SQRT(3.000/5.000)
GAUSPT(2,3) = 0.000
GAUSPT(3,3) = -GAUSPT(1,3)

111

GAUSWT(1,3) = 5.000/9.000
GAUSWT(2,3) 8.000/9.000
GAUSWT(3,3) = GAUSWT(1,3)

C --- GAUSSIAN QUADRATURE OF ORDER 4.

GAUSPT(1,4) = -0.8611363116
GAUSPT(2,4) = "0.3399810436
GAUSPT(3,4) = -GAUSPT(2,4)
GAUSPT(4,4) = -GAUSPT(1,4)
GAUSWT(1,4) = 0.3478548451
GAUSWT(2,4) = 0.6521451549
GAUSWT(3,4) = GAUSWT(2,4)
GAUSWT(4,4) GAUSWT(1,4)
RETURN

END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C

C SUBROUTINE: INPELM

C PURPOSE: READ ELEMENT DATA.
C

C
C
C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE INPELM (OUNIT , nodelm,
ndfnod)
IMPLICIT LOGICAL (A—Z)

integer

NDFNOD = 1
NDFELM = 4
NGPSTF = 2
NGPGRD = 2
NODELM = 4

WRITE(OUNIT,1000)

WRITE(OUNIT,1010) NODELM, ndfnod, ndfelm, ngpstf, ngpgrd
1000 FORMAT(' ',//3X,'ELEMENT DATAz')
1010 FORMAT(' ',

//10X,'NUMBER OF NODES PER ELEMENT, .............. NODELM =',I4,
/10X,'NUMBER OF D.O.F. IN EACH NODE, ............ NDFNOD =',I4,
/10X,'NUMBER OF D.O.F. IN EACH ELEMENT, ......... NDFELM =',I4,

//10X,'INTEGRATION RULE , ........................ NGPSTF =',I4,
/10X,'INTEGRATION POINTS , ...................... NGPGRD =',I4)

RETURN
END

ndfelm,

ngpstf, ngpgrd,

OUNIT , nodelm, ndfelm, ngpstf, ngpgrd, ndfnod

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C

C SUBROUTINE: INPLAY

C PURPOSE: INPUT ELEMENT LAYERS
C

C
C
C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

OUNIT
KZ

SUBROUTINE INPLAY (IUNIT ,
KX ,

, NELMSH, NLAYER,NLAMAX,NEMMAX,
, ZlLAY , ZZLAY)

112

IMPLICIT LOGICAL (A-Z)

INTEGER I,J, IUNIT , OUNIT , NELMSH, NLAYER, NLAMAX,NEMMAX
DOUBLE PRECISION KX , KZ , ZlLAY , ZZLAY
DIMENSION NLAYER(NEMMAX),

KX(nemmax,NLAMAX), KZ(nemmax,NLAMAX),

ZlLAY (nemmax,NLAMAX) , ZZLAY (nemmax, NLAMAX)

DO 80 I=l, NELMSH
READ(IUNIT,*) NLAYER(I)
DO 70 J=1, NLAYER(I)
READ(IUNIT,*) KX(I,J), KZ(I,J), ZlLAY(I,J), ZZLAY(I,J)
70 CONTINUE
80 CONTINUE

WRITE(OUNIT,1000)
DO 100 I=l, NELMSH

WRITE(OUNIT,1010)
WRITE(OUNIT,1020) I, NLAYER(I)
WRITE(OUNIT,1030)
DO 90 J=1, NLAYER(I)
WRITE(OUNIT,1040) J,KX(I,J),KZ(I,J),ZlLAY(I,J),ZZLAY(I,J)

90 CONTINUE
100 CONTINUE
C --- FORMAT STATEMENTS.
1000 FORMAT(' ,///3X,'NUMBERS AND CHARACTERISTICS OF LAYERS: ')

1020 FORMAT(' ,9X,I3,7X,I3)

1030 FORMAT(' ,/6X,'N LAYER',5X,'KX',12X,'KZ',12X,'ZlLAY',9X,
. 'ZZLAY',5X,/)

1040 FORMAT(' ',9X,I3,4X,E10.4,4X,E10.4,4X,E10.4,4X,E10.4)

I

1010 FORMAT(' ',//6X,'ELEMENT',3X,'N LAYER',/)
I
U

RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: INPMSH C
C PURPOSE: INPUT MESH PARAMETERS AND FORM THE MESH. C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE INPMSH (IUNIT , OUNIT , NEMMAX, NPEMAX, NNMMAX, nodelm,
ndfnod,DOM_X, DOM_Z, NELM_X, NELM_Z,
NELMSH, NDFMSH, nodmsh, node, BNDWTH, X, Z)

IMPLICIT LOGICAL (A-Z)

113

integer IUNIT , OUNIT , NEMMAX, NPEMAX,

NNMMAX,

nodelm,

ndfnod, NELM_X, NELM_Z, nodmsh, BNDWTH,

NDFMSH, node, IMESH, EORDER, NXINCR, NZINCR,
I , J , N, NELMSH, NW, NPE
DOUBLE PRECISION X , Z, DOM_X, DOM_Z, DX, DZ
DIMENSION NODE(NEMMAX,NPEMAX), X(NNMMAX), Z(NNMMAX),
DX(NNMMAX) , DZ(NNMMAX)

READ(IUNIT,*) IMESH, DOM_X, DOM_Z, NELM_X, NELM_Z

NODMSH =
NELMSH =

(NELM_X + 1)*(NELM_z + 1)
NELM_X * NELM_Z

IF (IMESH .EQ. 0) THEN
DO 10 I = 1, NELMSH

READ(IUNIT,*) (NODE(I,J), J=1,NODELM)
10 CONTINUE
DO 20 I = l, NODMSH
READ(IUNIT,*) X(I), Z(I)
20 CONTINUE
ELSE
EORDER = 1
NXINCR = NELM_X * EORDER
NZINCR = NELM_Z * EORDER
READ(IUNIT,*) (DX(I), I = l, NXINCR)
READ(IUNIT,*) (DZ(I), I = 1, NZINCR)
NPE = NODELM
CALL MESH (NPEMAX, NEMMAX, NNMMAX, EORDER, NODE , NPE ,
NODMSH, NELMSH, NELM_X, NELM_Z, X , Z ,
. DX , DZ)
END IF
NDFMSH = NDFNOD * NODMSH
BNDWTH = 0
DO 1 N = l,NELMSH
DO 2 I = l,NODELM
DO 3 J = l,NODELM
NW = ( ABS( NODE(N,I) - NODE(N,J) ) + 1 )
IF (BNDWTH .LT. NW) THEN
BNDWTH = NW
END IF
3 CONTINUE
2 CONTINUE
l CONTINUE

WRITE(OUNIT,1000)

WRITE(OUNIT,1010) NELM_X, NELM_Z, NODMSH, NDFMSH, BNDWTH

WRITE(OUNIT,1020)
WRITE(OUNIT,1030)
DO 90 N = l,NELMSH
WRITE(OUNIT,1040) N,
90 CONTINUE

(NODE(N,I), I=1,NODELM)

1000
1010

FORMAT(' ',//3X,'MESH DATA:')
FORMAT(' ',
//10X,'N. OF ELEMENTS IN THE MESH IN X DIRECTION,
/lOX,'N. OF ELEMENTS IN THE MESH IN Z DIRECTION,

114

.NELMSH =',I4,
.NELMSH =',I4,

/1OX,'NUMBER OF NODES IN THE MESH, .............. NODMSH =',I4,

/10X,'NUMBER OF D.O.F. IN THE MESH, ............. NDFMSH =',I4,
/10X,'HALF BAND WIDTH OF GLOBAL STIFFNESS MATRIX, BNDWTH =',I4)
1020 FORMAT(' ',//10X,'BOOLEAN (CONNECTIVITY) MATRIX: ')
1030 FORMAT(' ',/1OX,'ELEMENT NO.',6X,'NODE(1,2,3,4)',/)
1040 FORMAT(' ',11X,I3,10X,9(I3,1X))
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: MESH C
C PURPOSE: GENERATE THE CONNECTIVITY ARRAY NODE(I,J), THE C
C COORDINATES X(I), Z(I), AND MESH INFORMATION C
C (NODMSH, NELMSH, NODELM) FOR RECTANGULAR DOMAINS. C
C THE DOMAIN IS DIVIDED INTO RECTANGULAR ELEMENTS C
C WITH FOUR, EIGHT, OR NINE NODES (NUMELX BY NUMELY C
C NONUNIFORM MESH IN GENERAL). C
C C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE MESH (NPEMAX, NEMMAX, NNMMAX, IEL , NOD , NPE ,
NNM , NEM , NX , NZ , X , Z ,
DX , DZ )

IMPLICIT DOUBLE PRECISION (A-H,O-Z)

INTEGER NPEMAX, NEMMAX, NNMMAX, IEL , NOD , NPE ,
NNM , NEM , NX , NZ

DOUBLE PRECISION X , Z , DX , DZ

DIMENSION X(NNMMAX) , Z(NNMMAX) , NOD(NEMMAX,NPEMAX),
DX(NNMMAX), DZ(NNMMAX)

NEX1=NX+1

NEZl=Nz+1

NXX=IEL*NX

NZZ=IEL*NZ

NXX1=NXX+1

NZZ1=NZZ+1

NEM=NX*NZ
NNM=NXX1*NZZl—(IEL-1)*NX*NZ
IF(NPE.EQ.9)NNM=NXX1*NZZl

C --- GENERATE THE CONNECTIVITY ARRAY NOD(I,J).
KO=O

IF(NPE.EQ.9)K0=1
NOD(1,1) = 1

NOD(1,2) = IEL+1

NOD(1,3) = NXX1+(IEL-1)*NEX1+IEL+1
IF (NPE .EQ. 9)NOD(1,3)=4*NX+5
NOD(1,4) = NOD(1,3) - IEL

IF(NPE .EQ. 4)GO TO 200

NOD(1,5) = 2

NOD(1,6) = NXXl + (NPE-6)

NOD(1,7) = NOD(1,3) - l

NOD(1,8) = NXX1+1
IF (NPE .EQ. 9)NOD(1,9)=NXX1+2
200 IF(NZ .EQ. 1)GOTO 230

115

210
220
230

240

250

260

270

280

290

300

310

320
330

M = 1

DO 220 N = 2,NZ

L = (N-l)*NX + 1

DO 210 I = l,NPE

NOD(L,I) = NOD(M,I)+NXX1+(IEL-1)*NEX1+KO*NX
M=L

IF(NX .EQ .1)GO TO 270

D0 260 NI = 2,NX

DO 240 I = l,NPE

K1 = IEL

IF(I .EQ. 6 .OR. I .EQ. 8)K1=1+K0

NOD(NI,I) NOD(NI-1,I)+K1

M = NI

DO 260 NJ = 2,NZ

L = (NJ-1)*NX+NI

DO 250 J = l,NPE

NOD(L,J) = NOD(M,J)+NXX1+(IEL-l)*NEX1+KO*NX
M = L

GENERATE THE COORDINATES X(I) AND Z(I)

DX(NXX1)=0.0
DZ(NZZl)=0.0

ZC=0.0

IF (NPE .EQ. 9) GOTO 310
D0 300 NI = 1, NEZl

I = (NXX1+(IEL—1)*NEX1)*(NI-1)+1
J = (NI-l)*IEL+1

X(I) = 0.0

Z(I) = ZC

DO 280 NJ = l,NXX

I=I+1

X(I) = X(I-1)+DX(NJ)
Z(I) = ZC

IF(NI.GT.NZ .OR. IEL.EQ.1)GO TO 300
J = J+l

ZC = ZC+DZ(J-l)

I = I+l

X(I) = 0.0

Z(I) = ZC

DO 290 II = 1, NX

K = 2*II-1

I = I+1

X(I) = X(I—1)+DX(K)+DX(K+1)
Z(I) = ZC

ZC = ZC+DZ(J)
RETURN

DO 330 NI=1,NZZl
I=NXX1*(NI—1)
XC=0.0

DO 320 NJ=1,NXX1
I = I+1

X(I) = XC

Z(I) = ZC

XC = XC + DX(NJ)
ZC = ZC + DZ(NI)

RETURN
END

116

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C

SUBROUTINE: INPNBC
PURPOSE: INPUT NODAL BOUNDARY CONDITIONS.

C
C
C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

10

20

30

40

50

SUBROUTINE INPNBC (IUNIT , OUNIT , NBCMAX, NNMMAX, NODMSH, ndfnod,
NSGD , NSGF , ISGD , ISGF , X, Z , vsgd ,
VSGF )

IMPLICIT LOGICAL (A-Z)

INTEGER IUNIT , OUNIT , NBCMAX, NNMMAX, NODMSH, ndfnod,
NSGD , NSGF , ISGD , ISGF ,
I , J , K , IBC , NOD , DOF

DOUBLE PRECISION X , Z, VSGD , VSGF

DIMENSION ISGD(NBCMAX,2), ISGF(NBCMAX,2), IBC(4),
VSGD(NBCMAX) , VSGF(NBCMAX) , X(NNMMAX),
Z(NNMMAX)

READ(IUNIT,*) NSGD, NSGF
IF (NSGD .GT. 0) THEN

READ(IUNIT,*) ((ISGD(I,J),J=1,2),I=1,NSGD)
READ(IUNIT,*) (VSGD(I),I=1,NSGD)

END IF

IF (NSGF .GT. 0) THEN
READ(IUNIT,*) ((ISGF(I,J),J=1,2),I=1,NSGF)
READ(IUNIT,*) (VSGF(I),I=1,NSGF)

END IF

WRITE(OUNIT,1000)
WRITE(OUNIT,1010) NSGD
WRITE(OUNIT,1020)
DO 30 I = l,NODMSH
DO 10 K = l,NDFNOD
IBC(K) = 0
CONTINUE
DO 20 J = l,NSGD
NOD = ISGD(J,1)
DOF = ISGD(J,2)
IF (NOD .EQ. I) IBC(DOF) = DOF
CONTINUE
WRITE(OUNIT,1030) I, X(I), Z(I), (IBC(K), K = 1, NDFNOD)
CONTINUE
WRITE(OUNIT,1035)
WRITE(OUNIT,1070) (VSGD(I), I = 1, NSGD)

IF (NSGF .NE. 0) THEN
WRITE(OUNIT,1040) NSGF
WRITE(OUNIT,1050)

DO 60 I = l,NODMSH
DO 40 K = l,NDFNOD
IBC(K) = O
CONTINUE
DO 50 J = l,NSGF
NOD = ISGF(J,1)
DOF = ISGF(J,2)
IF (NOD .EQ. I) IBC(DOF) = DOF
CONTINUE
WRITE(OUNIT,1030) I, X(I), Z(I), (IBC(K), K = l, NDFNOD)

117

60 CONTINUE
WRITE(OUNIT,1060)
WRITE(OUNIT,1070) (VSGF(I), I = 1, NSGF)

END IF
C --- FORMAT STATEMENTS.
1000 FORMAT(' ',//3X,'NODAL BOUNDARY CONDITIONS: ',/)
1010 FORMAT(' ',//10X,'NUMBER OF SPECIFIED TEMPERATURE D.O.F =',I4)
1020 FORMAT(' ',//lOX,'NODE',4X,'X-COORD.',6X,'Z-COORD.'5X,

'SPECIFIED TEMPERATURE D.O.F.',/)
,9X,I3,4X,E10.4,4X,E10.4,3X,12(I2,1X))
,/10X,'VALUES OF THE SPECIFIED TEMPERATURE D.O.F.:')
,//1OX,'NUMBER OF SPECIFIED GENERALIZED FORCES =',I4)
,//10X,'NODE',4X,'X-COORD.',4X,'Z-COORD.',5X,
. 'SPECIFIED GENERALIZED FORCE',/)
1060 FORMAT(' ',/10X,'VALUES OF THE SPECIFIED GENERALIZED FORCES:')

1030 FORMAT('
1035 FORMAT('
1040 FORMAT('
1050 FORMAT('

1070 FORMAT(' ',/lOX,6(E10.4,2X))
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: INPOPT C
C PURPOSE: READ OPTION PARAMETERS FROM THE INPUT FILE. C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE INPOPT (IUNIT , OUNIT , PRTSTF, PRTGRD)
IMPLICIT LOGICAL (A-Z)
INTEGER IUNIT , OUNIT , PRTSTF, PRTGRD
CHARACTER*72 TITLE

READ(IUNIT,1000) TITLE

READ(IUNIT,*) PRTSTF
C READ(IUNIT,*) PRTSTF, PRTGRD
PRTGRD = 0

WRITE(OUNIT,1010)
WRITE(OUNIT,1020) TITLE

1000 FORMAT(A72)

1010 FORMAT(' ',//3X,
'RESULTS OF "TERM2D VERS. 1.0" - WRITTEN BY ANTONIO PANTANO')
1020 FORMAT(' ',//3X,A72)
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: PRTGRD C
C PURPOSE: WRITE THE SOLUTION GRADIENTS TO THE OUTPUT FILE. C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

subroutine prtdu (ounit , npemax, neqmax, nnmmax, edfmax, nemmax,
ngpmax, nodelm, ndfnod, ndfelm, nelmsh, ngpgrd,
node ,
gauspt, gu , eu , ex , sf ,

118

C -__

10

20

C __-

C ..-...

CA
CA

30

gdsf ,

IMPLICIT LOGICAL (A-Z)

integer ounit ,
ngpmax,
node ,
elmnum,
gdof

double precision gauspt,

gdsf ,

dimension

a0 ,
npemax,
nodelm,
idof ,
9U :
a0 ,
dudx ,

gu(neqmax)

a0(nemmax)
X(nnmmax) ,
gauspt(ngpmax,ngpmax),

WRITE(OUNIT,1000)
WRITE(OUNIT,1010)

I
sf(npemax) ,
I

a1 , a2 ,
neqmax, nnmmax,
ndfnod, ndfelm,
igp , enod ,
eu , ex ,
a1 , a2 ,
adudx , a I
eu(edfmax) I
gdsf(npemax) ,
a1(nemmax) ,

BEGIN LOOP OVER NUMBER OF ELEMENTS IN THE MESH.

DO 50 ELMNUM = l, NELMSH

edfmax, nemmax
nelmsh, ngpgrd

gnod , edof

sf ,
X r
xi , elenth

ex(npemax)

a2(nemmax)

node(nemmax,npemax)

DETERMINE THE LOCAL (ELEMENT) SOLUTION VECTOR AND NODAL

COORDINATES.

gdof = (node(elmnum,1) - 1)

do 10 edof = I, ndfelm
eu(edof) = gu(gdof+edof)

continue

do 20 enod = 1, nodelm
gnod = node(elmnum,enod)

ex(enod) = X(gnod)
continue
elenth = ex(nodelm) — ex(l)

BEGIN LOOP ON NUMBER OF GAUSS POINTS

DO 40 IGP = 1, NGPGRD
XI = GAUSPT(IGP,NGPGRD)

COMPUTE THE SOLUTION GRADIENTS AT THE GAUSS POINTS.

* ndfnod

CALL SHPLlD (NPEMAX, NODELM, XI
GDSF , JACOBN)

DUDX = 0.000
XC = 0.000
DO 30 IDOF = l, NODELM

(NGPGRD).

, ELENTH, SF

DUDX = DUDX + EU(IDOF) * GDSF(IDOF)
XC = XC + EX(IDOF) * SF(IDOF)

CONTINUE

A
ADUDX = A*DUDX

WRITE(OUNIT,1030) ELMNUM, XC, DUDX,

ADUDX

119

, DSF ,

A0(elmnum) + A1(elmnum)*XC + A2(elmnum)*XC*XC

I

I

I

40
50

C _--
1000

1010
1030

CONTINUE
CONTINUE

FORMAT STATEMENTS.

FORMAT(' ',/3X,'SOLUTION GRADIENTS: ',/)

FORMAT(' ',/10X,'ELEMENT NO. x—COORD. DU/Dx A*DU/DX',/)
FORMAT(' ', 11X,I3,6X,4(E10.4,2X))

RETURN

END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C

SUBROUTINE: PRTGU

PURPOSE: WRITE THE SOLUTION VECTOR TO THE OUTPUT FILE.

C
C
C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C ___

20
C _._-
1000

1010
1020

subroutine prtgu (OUNIT , NEQMAX, NNMMAX, NDFNOD, NODMSH,

X, Z )

IMPLICIT LOGICAL (A-Z)

INTEGER OUNIT , NEQMAX, NNMMAX, NDFNOD, NODMSH,

U1 I 02 l I I J

DOUBLE PRECISION GU , X, Z

DIMENSION GU(NEQMAX) , X(NNMMAX), Z(NNMMAX)

WRITE(OUNIT,1000)
OUTPUT OF SOLUTION VECTOR.

WRITE(OUNIT,1010)
U2 = 0
DO 20 I = 1, NODMSH
01 U2 + 1
02 = UZ + ndiOD
WRITE(OUNIT,1020) I, X(I), Z(I), (GU(J), J = U1, U2)
CONTINUE

FORMAT STATEMENTS.

FORMAT(' ',///3X,'SOLUTION: ',/)

FORMAT(' ',1OX,'NODE',7X,'X-COORD.',5X,'Z-COORD. '6X, 'U',/)

FORMAT(' ', 10X,I3,6X,E12.4,2X,E12.4,2X,E12.6)
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C
C
C

C

SUBROUTINE: SHPLlD_T C
PURPOSE: EVALUATE THE l-D LAGRANGIAN INTERPOLATION FUNCTIONS C
AND THEIR GLOBAL DERIVATIVES AT THE GAUSS POINTS. C

C

SUBROUTINE SHPL1D_T (NPEMAX, XI , ELENTH_X, DSF,
GDSF , JACOBN, SF2, IGP)

120

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

IMPLICIT LOGICAL (A-Z)
INTEGER NPEMAX, I, IGP
DOUBLE PRECISION XI , ELENTH_X, DSF , GDSF , JACOBN,
SF2
DIMENSION DSF(NPEMAX), GDSF(NPEMAX), SF2(2,2)
C --- LINEAR LAGRANGE INTERPOLATION FUNCTIONS FOR 2-NODED ELEMENTS.
SF2(1,IGP) = 0.5000 * (1.000 - XI)
SF2(2,IGP) = 0.5000 * (1.000 + XI)
DSF(I) = -0.5000
DSF(2) = 0.5000
C --- COMPUTE THE GLOBAL DERIVATIVES OF SF(I).
JACOBN = ELENTH_X * 0.5000
DO 10 I = 1, 2
GDSF(I) = DSF(I) / JACOBN
10 CONTINUE
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SOLVE: C
C THIS PROGRAM SOLVES A BANDED SYMMETRIC SYSTEM OF EQUATIONS. C
C THE BANDED MATRIX IS INPUT THROUGH BAND(NEQNS,NBW), AND C
C RHS IS THE RIGHT HAND SIDE (FORCE VECTOR) OF THE EQUATION. C
C NEQNS IS THE NO. OF EQUATIONS AND NBW IS THE HALF BAND WIDTH. C
C IN RESOLVING, IRES .GT. 0, LHS ELIMINATION IS SKIPPED. C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C CALL SOLVE (NEQMAX, BNDMAX, NDFMSH, BNDWTH,GSTIF, GF, IRES )
SUBROUTINE SOLVE (NRM , NCM , NEQNS , NBW , BAND, RHS, IRES)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)

INTEGER NRM , NCM , NEQNS , NBW , IRES

DOUBLE PRECISION BAND , RHS

DIMENSION BAND(NRM,NCM), RHS(NRM)

MEQNS=NEQNS-1

IF (IRES.GT.0) GO TO 40

DO 30 NPIV=1,MEQNS

NPIVOT=NPIV+1

LSTSUB=NPIV+NBW-1

IF (LSTSUB.GT.NEQNS) LSTSUB=NEQNS

DO 20 NROW=NPIVOT,LSTSUB
C --- INVERT ROWS AND COLUMNS FOR ROW FACTOR.

NCOL=NROW-NPIV+1

121

10
20
30

40

50
60

C ___

7O

80
90

FACTOR=BAND(NPIV,NCOL)/BAND(NPIV,1)
DO 10 NCOL=NROW,LSTSUB
ICOL=NCOL-NROW+1

JCOL=NCOL~NPIV+1
BAND(NROW,ICOL)=BAND(NROW,ICOL)-FACTOR*BAND(NPIV,JCOL)
RHS(NROW)=RHS(NROW)-FACTOR*RHS(NPIV)
CONTINUE

GO TO 70

DO 60 NPIV=1,MEQNS

NPIVOT=NPIV+1

LSTSUB=NPIV+NBW-1

IF (LSTSUB.GT.NEQNS) LSTSUB=NEQNS
DO 50 NROW=NPIVOT,LSTSUB
NCOL=NROW—NPIV+1
FACTOR=BAND(NPIV,NCOL)/BAND(NPIV,1)
RHS(NROW)=RHS(NROW)-FACTOR*RHS(NPIV)
CONTINUE

BACK SUBSTITUTION.

DO 90 IJK=2,NEQNS
NPIV=NEQNS-IJK+2
RHS(NPIV)=RHS(NPIV)/BAND(NPIV,1)
LSTSUB=NPIV-NBW+1

IF (LSTSUB.LT.1) LSTSUB=1
NPIVOT=NPIV-1

DO 80 JKI=LSTSUB,NPIVOT
NROW=NPIVOT-JKI+LSTSUB
NCOL=NPIV-NROW+1
FACTOR=BAND(NROW,NCOL)
RHS(NROW)=RHS(NROW)-FACTOR*RHS(NPIV)
CONTINUE

RHS(1)=RHS(l)/BAND(1,1)

RETURN

END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C
C
C

SUBROUTINE: STIFF_T

PURPOSE: COMPUTE THE ELEMENTAL STIFFNESS AND FORCE ARRAYS

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE STIFF_T(EDFMAX, NPEMAX, NGPMAX, ndfelm,
NGPstf, ELMNUM,

DSF, GDSF, EX, ESTIF, EF, GAUSPT, GAUSWT,
EZ, KX , KZ, ZlLAY, ZZLAY, NEMMAX,

NLAMAX, NLAYER,AK, FI_D, D_FI,
SF2)

IMPLICIT LOGICAL (A—Z)
integer EDFMAX, NPEMAX, NGPMAX, ndfelm,

NGPstf, I , J . igp:
NLAMAX, NLAYER, NEMMAX, ELMNUM, K

double precision DSF , GDSF , EX , ESTIF , EF
GAUSPT, GAUSWT,
xi , const , jacobn,

ELENTH_X, ELENTH_Z,

122

C
C
C
C
C
C

C ___
4

3

2

10

20

C ___-
21

E2, KX , Kz, ZlLAY, ZZLAY,
AK, FI_D, D_FI, SF2, K_LAY,
FI, D_FI_F, FI_N, H_LAYER, FXI, DXI

dimension dsf(npemax), gdsf(npemax),

EX(NPEMAX), ef(edfmax),

GAUSPT(NGPMAX,NGPMAX), GAUSWT(NGPMAX,NGPMAX) ,

ESTIF(EDFMAX,EDFMAX), EZ(NPEMAX), NLAYER(NEMMAX),

KX(nemmax,NLAMAX), KZ(nemmax,NLAMAX),

ZlLAY(nemmax,NLAMAX), ZZLAY(nemmax,NLAMAX),

AK(nemmax,NLAMAX), FI_D(nemmax),

D_FI(nemmax,NLAMAX), SF2(2,2), K_LAY(4,4,NLAMAX),
FI(2,3), D_FI_F(2), FXI(2,3), DXI(2)

ELENTH_X = EX(2) — EX(I)
ELENTH_Z = EZ(4) - EZ(I)

INITIALIZE THE ARRAYS.

DO 2 I= 1, 4
DO 3 J= 1, 4
DO 4 K= 1, NLAYER(ELMNUM)
K_LAY(I,J,K) = 0.000
CONTINUE
CONTINUE
CONTINUE

DO 20 I = 1, NDFELM
EF(I) = 0.000
DO 10 J = l, NDFELM
ESTIF(I,J) = 0.000
CONTINUE
CONTINUE

BEGIN GAUSS QUADRATURE LOOP FOR FULL INTEGRATION.
DO 21 IGP = 1, NGPstf
XI = GAUSPT(IGP,NGPstf)

CALL SHPLlD_T (NPEMAX, XI, ELENTH_X, DSF,
GDSF , JACOBN, SF2, IGP)

CONST = JACOBN*GAUSWT(IGP,NGPstf)
CONTINUE

DO 22 I= 1, 4
DO 23 J= 1, 4
DO 24 K= 1, NLAYER(ELMNUM)

CALL FI_AT_N (I, J, K, ELMNUM, ZlLAY, ZZLAY, AK,
NEMMAX, NLAMAX, npemax, FI_D, D_FI,
D_FI_F, FI, FI_N, FXI, DXI, SF2, gdsf)
H_LAYER = ( ZZLAY(ELMNUM,K) - ZlLAY(ELMNUM,K) ) / 2

K_LAY(I,J,K)= H_LAYER/3 * (
JACOBN * KX(ELMNUM,K) * ( 2*DXI(1)*DXI(2) ) *

123

24
23

( FI(1,1)*FI(2,1) + 4*FI(1,2)*FI(2,2)

6 * JACOBN * KZ(ELMNUM,K) *
( FXI(1,1)*FXI(2,1) + FXI(1,2)*FXI(2,2) )*
( D_FI_F(1) * D_FI_F(2) ) )

ESTIF(I,J) = ESTIF(I,J) + K_LAY(I,J,K)
CONTINUE
CONTINUE

22 CONTINUE

50 CONTINUE

RETURN
END

+ FI(1,3)*FI(2,3)

)+

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

000000

SUBROUTINE:

PURPOSE:

UBNDRY

IMPOSE THE PRESCRIBED DISPLACEMENT B.C.‘S ON THE
BANDED SYMMETRIC STIFFNESS MATRIX AND MODIFY

THE FORCE VECTOR.

000000

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

20

SUBROUTINE UBNDRY (NEQMAX, BNDMAX, NBCMAX, NDFNOD, NDFMSH, BNDWTH,

NSGD , ISGD , VSGD , GSTIF , GF )

IMPLICIT LOGICAL (A-Z)

INTEGER

NEQMAX, BNDMAX, NBCMAX, NDFNOD, NDFMSH, BNDWTH,

NSGD , ISGD , IBC , IDOF , BWMl , ROW

COL , K

DOUBLE PRECISION VSGD , GSTIF , GF , UVALUE

DIMENSION ISGD(NBCMAX,2), VSGD(NBCMAX), GSTIF(NEQMAX,BNDMAX),
GF(NEQMAX)
DO 40 IBC = 1, NSGD
IDOF = NDFNOD * (ISGD(IBC,1) - 1) + ISGD(IBC,2)
UVALUE = VSGD(IBC)
BWMl = BNDWTH - 1
ROW = IDOF - BNDWTH
DO 20 K - 1, BWMl
ROW = ROW + 1

IF (ROW .GT. 0) THEN

COL = IDOF - ROW + l
GF(ROW) = GF(ROW) - GSTIF(ROW,COL) * UVALUE
GSTIF(ROW,COL) = 0.000

END IF

CONTINUE

GSTIF(IDOF,1) = 1.000
GF(IDOF)

= UVALUE

124

I

ROW = IDOF

DO 30 K = 2, BNDWTH
ROW = ROW + 1

IF (ROW .LE. NDFMSH) THEN

GF(ROW) = GF(
GSTIF(IDOF,K)
END IF
3O CONTINUE
40 CONTINUE
RETURN
END

ROW) - GSTIF(IDOF,K) * UVALUE

= 0.000

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C
C
C

SUBROUTINE:
PURPOSE:

FI AT_N

CHOOSE THE FUNCTION FI AND XI TO BE USED IN THE

EVALUATION OF THE ACTUAL K_LAY

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE FI_AT_N (I, J, K, ELMNUM, ZlLAY, Z2LAY, AK,

IMPLICIT LOGICAL (A-Z)

NEMMAX, NLAMAX, npemax, FI_D, D_FI,
D_FI_F, FI, FI_N, FXI, DXI, SF2, gdsf)

double precision

dimension

integer I, J, K: R: Q,

ZlLAY,

NLAMAX, NEMMAX, ELMNUM, npemax
ZZLAY, AK, FI_D, D_FI,
FI, D FI F, FI_N, zz, FXI, DXI,

gdsflsFE

AK(nemmax,NLAMAX), FI_D(nemmax),
D_FI(nemmax,NLAMAX),

ZlLAY(nemmax,NLAMAX), Z2LAY(nemmax,NLAMAX),

FI(2,3), D_FI_F(2), FXI(2,3), DXI(2),

gdsf(npemax), SF2(2,2)

DO 90 R = 1,2

IF (I .EQ. 1 .OR. I .EQ. 4) THEN

DXI(l) = gdsf(l)
FXI(1,R) = SF2(1,R)
ELSE

DXI(l) = gdsf(2)
FXI(1,R) = SF2 (2,R)
END IF

IF (J .EQ. 1 .OR. J .EQ. 4) THEN
DXI(2) = gdsf(l)
FXI(2,R) = SF2(1,R)
ELSE
DXI(2) = gdsf(2)
FXI(2,R) = SF2(2,R)
END IF

9O CONTINUE

125

C
C
C
C
C
C

DO 110 R = 1,3

FI_N = 0.000
IF (R .EQ. 1) THEN
zz= ZlLAY(ELMNUM,K)-ZlLAY(ELMNUM,1)
ELSE IF (R .EQ. 2) THEN
zz=( (ZZLAY(ELMNUM,K)-ZlLAY(ELMNUM,1)) +
(ZlLAY(ELMNUM,K)-ZILAY(ELMNUM,l)) )/2
ELSE IF (R .EQ. 3) THEN
zz= ZZLAY(ELMNUM,K)-ZILAY(ELMNUM,1)
END IF

FI_N = 22
Do 100 Q = 2, K
FI_N = FI_N +(zz - (ZZLAY(ELMNUM,Q-l)-ZlLAY(ELMNUM,1)))
*AK(ELMNUM,Q—1)

100 CONTINUE

IF (I .LT. 3) THEN

D_FI_F(I) = -D_FI(ELMNUM,K)
FI(1,R)= (l - FI_N/FI_D(ELMNUM) )
ELSE

D_FI_F(l) = D_FI(ELMNUM,K)

FI(1,R) = FI_N/FI_D(ELMNUM)
END IF

IF (J .LT. 3) THEN

D_FI_F(2) = ~DflFI(ELMNUM,K)
FI(2,R) = (1 - FI_N/FI_D(ELMNUM) )
ELSE
D_FI_F(2) = D_FI(ELMNUM,K)
FI(2,R) = FI_N/FI_D(ELMNUM)
END IF

110 CONTINUE

RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C

SUBROUTINE: PR_L_T
PURPOSE: EVALUATE FUNCTIONS FI AND ITS DERIVATIVES

C
C
C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE PR_L_T(OUNIT, NPEMAX,NLAMAX,NELMSH, NEMMAX,
NLAYER, NEQMAX, ZlLAY, Z2LAY,
AK, FI_D, NODE, GUL, GF)

IMPLICIT LOGICAL (A-Z)

integer OUNIT , NPEMAX, NLAMAX, NELMSH,NEMMAX,
NEQMAX, NLAYER, I , J, K1, K2, NODE

double precision ZlLAY, ZZLAY, AK, FI_D,
GUL, GF, T_B,

126

T_T, FI_N, FI_T

DIMENSION NLAYER(NEMMAX), NODE(NEMMAX,NPEMAX),
ZlLAY(nemmax,NLAMAX) , ZZLAY(nemmax,NLAMAX) ,
AK(nemmax,NLAMAX) , FI_D (nemmax) ,

GUL(nemmax,2,NLAMAX), GF(NEQMAX)

DO 10 I = 1, NELMSH
DO 20 J = 1, 2

IF (J .EQ. 1) THEN
T_E = GF(NODE(I,1))
T_T = GF(NODE(I,4))
FI_T = (ZZLAY(I,1)- ZlLAY(I,1)) / FI_D(I)
GUL(I,1,1) = T_B * (1 - FI_T) + T_T * FI_T
FI_T = 0.0000

DO 30 K1 = 2, NLAYER(I)

FI_N = ZZLAY(I,K1) - ZlLAY(I,1)
DO 40 K2 = 2, K1
FI_N = FI_N +(22LAY(I,K1) - ZZLAY(I,K2-1))*AK(I,K2-l)

40 CONTINUE

FI_T = FI_N / FI_D(I)

GUL(I,1,K1) = T_B * (1 - FI_T) + T_T * FI_T
30 CONTINUE

ELSE

T_B = GF(NODE(I,2))
T_T = GF(NODE(I,3))
FI_T = (ZZLAY(I,1)- ZlLAY(I,1))/FI_D(I)
GUL(I,2,1) = T_B * (1 - FI_T) + T_T * FI_T
FI_T = 0.0000

DO 50 K1 = 2, NLAYER(I)

FI_N = Z2LAY(I,K1) - ZlLAY(I,1)
DO 60 K2 = 2, K1
FI_N = FI_N +(Z2LAY(I,K1) - ZZLAY(I,K2-l))*AK(I,K2-1)

60 CONTINUE

FI_T = FI_N / FI_D(I)

GUL(I,2,K1) = T_B * (1 - FI_T) + T_T * FI_T
50 CONTINUE

END IF
20 CONTINUE
10 CONTINUE

WRITE(OUNIT,1000)
DO 100 I=1, NELMSH

WRITE(OUNIT,1010)
WRITE(OUNIT,1020) I, NLAYER(I)
WRITE(OUNIT,1030)
DO 90 J=1, NLAYER(I)
WRITE(OUNIT,1040) J, Z2LAY(I,J), GUL(I,1,J), GUL(I,2,J)

9O CONTINUE
100 CONTINUE

127

C ___

1000
1010
1020
1030

1040

FORMAT STATEMENTS.

FORMAT(' ',///3X,'SOLUTION AT EACH LAYER: ')

FORMAT(' ',//6X,'ELEMENT',3X,'LAYER',/)

FORMAT(' ',9X,I3,7X,I3)

FORMAT(' ',/6X,'LAYER N.',5X,'Z LAYER',6X,'SIDE 1-4',6X,

'SIDE 2-3',/)
FORMAT(' ',9X,I3,4X, E12.4,4X, E12.6,4X, E12.6)

RETURN
END

128

APPENDIX B FORTRAN CODE, SINLOAD

129

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C S Y N L O A D C
C C
C C
C This program solves the two—dimensional the steady—state C
C heat conduction problem in 2-d for laminated plates C
C subjected to a synusoidal temperature distribution C
C on the top and bottom surfaces. C
C C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C INPUT: C
C C
C > TITLE (A70) C
C > NMATER NLAYER C
C C
C DO LOOP: I = 1 TO NMATER C
C > MATNAM(I) C
C > K1(I) K3(I) C
C END OF LOOP C
C C
C DO LOOP: I = 1 TO NLAYER C
C > MATNUM(I) Z(I) C
C END OF LOOP C
C C
C > XLENTH C
C C
C > T_B T_T c
C DO LOOP: I = 1 TO P_SOL C
C > PX(I) C
C END OF LOOP C
C C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

implicit none

INTEGER IUNIT , OUNIT , LAYMAX, matmax, neqmax
PARAMETER (IUNIT = 10, OUNIT = 12, laymax = 20,
matmax = 10, NEQMAX = 2 * laymax)
integer nmater, nlayer, matnum, P_SOL,
ipivot

double precision

character

dimension

Z, K1, K3, KlL, K3L, L, XLENTH, T_B, T_T, PX,
tthick, a , f

matnam*50

K3(matmax),

Z(laymax+1), K1(matmax),
K1L(laymax), K3L(laymax), L(laymax),
matnum(laymax), MATNAM(MATMAX),

f(neqmax), a(neqmax,neqmax), ipivot(neqmax)

PX(5),

open(unit=10,file='Synload.data',status='old')
open(unit=12,file='Synload.out',status='unknown')

130

C --- PREPROCESSOR ---

C ____________________
call inpdat (iunit , ounit , laymax, matmax, nmater, nlayer,
matnum, matnam, P_SOL,
z, K1, K3, XLENTH, T_B, T_T, PX,
tthick)
C _________________
C --- PROCESSOR ---
C _________________
call Solve (ounit , neqmax,1aymax, matmax,
nlayer, matnum, P_SOL, ipivot, Z, K1, K3,
KlL, K3L, L, XLENTH, T_B, T_T, PX, a, f)
C _____________________
C --- POSTPROCESSOR ---
C _____________________
STOP
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: INPDAT C
C PURPOSE: Read input data from a file. C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE INPDAT (iunit , ounit , laymax, matmax, nmater,
nlayer, matnum, matnam, P_SOL,
Z, K1, K3, XLENTH, T_B, T_T, PX,
tthick)
IMPLICIT none

integer iunit , ounit , laymax, matmax, nmater, nlayer,
matnum, P_SOL, i

double precision Z, K1, K3, XLENTH, T_B, T_T, PX,

tthick
CHARACTER TITLE*70 , MATNAM*50
dimension Z(laymax+l), K1(matmax), K3(matmax),

PX(5), matnum(laymax), MATNAM(MATMAX)

C --- READ TITLE AND CONTROL DATA.

READ(iunit,5000) TITLE
READ(iunit,*) nmater, nlayer

C --- READ MATERIAL PROPERTIES.

DO 10 I = 1,nmater
READ(iunit,5010) MATNAM(I)

131

READ(iunit,*) K1(I), K3(I)
10 CONTINUE

C --- READ LAMINATION SCHEME.

tthick = 0.000
DO 20 I = 1,nlayer
READ(iunit,*) MATNUM(I), Z(I)
20 CONTINUE
tthick = Z(nlayer)
READ(iunit,*) xlenth
READ(iunit,*) T_B, T_T
READ(iunit,*) P_SOL
DO 30 I = 1, P_SOL
READ(iunit,*) PX(I)
30 CONTINUE

C --- PRINT INPUT DATA.

WRITE(Ounit,4900)
WRITE(ounit,5020) TITLE
WRITE(ounit,5025)
WRITE(ounit,5030)
DO 40 I = 1, NMATER
WRITE(ounit,5040) I,MATNAM(I), K1(I), K3(I)
4O CONTINUE
WRITE(ounit,5070)
WRITE(ounit,5080)
DO 50 I = l,NLAYER
WRITE(Ounit,5090) I, MATNUM(I), Z(I)
50 CONTINUE
WRITE(Ounit,5100) tthick, xlenth
WRITE(ounit,5200) T_B, T_T
WRITE(ounit,5300)
DO 60 I = 1, P_SOL
WRITE(ounit,5400) PX(I)
60 CONTINUE

C --- FORMAT STATEMENTS.

4900 FORMAT(' ',//3X,'RESULTS FROM PROGRAM "SYNLOAD" - WRITTEN BY ',

. 'ANTONIO PANTANO')

5000 FORMAT(A70)

5010 FORMAT(A50)

5020 FORMAT(' ',//3x,A70)

5025 FORMAT(' ',//3x,'INPUT DATA:')

5030 FORMAT(' ',//6x,'MATERIAL PROPERTIES:')
v

5040 FORMAT(' ,/10X,'MATERIAL NO. ',Il,': ',A50,
//10X,'K1 = ',E12.4,
/10X,'K3 = ',E12.4)

5070 FORMAT(' ',//6X,'LAMINATION SCHEME:')

5080 FORMAT(' ',/10X,'LAYER',3X,'MATERIAL',4X,'Z LAYER',/)
5090 FORMAT(' ',10X,I2,8X,Il,6X,E12.4,6X,E12.4)

I

5100 FORMAT(' ,/10X,'TOTAL THICKNESS OF LAMINATE = ',E12.4,
. /10X,'X-DIMENSION OF LAMINATE = ',E12.4)
5200 FORMAT(' ',/10X,'TEMPERATURE AT THE BOTTOM = ',E12.4,
. /10X,'TEMPERATURE AT THE TOP = ',E12.4)
5300 FORMAT(' ',//6X,'X AT WHICH THE SOLUTION IS REQUIRED:')
5400 FORMAT(' ',/10X,'COORDINATE X = ',E12.4)
RETURN
END

132

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C Subroutine:
C

Solve

C
C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

subroutine Solve
nlayer,

(ounit , neqmax, laymax, matmax,
matnum, P_SOL,

ipivot, Z, K1, K3,

KlL, K3L, L, XLENTH, T_B, T_T, PX, a, f)

implicit none

ounit , laymax, neqmax,
matmax, nlayer,
i! j IkI

integer

double precision X, Z, K1, K3, KlL, K3L,
Tr arprin

dimension

matnum(1aymax),

matnum,
neq, ierror, ii

Z(laymax+1), K1(matmax),
K1L(laymax), K3L(laymax),

ipivot,

P_SOL,

L, XLENTH, T_B,T_T, PX,

K3(matmax),

L(laymaX), PX(5),

a(neqmax,neqmax), f(neqmax), ipivot(neqmax)

DO 5 I = l, NLAYER

K1L(I) = K1(MATNUM(I))

K3L(I) = K3(MATNUM(I))
5 continue

pi = 4.000 * datan(1.0d0)

p = pi / xlenth

DO 7 I = 1, NLAYER

L(I) = P * dsqrt(K1L(I)/K3L(I))

7 continue
neq = 2*n1ayer

write(ounit,lOOO)

C LOOP FOR X = 1 TO P_SOL
DO 25 K = 1, P_SOL
X = PX(K)

C -—- Initialize the arrays.

do 20 i = 1, neq
f(i) = 0.000
do 10 j = 1, neq
a(i,j) = 0.000
10 continue
20 continue

A(1,1) = 1

133

40
30

C ___...

C _._.-

90

1010

1100

1200
1300

25

A(1,2) = 1
A(NEQ,NEQ-l) = dexp( - L(nlayer) * Z(nlayer))
A(NEQ,NEQ) dexp( L(nlayer) * Z(nlayer))

DO 30 I = 1, (nlayer-l)
DO 40 J = 1, 4
IF (J .EQ. 1) THEN
II = I*2
A(II,II—1)= dexp( - L(I) * Z(I))
A(II+1,II-1)= -K3L(I)*L(I)*dexp( - L(I) * Z(I))
ELSE IF (J .EQ. 2) THEN
II = I*2
A(II,II)= dexp( L(I) * Z(I))
A(II+1,II)= K3L(I)*L(I)*dexp( L(I) * Z(I))
ELSE IF (J .EQ. 3) THEN
II = I*2
A(II,II+1)= — dexp( - L(I+l) * Z(I))
A(II+1,II+1)= K3L(I+l)*L(I+1)*dexp( — L(I+l) * Z(I))
ELSE IF (J .EQ. 4) THEN

II = I*2

A(II,II+2)= - dexp( L(I+1) * Z(I))

A(II+1,II+2)= - K3L(I+1)*L(I+l)*dexp( L(I+l) * Z(I))
END IF

continue

continue

F(1) = T_B * DSIN( x * p )

F(NEQ) = T T * DSIN( X * p )

Call a (incore) LINPACK subroutine to solve the equations.

call genfac (a , neqmax, neq , ipivot, ierror)
if (ierror .ne. 0) write(ounit,1010) ierror
call genslv (a , neqmax, neq , ipivot, f )

Write the TEMPERATURE (at x=PX(K)) to the output file.
write(ounit,1100) X
write(ounit,1200)

do 90 i = 1, nlayer-l

IF (I .EQ. 1) THEN

T = f(1)*dexp(- L(l) * Z(I)) + f(2)*dexp(L(l) * Z(I))
ELSE

J = I*2-1

T = f(J)*dexp(- L(I) * Z(I)) + f(J+1)*dexp(L(I) * Z(I))
END IF

write(ounit,1300) i, Z(i), T

continue

format(' ',/3x,'Error in factorization of coefficient matrix.',

/3x,'Leading minor of order ',i4,

' is not positive definite.')
format(' ',//3x,'Temperatures at X=',E12.4)
format(' ',/7x,'Layer', 8x, '2', 15x, 'T',/)
format(' ',/8x, 14, 8x, e10.5, 6x, e11.6)

continue
1000 format(' ',///3x,'SOLUTION')

134

return
end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

THESE SUBROUTINE
USED TO SOLVE TH
THE ROUTINES "DG
"GENSLV", RESPEC
ARE PRIMARILY CO

SUBROUTINE: GEN
PURPOSE: FAC
ELI

ON ENTRY
A DOUBL

LDA INTEG
THE L
N INTEG
THE O
ON RETURN
A AN UP
WHICH
L IS
TRIAN
IPVT INTEG
AN IN

INFO INTEG
= 0
= K

LINPACK. THIS VE
CLEVE MOLER, UNI

MODIFIED COSMETI

00000000000000000000000000000000000000000000000000

G E N S O L V

S ARE MODIFIED VERSIONS OF LINPACK ROUTINES

E DOUBLE PRECISION GENERAL SYSTEM A * X = B.
EFA" AND "DGESL" ARE RENAMED "GENFAC" AND
TIVELY. THE MODIFICATIONS TO THE SUBROUTINES
SMETIC, MAKING THEM UNIFORM WITH OTHER PROGRAMS.

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

FAC

TOR A DOUBLE PRECISION MATRIX BY GAUSSIAN
MINATION.

E PRECISION(LDA,LDA)

THE MATRIX TO BE FACTORED.

ER
EADING DIMENSION OF THE ARRAY A .

ER
RDER OF THE MATRIX A .

PER TRIANGULAR MATRIX AND THE MULTIPLIERS
WERE USED TO OBTAIN IT.

THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE

A PRODUCT OF PERMUTATION AND UNIT LOWER
GULAR MATRICES AND U IS UPPER TRIANGULAR.

ER(LDA)
TEGER VECTOR OF PIVOT INDICES.

ER

NORMAL VALUE.

IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR
CONDITION FOR THIS SUBROUTINE, BUT IT DOES
INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO
IF CALLED. USE RCOND IN DGECO FOR A RELIABLE
INDICATION OF SINGULARITY.

RSION DATED 08/14/78
VERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.

CALLY BY RON AVERILL, 03/14/91.

00000000000000000000000000000000000000000000000000

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE GENFAC (A , LDA , N , IPVT , INFO )
INTEGER LDA , N , IPVT ,INFO ,
IDAMAX, J , K , KPl , L , NMl

DOUBLE PRECISION

135

DIMENSION A(LDA,LDA) , IPVT(LDA)

C --- GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING.
INFO = 0
NM1 = N - 1

IF (NM1 .LT. 1) GO TO 70
DO 60 K = 1, NM1
KP1 = K + 1

C --- FIND L = PIVOT INDEX.

L = IDAMAX(N-K+1,A(K,K),1) + K - 1

IPVT(K) = L

C --- ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED.

IF (A(L,K) .EQ. 0.0D0) GO TO 40

C --- INTERCHANGE IF NECESSARY.

IF (L .EQ. K) GO TO 10
T = A(L,K)
A(L,K) = A(K,K)
A(K,K) T
10 CONTINUE

C --- COMPUTE MULTIPLIERS.
T = -1.0D0/A(K,K)
CALL DSCAL(N-K,T,A(K+1,K),1)
C --- ROW ELIMINATION WITH COLUMN INDEXING.

DO 30 J = KP1, N

T = A(L,J)
IF (L .EQ. K) GO TO 20
A(L,J) = A(K,J)
A(K,J) = T
20 CONTINUE

CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
30 CONTINUE
GO TO 50
40 CONTINUE
INFO = K
50 CONTINUE
60 CONTINUE
70 CONTINUE
IPVT(N) = N
IF (A(N,N) .EQ. 0.0D0) INFO = N

RETURN

END

136

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: GENSLV C
C PURPOSE: SOLVE THE DOUBLE PRECISION SYSTEM A * X B USING C
C THE FACTORS COMPUTED BY GENFAC. C
C C
C ON ENTRY C
C C
C A DOUBLE PRECISION(LDA, N) C
C THE OUTPUT FROM GENFAC. C
C C
C LDA INTEGER C
C THE LEADING DIMENSION OF THE ARRAY A . C
C C
C N INTEGER C
C THE ORDER OF THE MATRIX A . C
C C
C IPVT INTEGER(N) C
C THE PIVOT VECTOR FROM GENFAC. C
C C
C B DOUBLE PRECISION(N) C
C THE RIGHT HAND SIDE VECTOR. C
C C
C JOB INTEGER C
C = 0 TO SOLVE A*X B , C
C = NONZERO TO SOLVE TRANS(A)*X = B WHERE C
C TRANS(A) IS THE TRANSPOSE. C
C (HARDWIRED TO ZERO IN THE PROGRAM) C
C C
C ON RETURN C
C C
C B THE SOLUTION VECTOR X C
C C
C ERROR CONDITION C
C C
C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A C
C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY C
C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER C
C SETTING OF LDA . C
C C
C LINPACK. THIS VERSION DATED 08/14/78 C
C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C
C C
C MODIFIED COSMETICALLY BY RON AVERILL, 03/14/91. C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE GENSLV (A , LDA , N , IPVT , B )

INTEGER LDA , N , IPVT , JOB,

K , KB , L , NM1

DOUBLE PRECISION A , B , DDOT , T

DIMENSION A(LDA,LDA) , B(LDA), IPVT(LDA)

JOB = 0

NM1 = N - 1
IF (JOB .NE. 0) GO TO 50

C --- IF JOB = 0 , SOLVE A*X = B.
C --- FIRST SOLVE L*Y = B.

137

IF (NM1 .LT. 1) GO TO 30
DO 20 K = 1, NM1
L = IPVT(K)
T = B(L)
IF (L .EQ. K) GO TO 10
B(L) = B(K)
B(K) - T
10 CONTINUE
CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1)
20 CONTINUE
30 CONTINUE

C --- NOW SOLVE U*X = Y.

DO 40 KB = l, N
K = N + 1 — KB
B (K) = B (K) /A(K, K)
T = -B(K)
CALL DAXPY(K-1,T,A(1,K),1,B(1),1)
4O CONTINUE
GO TO 100
50 CONTINUE

C --- IF JOB = NONZERO, SOLVE TRANS(A) * X = B.
C --- FIRST SOLVE TRANS(U)*Y = B.

DO 60 K = l, N
T = DDOT(K-1,A(1,K),1,B(1),1)
B(K) = (B(K) - T)/A(K,K)
60 CONTINUE

C --- NOW SOLVE TRANS(L)*X = Y.

IF (NM1 .LT. 1) GO TO 90
DO 80 KB = 1, NM1

K = N - KB
B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1)
L = IPVT(K)
IF (L .EQ. K) GO TO 70
T = B(L)
B(L) = B(K)
B(K) = T
70 CONTINUE

80 CONTINUE
90 CONTINUE
100 CONTINUE
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: DAXPY C
C PURPOSE: MULTIPLY A CONSTANT TIMES A VECTOR PLUS A VECTOR. C
C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C
C AUTHOR: JACK DONGARRA, LINPACK, 3/11/78. C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE DAXPY (N , DA , DX , INCX , DY , INCY )
DOUBLE PRECISION DX(l) , DY(1) , DA

INTEGER I , INCX , INCY , M , MP1 ,

138

 

N

IF (N .LE. 0) RETURN
IF (DA .EQ. 0.0D0) RETURN
IF (INCX .EQ. 1 .AND. INCY .EQ. 1) GO TO 20

C --- CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL TO 1.
IX = l
IY = 1
IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
DO 10 I = 1,N
DY(IY) = DY(IY) + DA*DX(IX)

IX = IX + INCX
IY = IY + INCY
10 CONTINUE

RETURN
C --- CODE FOR BOTH INCREMENTS EQUAL TO 1.
C --- CLEAN-UP LOOP.

20 M = MOD(N,4)
IF (M .EQ. 0) GO TO 40
DO 30 I = 1,M
DY(I) = DY(I) + DA*DX(I)
30 CONTINUE
IF (N .LT. 4) RETURN
40 MP1 = M + 1
DO 50 I = MP1,N,4
DY(I) = DY(I) + DA*DX(I)
DY(I + 1) = DY(I + 1) + DA*DX(I + l)

DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
50 CONTINUE
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C FUNCTION: DDOT C
C PURPOSE: FORM THE DOT PRODUCT OF TWO VECTORS. C
C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C
C AUTHOR: JACK DONGARRA, LINPACK, 3/11/78. C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

DOUBLE PRECISION FUNCTION DDOT (N , DX , INCX , DY ,

INCY )

DOUBLE PRECISION DX(1) , DY(1) , DTEMP

INTEGER I , INCX , INCY , IX , IY , M ,
MP1 , N

DDOT = 0.0D0
DTEMP = 0.000
IF (N .LE. 0) RETURN
IF (INCX .EQ. 1 .AND. INCY .EQ. 1) GO TO 20
C --- CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL TO 1.

IX = 1

139

IY = 1
IF (INCX .LT. 0) IX
IF (INCY .LT. 0) IY
DO 10 I = 1,N
DTEMP = DTEMP + DX(IX)*DY(IY)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
DDOT = DTEMP

(-N+l)*INCX + 1
(-N+l)*INCY + l

RETURN
C --- CODE FOR BOTH INCREMENTS EQUAL TO 1.
C --- CLEAN-UP LOOP.

20 M = MOD(N,5)
IF (M .EQ. 0) GO TO 40
DO 30 I = 1,M
DTEMP DTEMP + DX(I)*DY(I)
30 CONTINUE
IF (N .LT. 5) GO TO 60
40 MP1 = M + 1
DO 50 I = MP1,N,5
DTEMP DTEMP + DX(I)*DY(I) + DX(I + l)*DY(I + l) +
* DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
50 CONTINUE
60 DDOT = DTEMP
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: DSCAL C
C PURPOSE: SCALE A VECTOR BY A CONSTANT. C
C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. C
C AUTHOR: JACK DONGARRA, LINPACK, 3/11/78. C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE DSCAL (N , DA , DX , INCX )
DOUBLE PRECISION DA , DX(1)
INTEGER I , INCX , M , MP1 , N , NINCX

IF (N .LE. 0) RETURN
IF (INCX .EQ. 1) GO TO 20

C --- CODE FOR INCREMENT NOT EQUAL TO 1.

NINCX = N*INCX
DO 10 I = 1,NINCX,INCX
DX(I) = DA*DX(I)
10 CONTINUE

RETURN
C --- CODE FOR INCREMENT EQUAL TO 1.
C --- CLEAN-UP LOOP.

20 M = MOD(N,5)
IF (M .EQ. 0) GO TO 40
DO 30 I = 1,M
DX(I) = DA*DX(I)
30 CONTINUE

140

4O

50

IF (N .LT. 5) RETURN

MP1 = M + 1

DO 50 I = MP1,N,5
DX(I) = DA*DX(I)
DX(I + 1) = DA*DX(I + 1)
DX(I + 2) = DA*DX(I + 2)
DX(I + 3) = DA*DX(I + 3)
DX(I + 4) = DA*DX(I + 4)

CONTINUE

RETURN

END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C
C
C

C ___

C ___.

20

30

FUNCTION: IDAMAX

PURPOSE: FIND THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE.
AUTHOR: JACK DONGARRA, LINPACK, 3/11/78.

INTEGER FUNCTION IDAMAX ( N , DX
DOUBLE PRECISION DX(1) , DMAX
INTEGER I , INCX , IX
IDAMAX = 0

IF (N .LT. l) RETURN

IDAMAX = 1

IF (N .EQ. 1) RETURN
IF (INCX .EQ. 1) GO TO 20

CODE FOR INCREMENT NOT EQUAL TO 1.

IX = 1
DMAX = DABS(DX(1))
IX = IX + INCX
DO 10 I = 2,N
IF (DABS(DX(IX)) .LE. DMAX) GO TO 5
IDAMAX = I
DMAX = DABS(DX(IX))
IX = IX + INCX
CONTINUE
RETURN

CODE FOR INCREMENT EQUAL TO 1.

DMAX = DABS(DX(1))

DO 30 I = 2,N
IF (DABS(DX(I)) .LE. DMAX) GO TO 30
IDAMAX = I
DMAX = DABS(DX(1))

CONTINUE

RETURN

END

141

I

INCX

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

)

C
C
C
C
C
C

APPENDIX C FORTRAN CODE, THERMN L

142

 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C THERM 2 D - NON LINEAR C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

I N P U T

00

TITLE

PRTSTF, PRTLAY, INPLAY, PRTITR
NLSTEP, ITRMAx, EPSILN

IMESH, DOM_X, DOM_Z, NELM_X, NELM_Z

VVVV

IF (IMESH = 0) THEN
DO LOOP: I = 1 TO NELMSH
> NODE(I,J), J = 1 TO NODELM
END OF LOOP
> X(I), Z(I) I = 1 TO NODMSH

ELSE
DO LOOP: I = 1 TO NELM_X
> DX(I)

END OF LOOP
DO LOOP: I = 1 TO NELM_Z

> DZ(I)

END OF LOOP
END IF
> NSGD NSGF
> ISGD(I,J), J = 1 TO 2, I = 1 TO NSGD
> VSGD(I), I = 1 TO NSGD
> ISGF(I,J), J = 1 TO 2, I = 1 TO NSGF
> VSGF(I), I = 1 TO NSGF

IF (INPLAY = 0) THEN
DO LOOP I=1, NELMSH
> NLAYER(I)
DO LOOP J=1, NLAYER(I)
> KX(I,J), KZ(I,J), ZlLAY(I,J), ZZLAY(I,J)
END OF LOOP
END OF LOOP
ELSE
DO LOOP I=1, NELM_Z
> NLATMP
DO 70 K = 1, NELM_X
E = E + 1
NLAYER(E) NLATMP
DO LOOP J = 1, NLATMP
> KXOTMP, KZOTMP, LBlTMP, LBZTMP, ZlLTMP, z2LTMP

KXO(E.J) = KXOTMP
KZO(E,J) = KZOTMP
LAMBDl(E,J) = LBlTMP
LAMBD2(E,J) = LBZTMP
ZlLAY(E,J) = ZlLTMP
ZZLAY(E,J) = z2LTMP

END OF LOOP
END OF LOOP
END OF LOOP
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

V A R I A B L E S

00000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000

PARAMETERS:

143

000000000000000000000000000000000000000000000000000000000000000

IUNIT
OUNIT

BNDMAX

EDFMAX
NBCMAX
NDFMAX
NEMMAX
NEQMAX

NGPMAX
NNMMAX
NPEMAX
NLAMAX

PRTSTF

PRTLAY

PRTITR

NLSTEP
DLAM

ITRMAX
FLAG_D

EPSILN

NDFELM
NODELM
NLAYER

MESH:

BNDWTH

EX

IMESH

INPLAY

NDFMSH
NDFNOD
NELMSH
NODE

INTEGER, UNIT NUMBER OF INPUT FILE.
INTEGER, UNIT NUMBER OF OUTPUT FILE.

INTEGER, COLUMN DIMENSION OF GLOBAL STIFFNESS:
GSTIF IS BANDED IF IMASS = 0,

GSTIF IS SQUARE IF IMASS = 1.

INTEGER, MAXIMUM NUMBER OF D.O.F. PER ELEMENT.
INTEGER, MAXIMUM NUMBER OF BOUNDARY CONDITIONS.
INTEGER, MAXIMUM NUMBER OF D.O.F. PER NODE.
INTEGER, MAXIMUM NUMBER OF ELEMENTS IN THE MESH.
INTEGER, MAXIMUM NUMBER OF EQUATIONS (OR MAXIMUM
NUMBER OF D.O.F. IN THE MESH).

INTEGER, MAXIMUM NUMBER OF GAUSS POINTS.
INTEGER, MAXIMUM NUMBER OF NODES IN THE MESH.
INTEGER, MAXIMUM NUMBER OF NODES PER ELEMENT.
INTEGER, MAXIMUM NUMBER OF LAYERS PER ELEMENT.

CONTROL VARIABLES:

INTEGER, FLAG FOR PRINTING STIFFNESS MATRICES:
PRTSTF = 0 FOR NO PRINT,

PRTSTF = 1 FOR PRINTING STIFFNESS OF ELEMENT ONE,
PRTSTF = 2 FOR PRINTING GLOBAL STIFFNESS.
INTEGER, FLAG FOR PRINTING LAYERS PROPERTY:
PRTSTF = 0 FOR NO PRINT,

PRTSTF = 1 FOR PRINTING.

000000000000000000000000000

INTEGER, FLAG FOR PRINTING SOLUTIONS AT EACH ITERATION

AT THE CENTRAL NODE OF THE MESH:
PRTITR = 0 FOR NO PRINT,
PRTITR = 1 FOR PRINTING.
INTEGER, NUMBER OF LOAD STEP.
DOUBLE PRECISION, LOAD INCREMENT SCALE.
INTEGER, MAXIMUM NUMBER OF ITERATIONS.
INTEGER, FLAG FOR DISPLACEMENT OR FORCE B.C.'S:
FLAG_D = 0 FOR MIXED FORCE AND DISPLACEMENT B.C.'S
FLAG_D = 1 FOR ONLY DISPLACEMENT B.C.'S
DOUBLE PRECISION, TOLLERANCE.

ELEMENT PROPERTIES:

INTEGER, NUMBER OF D.O.F. PER ELEMENT.
INTEGER, NUMBER OF NODES PER ELEMENT.
INTEGER, NUMBER OF LAYERS PER ELEMENT.

INTEGER, HALF BANDWIDTH OF GLOBAL STIFFNESS MATRIX
FOR BENDING PROBLEMS. .
DOUBLE PRECISION(NPEMAX), EX(I) IS THE X-COORDINATE

OF THE I-TH NODE OF THE ELEMENT UNDER CONSIDERATION.

INTEGER, FLAG TO INDICATE HOW MESH WILL BE INPUT:

IMESH = 0 IF THE MESH IS INPUT BY HAND,

IMESH = 1 IF THE MESH IS COMPUTED BY THE PROGRAM.

INTEGER, FLAG TO INDICATE HOW LAYERS WILL BE INPUT:

IMESH = 0 LAYERS ARE INPUT ELEMENT BY ELEMENT,

IMESH = 1 LAYERS ARE INPUT BY GROUP OF HORIZZONTAL
ELEMENTS.

INTEGER, NUMBER OF D.O.F. IN THE MESH.

INTEGER, NUMBER OF D.O.F. PER NODE.

INTEGER, NUMBER OF ELEMENTS IN THE MESH.

INTEGER(NEMMAX,NPEMAX), CONNECTIVITY MATRIX:

NODE(I,J) IS THE J-TH NODE OF THE I-TH ELEMENT.

144

00000000000000000000000000000000000

NODMSH
X

NELM_X
NELM_Z
DOM_X
DOM_Z
DX

DZ

KX

KZ
ZlLAY

Z2LAY

INTEGER, NUMBER OF NODES IN THE MESH.

DOUBLE PRECISION(NNMMAX), X(I) IS THE GLOBAL
COORDINATE OF THE I-TH NODE.

INTEGER, NUMBER OF ELEMENTS IN THE MESH IN X DIR.
INTEGER, NUMBER OF ELEMENTS IN THE MESH IN Z DIR.
DOUBLE PRECISION, DOMAIN IN X DIRECTION

DOUBLE PRECISION, DOMAIN IN Z DIRECTION

DOUBLE PRECISION, DIMENSION OF ELEMENTS

IN X DIRECTION STARTING FROM THE LEFT

DOUBLE PRECISION, DIMENSION OF ELEMENTS

IN Z DIRECTION STARTING FROM THE BOTTOM

DOUBLE PRECISION, TERMAL CONDUCTIVITY COEFFICIENT OF

THE LAYER ALONG X DIRECTION

DOUBLE PRECISION, TERMAL CONDUCTIVITY COEFFICIENT OF

THE LAYER ALONG Z DIRECTION

DOUBLE PRECISION, Z COORDINATE OF THE BOTTOM OF THE

LAYER
DOUBLE PRECISION, Z COORDINATE OF THE TOP OF THE
LAYER

BOUNDARY CONDITIONS:

ISGD
ISGF
NSGD

NSGF

VSGF

INTEGER(NBCMAX,2), ISGD(I,J) INDICATES THAT THE J-TH

DISPLACEMENT D.O.F. OF THE I-TH NODE IS SPECIFIED.

INTEGER(NBCMAX,2), ISGF(I,J) INDICATES THAT THE J-TH

FORCE COMPONENT OF THE I-TH NODE IS SPECIFIED.

INTEGER, NUMBER OF DISPLACEMENT BOUNDARY CONDITIONS.
INTEGER, NUMBER OF NODAL FORCE BOUNDARY CONDITIONS.

DOUBLE PRECISION(NBCMAX), VSGD(I) IS THE SPECIFIED
VALUE OF THE I-TH DISPLACEMENT BOUNDARY CONDITION.
DOUBLE PRECISION(NBCMAX), VSGF(I) IS THE SPECIFIED
VALUE OF THE I-TH NODAL FORCE BOUNDARY CONDITION.

FINITE ELEMENT MODEL:

DSF

EF
ESTIF
EU
GAUSPT
GAUSWT
GDSF
GF
GSTIF
NGPSTF
NGPGRD

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C VSGD
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C SF
C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

DOUBLE PRECISION(NPEMAX), DSF(I) IS THE LOCAL (XI)
DERIVATIVE OF THE I-TH SHAPE FUNCTION.
DOUBLE PRECISION(EDFMAX), ELEMENT FORCE VECTOR.

DOUBLE PRECISION(EDFMAX,EDFMAX), ELEMENT STIFFNESS.
DOUBLE PRECISION(EDFMAX), ELEMENT EXTENSION VECTOR.

DOUBLE PRECISION(NGPMAX,NGPMAX) GAUSPT(I,J) IS THE
I-TH GAUSS POINT OF J-TH GAUSS RULE.

DOUBLE PRECISION(NGPMAX,NGPMAX) GAUSWT(I,J) IS THE
I-TH GAUSS WEIGHT OF J-TH GAUSS RULE.

DOUBLE PRECISION(NPEMAX), GDSF(I) IS THE GLOBAL (X)
DERIVATIVE OF THE I-TH SHAPE FUNCTION.

DOUBLE PRECISION(NEQMAX), GLOBAL FORCE VECTOR,
WHICH, ON RETURN, CONTAINS THE SOLUTION VECTOR;
ALSO USED AS A DUMMY VECTOR IN EIGENPROBLEMS.
DOUBLE PRECISION(NEQMAX,BNDMAX), GLOBAL STIFFNESS.
INTEGER, NUMBER OF GAUSS POINTS TO BE USED TO
INTEGRATE THE STIFFNESS MATRIX.

INTEGER, NUMBER OF GAUSS POINTS AT WHICH THE
GRADIENTS ARE TO BE EVALUATED.

INTEGER(NPEMAX), SF(I) IS THE I-TH SHAPE FUNCTION.

IMPLICIT LOGICAL (A-Z)

INTEGER

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

IUNIT , OUNIT , NEMMAX, NPEMAX, NDFMAX, NGPMAX,

NNMMAX, NBCMAX, NEQMAX, EDFMAX, BNDMAX,
NELM_X, NELM_Z, NLAMAX,

145

NLSTEP, ITRMAX, FLAG_D
PARAMETER (IUNIT = 10, OUNIT = 12, NEMMAX = 400,
NPEMAX = 4, NDFMAX = 1, NGPMAX = 5,
BNDMAX = 400, NLAMAX = 40,
EDFMAX = NDFMAX * NPEMAX,
NNMMAX = NEMMAX * (NPEMAX - 1) + 1,
NEQMAX = NNMMAX * NDFMAX,
NBCMAX = NEMMAX*4)
INTEGER PRTSTF, PRTGRD, nodelm, ndfelm, ngpstf, ngpgrd,
nelmsh, ndfmsh, nodmsh, ndfnod, node , bndwth,
nsgd , nsgf , isgd , isgf, PRTLAY, INPLAY,
NLAYER, FLAGNR, PRTITR
DOUBLE PRECISION X , vsgd, vsgf , ex , DOM_X, DOM_Z,
GAUSPT, GAUSWT, ESTIF , EF , GSTIF ,
GF , DSF , GDSF ,
Z , KX , KZ, ZlLAY, ZZLAY, EZ,
AK , FI_D , D_FI_N, D_FI, SF2, GUL,
LAMBDI, LAMBDZ, DLAM, EPSILN, GU, GUPLS,
KXO , KZO, T_LAY
DIMENSION NODE(NEMMAX,NPEMAX), ISGD(NBCMAX,2),
ISGF(NBCMAX,2), X(NNMMAX) , EX(NPEMAX) ,
VSGD(NBCMAX) , VSGF(NBCMAX) ,
EF(EDFMAX) , GF(NEQMAX) ,
DSF(NPEMAX) , GDSF(NPEMAX) ,
GAUSPT(NGPMAX,NGPMAX), GAUSWT(NGPMAX,NGPMAX),
ESTIF(EDFMAX,EDFMAX) ,
GSTIF(NEQMAX,BNDMAX),
Z(NNMMAX), NLAYER(NEMMAX),
KX(nemmax,NLAMAX), KZ(nemmax,NLAMAX),
ZlLAY(nemmax,NLAMAX), ZZLAY(nemmax,NLAMAX),
EZ(NPEMAX), AK(nemmax,NLAMAX), FI_D(nemmax),
D_FI_N(nemmax,NLAMAX), D_FI(nemmax,NLAMAX),
SF2(2,2),GUL(nemmax,2,NLAMAX),
LAMBDl(nemmax,NLAMAX), LAMBD2(nemmax,NLAMAX),
GU(NEQMAX), GUPLS(NEQMAX), T_LAY(nemmax,NLAMAX),
KXO(nemmax,NLAMAX), KZO(nemmax,NLAMAX)
Open(unit=10,file='TERMNL.data',status='old')
open(unit=12,fi1e='TERMNL.out', status='unknown')
C ____________________
PREPROCESSOR ---
C ____________________
INPUT OPTIONS.
CALL INPOPT (IUNIT, OUNIT, PRTSTF, PRTGRD, PRTLAY, INPLAY, PRTITR)

C ___

INPUT ELEMENT DATA.

146

call INPELM (OUNIT , NODELM, NDFELM, NGPSTF, NGPGRD,
NDFNOD)

INPUT NON LINEAR ANALYSIS DATA

call INP_NL (IUNIT , OUNIT , NLSTEP, DLAM, ITRMAX, EPSILN)

INPUT MESH PARAMETERS AND FORM THE MESH.

call INPMSH (IUNIT , OUNIT , NEMMAX, NPEMAX, NNMMAX, NODELM,
NDFNOD, DOM_X, DOM_Z, NELM_X, NELM_Z,
NELMSH, NDFMSH, NODMSH, NODE, BNDWTH, x, z )

INPUT NODAL BOUNDARY CONDITIONS.

call INPNBC (IUNIT , OUNIT , NBCMAX, NNMMAX, NODMSH, NDFNOD,
NSGD , NSGF , ISGD , ISGF , x , z, VSGD,
VSGF, FLAG_D)

INPUT ELEMENT LAYERS

call INLAYR (IUNIT , OUNIT , NELMSH, NLAYER, NLAMAX, NEMMAX,

Kxo, Kzo, ZlLAY, z2LAY, LAMBDl, LAMED2, PRTLAY,

NELM_X, NELM_Z, INPLAY)

INITIALIZE THE ARRAYS WHICH CONTAIN THE GAUSS POINTS AND WEIGHTS.

CALL GQUAD (NGPMAX, GAUSPT, GAUSWT)

IF (FLAG_D .EQ. 0) THEN

CALL PROCSR_F (OUNIT , NEMMAX, NPEMAX, NGPMAX, NNMMAX,
NEQMAX, BNDMAX, NBCMAX, EDFMAX,

ndfnod, ndfelm, nodelm,
ndfmsh, nelmsh, node , ngpstf, prtstf, bndwth,
nsgd , isgd , nsgf , isgf , X , EX,
GAUSPT, GAUSWT, dsf , GDSF , ESTIF,
EF , GSTIF , gf , vsgd , VSGF ,
Z , EZ, KX , KZ, ZlLAY, ZZLAY,
NLAMAX, NLAYER,AK, FI_D, D_FI_N, D_FI, SF2,
NLSTEP, DLAM, ITRMAX, EPSILN,
GU, GUPLS, FLAGNR, KXO, KZO, T_LAY,
NELM_X, NELM_Z, LAMBDl, LAMBDZ, PRTITR)

ELSE

CALL PROCSR_D (OUNIT , NEMMAX, NPEMAX, NGPMAX, NNMMAX,
NEQMAX, BNDMAX, NBCMAX, EDFMAX,

ndfnod, ndfelm, nodelm,
ndfmsh, nelmsh, node , ngpstf, prtstf, bndwth,
nsgd , isgd , X , EX,
GAUSPT, GAUSWT, dsf , GDSF , ESTIF,
EF , GSTIF , gf , vsgd ,
Z , EZ, KX , KZ, ZlLAY, ZZLAY,
NLAMAX, NLAYER,AK, FI_D, D_FI_N, D_FI, SF2,
NLSTEP, DLAM, ITRMAX, EPSILN,
GU, GUPLS, FLAGNR, KXO, KZO, T_LAY,

147

NELM_X, NELM_Z, LAMBDl, LAMBD2, PRTITR)

END IF
C _____________________
C -—— POSTPROCESSOR —--
C _____________________
CAll PRTGU (OUNIT , NEQMAX, NNMMAX, NDFNOD, NODMSH, GU ,
x, z )
CAll PR_L_T (OUNIT, NPEMAX, NLAMAX, NELMSH,
NEMMAX, NLAYER, NEQMAX,
ZlLAY, Z2LAY, AK, FI_D, NODE, GUL, GU, T_LAY)
STOP
END
C ---------------------------------------------------------------------- C
C C
C SUBROUTINES C
C C
C ---------------------------------------------------------------------- c

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: Procsr_F C
C PURPOSE: DRIVER ROUTINE TO SOLVE THE EQUATIONS. C
C FORCE STEPPING C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE PROCSR_F (OUNIT , NEMMAX, NPEMAX, NGPMAX, NNMMAX,

IMPLICIT LOGICAL

integer

NEQMAX,

ndfnod,
ndfmsh,
nsgd ,
GAUSPT,
EF ,
Z ,
NLAMAX,
NLSTEP,

ndfelm,
nelmsh,
isgd ,
GAUSWT,
GSTIF ,
EZ, KX

BNDMAX, NBCMAX,

nodelm,
node
nsgf
dsf
gf

, KZ,

\ ‘ ‘ ‘

EDFMAX,

ngpstf, prtstf, bndwth,

isgf ,
GDSF ,
vsgd ,

ZlLAY,

X :
ESTIF,
VSGF ,

Z2LAY,

NLAYER,AK, FI_D, D_FI_N, D_FI,

DLAM, ITRMAX,

EPSILN,

GU, GUPLS, FLAGNR, KXO, KZO, T_LAY,
NELM_X, NELM_Z, LAMBDl, LAMBD2, PRTITR)

(A-Z)

OUNIT ,
BNDMAX,
ndfmsh,
nsgd ,
ELMNUM,
NLSTEP,
FLAGNR,

NEMMAX,
NBCMAX,
nelmsh,
isgd ,
ENOD ,
ITRMAX,

NPEMAX,
EDFMAX,
node ,
nsgf ,
GNOD ,
ILSTEP,

NGPMAX,
ndfnod,
ngpstf,
isgf ,
IRES ,
ICNVRG,

NNMMAX,
ndfelm,
prtstf,

1:
NLAMAX,
PRTITR,

EX,

SF2,

NEQMAX,
nodelm,
bndwth,

J,
NLAYER,

ITER, F_K0, NELM_X, NELM_Z, CEN, Flag_P

148

double precision X , EX , GAUSPT, GAUSWT, dsf ,
GDSF , ESTIF , EF , GSTIF , gf , vsgd ,
VSGF , Z, EZ, KX , KZ, ZlLAY, ZZLAY,
AK, FI_D, D_FI_N, D_FI, SF2,
LAMBDA, LAMBDl, LAMBD2, DLAM, EPSILN, GU,GUPLS,
ERROR, KXO, KZO, T_LAY, IVSGD, GREF,
GF2, GUZ, UINC, UOLD, DELTA_D

DIMENSION ISGD(NBCMAX,2), ISGF(NBCMAX,2),

X(NNMMAX) , EX(NPEMAX) , EF(EDFMAX) ,
DSF(NPEMAX) , GDSF(NPEMAX) ,

VSGD(NBCMAX) , VSGF(NBCMAX) , GF(NEQMAX) ,
GAUSPT(NGPMAX,NGPMAX) , GAUSWT(NGPMAX,NGPMAX),
ESTIF(EDFMAX,EDFMAX) ,

GSTIF(NEQMAX,BNDMAX) , NODE(NEMMAX,NPEMAX),
Z(NNMMAX) , EZ(NPEMAX), NLAYER(NEMMAX),

KX(nemmax,NLAMAX), KZ(nemmax,NLAMAX),

ZlLAY(nemmax,NLAMAX), ZZLAY(nemmax,NLAMAX),

AK(nemmax,NLAMAX), FI_D(nemmax), D_FI_N(nemmax,NLAMAX),

D_FI(nemmax,NLAMAX), SF2(2,2),

LAMBD1(nemmax,NLAMAX), LAMBD2(nemmax,NLAMAX),

GU(NEQMAX), GUPLS(NEQMAX), T_LAY(nemmax,NLAMAX),
KXO(nemmax,NLAMAX), KZO(nemmax,NLAMAX),

IVSGD(NBCMAX), GREF(NEQMAX),

DELTA_D(NEQMAX)

C ----INITIALIZE LOAD SCALE
LAMBDA = 0
C ----INITIALIZE SOLUTION VECTOR

DO 5 I = 1, NEQMAX
GU(I) = 0.000
GUPLS(I) = 0.000

5 CONTINUE

DO 43 I = 1, NBCMAX
IVSGD(I) = VSGD(I)
43 CONTINUE

C ----- COMPUTE REFERENCE LOAD

IF (NSGF .GT. 0) THEN
CALL FBNDRY (NEQMAX, NBCMAX, NDFNOD, NSGF, ISGF, VSGF, GF)

END IF

DO 57 I = l, NEQMAX
GREF(I) = GF(I)
57 CONTINUE

DO 8 I = 1, NDFMSH

GF(I) 0.000
8 CONTINUE

149

C ----LOOP ON N. OF LOAD STEP
DO 500 ILSTEP = 1, NLSTEP
C ---- INITIALIZE COUNTERS AND FLAG
ICNVRG =
ITER =

ERROR =
FLAGNR =

0
0
1.0
0
C ---- ITERATION LOOP

1 ITER = ITER + 1

C ---- CHECK FOR ITERMAX

IF (ITER .GT. ITRMAX) STOP

NDFMSH
0.000

DO 88 I = 1,
GF(I) =
88 CONTINUE
IF (FLAGNR .EQ. 0) THEN
DO 21 I = 1, NDFMSH
DO 11 J = l, BNDWTH
GSTIF(I,J) = 0.000
CONTINUE
CONTINUE
END IF

11
21

C ___.

F_KO = 0

IF (ILSTEP .EQ. 1
F_KO = 1

END IF

.AND. ITER .EQ.

NLAMAX,
NEQMAX,

NELMSH,

CALL CALC_T (NPEMAX,

NLAYER,

1)

INITIALIZE THE GLOBAL STIFFNESS ARRAY AND FORCE VECTOR.

VERIFY IF WE ARE AT THE FIRST ITERATION OF THE FIRST LOAD STEP.

THEN

NEMMAX,

ZlLAY, ZZLAY, AK,

FI_D, NODE, GU, T_LAY, KXO, KZO, F_K0,

LAMBDl, LAMED2)

CALL CR_FI(NLAMAx, NELMSH, NEMMAX,
NLAYER, Kz, ZlLAY, z2LAY,
AK, FI_D, D_FI_N, D_FI)

C __-
C ___

C -.__.

BEGIN DO-LOOP ON NUMBER OF ELEMENTS TO CALCULATE THE ELEMENT
MATRICES AND TO ASSEMBLE THE GLOBAL MATRICES.
NELMSH

DO 50 ELMNUM = l,

DETERMINE THE LOCAL (ELEMENT) DATA.

DO 30 ENOD = 1, NODELM

150

GNOD = NODE(ELMNUM,ENOD)
EX(ENOD) = X(GNOD)
EZ(ENOD) = Z(GNOD)

3O CONTINUE

CALL STIFF_T (EDFMAX, NPEMAX, NGPMAX, ndfelm,
NGPstf, ELMNUM,
DSF, GDSF, EX, ESTIF, EF, GAUSPT, GAUSWT,
EZ, KX, KZ, ZlLAY, ZZLAY, NEMMAX,
NLAMAX, NLAYER, AK, FI_D, D_FI, SF2)

C --- PRINT STIFFNESS MATRIX AND FORCE VECTOR FOR ELEMENT ONE.

IF (PRTSTF .EQ. 1 .AND. ELMNUM .EQ. 1 .AND. F_KO .EQ. 1) THEN
WRITE(OUNIT,1000)
DO 40 I = 1, NDFELM
WRITE(OUNIT,1010) I
WRITE(OUNIT,1020) (ESTIF(I,J), J = 1, NDFELM)
4O CONTINUE
WRITE(OUNIT,1030)
WRITE(OUNIT,1020) (EF(I), I = l, NDFELM)

END IF
C --- ASSEMBLE THE GLOBAL MATRICES IN BANDED FORM.
CALL ASMBLE (NEQMAX, BNDMAX, EDFMAX, NEMMAX, NPEMAX, NODELM,
NDFNOD, NODE , ESTIF , EF , GSTIF , GF ,
ELMNUM)
C --- END OF LOOP TO COMPUTE GLOBAL STIFFNESS AND FORCE MATRICES.
50 CONTINUE
C ----COMPUTE RESIDUAL LOAD VECTOR

DO 14 I = 1, NEQMAX
GF(I) = (LAMBDA + DLAM) * GREF(I)
14 CONTINUE

C --- PRINT GLOBAL STIFFNESS MATRIX AND FORCE VECTOR. --
C --- IMPOSE PRIMARY BOUNDARY CONDITIONS.
CALL UBNDRY (NEQMAX, BNDMAX, NBCMAX, NDFNOD, NDFMSH, BNDWTH,

NSGD , ISGD , VSGD , GSTIF , GF )

C --* SOLVE THE SYSTEM OF EQUATIONS.

IRES = 0
CALL SOLVE (NEQMAX, BNDMAX, NDFMSH, BNDWTH, GSTIF, GF, IRES)

C --- COMPUTE THE INCREMENT.
DO 113 I = 1, NEQMAX
DELTA_D(I) = GF(I) - GU(I)
113 CONTINUE

151

 

C ----- CHECK FOR CONVERGENGE

GF2 = 0.0000
GU2 = 0.0000
DO 63 I = 1, NEQMAX
GF2 = GF2 + DELTA_D(I)*DELTA_D(I)
GU2 = GUZ + GU(I)*GU(I)
63 CONTINUE

IF (ITER .GT. 1) THEN
UINC = GF2
UOLD = GU2
ERROR = DSQRT(UINC/UOLD)
IF (ERROR .LT. EPSILN) THEN

ICNVRG = 1
END IF
END IF
C ----- UPDATE SOLUTION

DO 73 I = 1, NEQMAX
GU(I) = GF(I)
73 CONTINUE

IF (ICNVRG .EQ. 1) THEN
DO 83 I = 1, NEQMAX
GUPLS(I) = GU(I)
83 CONTINUE
END IF

C ----- PRINT VALUE OF T AT THE CENTER OF THE MESH OR NEAR

IF (Flag_P .EQ. 0) THEN
WRITE(OUNIT,2000)
WRITE(OUNIT,2030)

end if
Flag_P = 1
IF (PRTITR .EQ. 1 .and. ICNVRG .EQ. I ) THEN
CEN = (NELM_X+1)*(NELM_Z/2) + (NELM_X/2+1)
WRITE(OUNIT,2040) GU(CEN)
2000 FORMAT(' ',/3X,'SOLUTION AT THE CENTER OF THE MESH: ')
2030 FORMAT(' ',8X,'T (CENTER)')
2040 FORMAT(' ',SX, E14.8,6X)
END IF
C ----INCREMENT LOAD

IF (ICNVRG .EQ. 1) THEN

C---- CONVERGENCE OBTAINED, NEXT LOAD STEP

LAMBDA = LAMBDA + DLAM
GOTO 500

ELSE

152

C---- CONVERGENCE NOT OBTAINED, NEXT ITERATION
GOTO 1
END IF
500 CONTINUE
1000 FORMAT(' ',//3X,'STIFFNESS MATRIX FOR ELEMENT ONE:
1010 FORMAT(' ',/5X,'ROW: ',I2,/)
1020 FORMAT(' ',5X,6E12.4)
1030 FORMAT(' ',//3X,'FORCE VECTOR FOR ELEMENT ONE:
1040 FORMAT(' ',//3X,'GLOBAL STIFFNESS MATRIX:
1050 FORMAT(' ',//3X,'GLOBAL FORCE VECTOR:
RETURN
END

"/)

':/)

',/)

')

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: Procsr_D C
C PURPOSE: DRIVER ROUTINE TO SOLVE THE EQUATIONS. C
C DISPLACEMENT STEPPING C
C C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE PROCSR_D (OUNIT , NEMMAX, NPEMAX, NGPMAX, NNMMAX,

NEQMAX,

ndfnod,
ndfmsh,
nsgd ,
GAUSPT,
EF ,
Z I
NLAMAX,
NLSTEP,
GU,
NELM_X,

IMPLICIT LOGICAL (A-Z)

ndfelm,
nelmsh,
isgd ,
GAUSWT,
GSTIF ,
EZ, KX

BNDMAX, NBCMAX,

nodelm,
node ,
X ,
dsf ,
9f I

, KZ,

ngpstf,
EX,
GDSF ,
vsgd ,
ZlLAY,

EDFMAX,

prtstf,
ESTIF,

ZZLAY,

bndwth,

GUPLS,

NLAYER,AK, FI_D, D_FI_N, D_FI, SF2,
DLAM, ITRMAX, EPSILN,

FLAGNR, Kxo, KZO, T_LAY,

NELM_z, LAMBDl, LAMBD2, PRTITR)

integer

double precision

GU,GUPLS,

OUNIT ,
BNDMAX,
ndfmsh,
nsgd ,
ELMNUM,
NLSTEP,
FLAGNR,

X
GDSF ,

z, Ez,
FI_D,

AK,

NEMMAX,
NBCMAX,
nelmsh,
isgd ,
ENOD ,
ITRMAX,

NPEMAX,
EDFMAX,
node ,

I,
GNOD ,
ILSTEP,

NGPMAX,
ndfnod,
ngpstf,

J.
IRES ,
ICNVRG,

NNMMAX,
ndfelm,
prtstf,

NLAMAX,
PRTITR,

NEQMAX,
nodelm,
bndwth,

NLAYER,

ITER, F_K0, NELM_X, NELM_Z, CEN, F1ag_P

, EX

ESTIF ,
KX

I

EF ,
r KZ,

153

GAUSPT,

GSTIF ,
ZlLAY,

GAUSWT,

LAMBD2,

dsf

9f I
ZZLAY,

vsgd ,

D_FI_N, D_FI, SF2, DELTA_D, ERROR,

LAMBDA, LAMBDl, DLAM, EPSILN,

KXO, KZO, T_LAY, IVSGD,
UINC, UOLD, GF2, GU2

DIMENSION ISGD(NBCMAX,2),

X(NNMMAX) , EX(NPEMAX) , EF(EDFMAX) ,
DSF(NPEMAX) , GDSF(NPEMAX) ,
VSGD(NBCMAX) , GF(NEQMAX) ,

GAUSPT(NGPMAX,NGPMAX) , GAUSWT(NGPMAX,NGPMAX),
ESTIF(EDFMAX,EDFMAX) ,
GSTIF(NEQMAX,BNDMAX) , NODE(NEMMAX,NPEMAX),
Z(NNMMAX) , EZ(NPEMAX), NLAYER(NEMMAX),
KX(nemmax,NLAMAX), KZ(nemmax,NLAMAX),
ZlLAY(nemmax,NLAMAX), ZZLAY(nemmax,NLAMAX),
AK(nemmax,NLAMAX), FI_D(nemmax), D_FI_N(nemmax,NLAMAX),
D_FI(nemmax,NLAMAX), SF2(2,2),
LAMBDl(nemmax,NLAMAX), LAMBD2<nemmax,NLAMAX),
GU(NEQMAX), GUPLS(NEQMAX), T_LAY(nemmax,NLAMAX),
KXO(nemmax,NLAMAX), KZO(nemmax,NLAMAX),
IVSGD(NBCMAX), DELTA_D(NEQMAX)

C ----INITIALIZE LOAD SCALE
LAMBDA = 0
C ~---INITIALIZE SOLUTION VECTOR

DO 5 I = 1, NEQMAX
GU(I) = 0.000
GUPLS(I) = 0.000

5 CONTINUE

DO 88 I =
GF(I)
88 CONTINUE

1, NEQMAX
= 0.000

C ----EBC'S COPIED IN VECTOR IVSGD
DO 43 I = 1, NBCMAX
IVSGD(I) = VSGD(I)
43 CONTINUE

C ----LOOP ON N. OF LOAD STEP

DO 500 ILSTEP = 1, NLSTEP

C ---- INITIALIZE COUNTERS AND FLAG
ICNVRG = 0
ITER = 0
ERROR = 1.0
FLAGNR = 0
C ---- ITERATION LOOP

1 ITER = ITER + 1

154

C ___-
C _-_
8
11
21
C ___
C _-_
C ___
C _--
30
C ___

CHECK FOR ITERMAX

IF (ITER .GT. ITRMAX) STOP

INITIALIZE THE GLOBAL STIFFNESS ARRAY AND FORCE VECTOR.

DO 8 I = 1, NDFMSH
GF(I) = 0.000
CONTINUE

IF (FLAGNR .EQ. 0) THEN
DO 21 I = 1, NDFMSH
DO 11 J = 1, BNDWTH
GSTIF(I,J) = 0.000
CONTINUE
CONTINUE
END IF

VERIFY IF WE ARE AT THE FIRST ITERATION OF THE FIRST LOAD STEP.

F_KO = 0
IF (ILSTEP .EQ. 1 .AND. ITER .EQ. 1) THEN
F_KO = 1

END IF

CALL CALC_T (NPEMAX, NLAMAX, NELMSH, NEMMAX,
NLAYER, NEQMAX, KX, KZ, ZlLAY, Z2LAY, AK,
FI_D, NODE, GU, T_LAY, KXO, KZO, F_K0,
LAMBDl, LAMBDZ)

call CR_FI(NLAMAX, NELMSH, NEMMAX,
NLAYER, Kz, ZlLAY, ZZLAY,
AK, FI_D, D_FI_N, D_FI)

BEGIN DO-LOOP ON NUMBER OF ELEMENTS TO CALCULATE THE ELEMENT
MATRICES AND TO ASSEMBLE THE GLOBAL MATRICES.

DO 50 ELMNUM = 1, NELMSH

DETERMINE THE LOCAL (ELEMENT) DATA.

DO 30 ENOD = 1, NODELM
GNOD = NODE(ELMNUM,ENOD)
EX(ENOD) = X(GNOD)
EZ(ENOD) = Z(GNOD)

CONTINUE

call stiff_T (EDFMAX, NPEMAX, NGPMAX, ndfelm,
NGPstf, ELMNUM,
DSF , GDSF , EX , ESTIF , EF,
GAUSPT, GAUSWT,
EZ, KX , KZ, ZlLAY, ZZLAY, NEMMAX,
NLAMAX, NLAYER,AK, FI_D, D_FI, SF2)

PRINT STIFFNESS MATRIX AND FORCE VECTOR FOR ELEMENT ONE.

155

 

40

IF (PRTSTF .EQ. 1 .AND. ELMNUM .EQ. 1 .AND. F_KO .EQ. 1) THEN
WRITE(OUNIT,1000)
DO 40 I = 1, NDFELM
WRITE(OUNIT,1010) I
WRITE(OUNIT,1020) (ESTIF(I,J), J = 1, NDFELM)
CONTINUE
WRITE(OUNIT,1030)
WRITE(OUNIT,1020) (EF(I), I = 1, NDFELM)
END IF

C --- ASSEMBLE THE GLOBAL MATRICES IN BANDED FORM.
CALL ASMBLE (NEQMAX, BNDMAX, EDFMAX, NEMMAX, NPEMAX, NODELM,
NDFNOD, NODE , ESTIF , EF , GSTIF , GF ,
ELMNUM)
C --- END OF LOOP TO COMPUTE GLOBAL STIFFNESS AND FORCE MATRICES.
50 CONTINUE
C --- PRINT GLOBAL STIFFNESS MATRIX AND FORCE VECTOR. -- VEDI PRECEDENTI
VERSIONI
C ----COMPUTE NEW DISPLACEMENTS VECTOR
DO 13 I = 1, NBCMAX
VSGD(I) = IVSGD(I)*(LAMBDA + DLAM)
13 CONTINUE
C --- IMPOSE PRIMARY BOUNDARY CONDITIONS.
CALL UBNDRY (NEQMAX, BNDMAX, NBCMAX, NDFNOD, NDFMSH, BNDWTH,
NSGD , ISGD , VSGD , GSTIF , GF )
C --- SOLVE THE SYSTEM OF EQUATIONS.
IRES = 0
CALL SOLVE (NEQMAX, BNDMAX, NDFMSH, BNDWTH, GSTIF, GF, IRES)
C --- COMPUTE THE INCREMENT.
DO 113 I = 1, NEQMAX
DELTA_D(I) = GF(I) - GU(I)
113 CONTINUE
C ----- CHECK FOR CONVERGENGE
GF2 = 0.0000
GU2 = 0.0000
DO 63 I = 1, NEQMAX
GF2 = GF2 + DELTA_D(I)*DELTA_D(I)
GU2 = GU2 + GU(I)*GU(I)
63 CONTINUE

IF (ITER .GT. 1) THEN

156

UINC = GF2
UOLD = GU2
ERROR = DSQRT(UINC/UOLD)
IF (ERROR .LT. EPSILN) THEN
ICNVRG = 1
END IF
END IF

C ----- UPDATE SOLUTION

DO 73 I = 1, NEQMAX
GU(I) = GF(I)
73 CONTINUE

IF (ICNVRG .EQ. 1) THEN
DO 83 I = 1, NEQMAX
GUPLS(I) = GU(I)
83 CONTINUE
END IF

C ----- PRINT VALUE OF T AT THE CENTER OF THE MESH OR NEAR
IF (F1ag_P .EQ. 0) THEN

WRITE(OUNIT,2000)
WRITE(OUNIT,2030)

end if
F1ag_P = 1
IF (PRTITR .EQ. 1 .and. ICNVRG .EQ. 1 ) THEN
CEN = (NELM_X+1)*(NELM_Z/2) + (NELM_X/2+1)
WRITE(OUNIT,2040) GU(CEN)
2000 FORMAT(' ',/3X,'SOLUTION AT THE CENTER OF THE MESH: ')
2030 FORMAT(' ',8X,'T (CENTER)')
2040 FORMAT(' ',5X, E14.8,6X)
END IF
C ---—INCREMENT LOAD

IF (ICNVRG .EQ. 1) THEN

C---- CONVERGENCE OBTAINED, NEXT LOAD STEP

LAMBDA = LAMBDA + DLAM
GOTO 500

ELSE

C---- CONVERGENCE NOT OBTAINED, NEXT ITERATION
GOTO 1
END IF

500 CONTINUE

1000 FORMAT(' ',//3X,'STIFFNESS MATRIX FOR ELEMENT ONE: ')

157

1010 FORMAT(' ',/5X,'ROW: ',I2,/)

1020 FORMAT(' ',5X,6E12.4)

1030 FORMAT(' ',//3X,'FORCE VECTOR FOR ELEMENT ONE: ',/)
1040 FORMAT(' ',//3X,'GLOBAL STIFFNESS MATRIX: ',/)

1050 FORMAT(' ',//3X,'GLOBAL FORCE VECTOR: ',/)

RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: CALC_T C
C PURPOSE: EVALUATE TEMPERATURE AT THE CENTER OF EACH LAYER C
C IN EACH ELEMENT C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE CALC_T(NPEMAX, NLAMAX, NELMSH, NEMMAX,
NLAYER, NEQMAX, KX, Kz, ZlLAY, ZZLAY, AK,
FI_D, NODE, GU, T_LAY, KXO, KZO, F_K0,
LAMBDl, LAMBDZ)

IMPLICIT LOGICAL (A-Z)

integer NPEMAX, NLAMAX, NELMSH,NEMMAX,
NEQMAX, NLAYER,
I, L, K1, K2, NODE, F_KO

double precision KX, KZ, ZlLAY, Z2LAY, AK, FI_D,
FI_N, FI_T, GU, T_LAY,
KXO, KZO, XIl, XI2, X3, LAMBDl, LAMBDZ

DIMENSION NLAYER(NEMMAX), NODE(NEMMAX,NPEMAX),
KX(nemmax,NLAMAX), KZ(nemmax,NLAMAX),
ZlLAY(nemmax,NLAMAX), ZZLAY(nemmax,NLAMAX),

AK(nemmax,NLAMAX), FI_D(nemmax),
GU(NEQMAX), T_LAY(nemmax,NLAMAX),

KXO(nemmax,NLAMAX), KZO(nemmax,NLAMAX),
LAMBDl(nemmax,NLAMAX), LAMBD2(nemmax,NLAMAX)

IF (F_K0 .EQ. 0) THEN
DO 10 I = 1, NELMSH
C ----- TEMPERATURE AT THE CENTER OF LAYER 1

FI_T = ( (z2LAY(1,1)- ZlLAY(I,1)) * 0.5 ) / FI_D(I)

x11 = 0.5

x12 = 0.5

T_LAY(I,1) = GU(NODE(I,1))*(1 - FI_T)*XI1 +
GU(NODE(I,2))*(1 - FI_T)*XI2 +
GU(NODE(I,3))*(FI_T)*XI2 +
GU(NODE(I,4))*(FI_T)*XI1

FI_T = 0.0000
DO 30 K1 = 2, NLAYER(I)

158

 

 

FI_N = (Z2LAY(I,K1) - ZlLAY(I,1)) - (ZZLAY(I,K1) - Z1LAY(I,K1))/2

x3 = FI_N
DO 40 K2 = 2, K1
FI_N = FI N + ( x3 - ( Z2LAY(I,K2-1)-Z1LAY(I,1) ) )*AK(I,K2-1)

40 CONTINUE

C ----- COMPUTE TEMPERATURE AT THE CENTER OF LAYER 2,3,. .,N
FI_T = FI_N / FI_D(I)
x11 = 0.5
x12 = 0.5

T_LAY(I,K1) = GU(NODE(I,1))*(1 - FI_T)*XIl +
GU(NODE(I,2))*(1 - FI_T)*X12 +
GU(NODE(I,3))*(FI_T)*XI2 +
GU(NODE(I,4))*(FI_T)*XI1

30 CONTINUE

10 CONTINUE

C---— EVALUATE THE NEW KX AND KZ

DO 101 I = 1, NELMSH
DO 151 L = 1, NLAYER(I)
KX(I,L)=KXO(I,L) + LAMBD1(I,L)*T_LAY(I,L)+
. LAMBD2(I,L)*T_LAY(I,L)*T_LAY(I,L)
KZ(I,L)=KZO(I,L) + LAMBD1(I,L)*T_LAY(I,L)+
. LAMBD2(I,L)*T_LAY(I,L)*T_LAY(I,L)
151 CONTINUE '
101 CONTINUE

ELSE

DO 200 I = 1, NELMSH
DO 250 L = 1, NLAYER(I)
KX(I,L) = KXO(I,L)
KZ(I,L) = KZO(I,L)
250 CONTINUE
200 CONTINUE

END IF

RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: ASMBLE C
C PURPOSE: ASSEMBLE THE GLOBAL MATRICES IN BANDED FORM. C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE ASMBLE (NEQMAX, BNDMAX, EDFMAX, NEMMAX, NPEMAX, NODELM,
NDFNOD, NODE , ESTIF , EF , GSTIF , GF ,
ELMNUM)

IMPLICIT LOGICAL (A-Z)

159

 

10
20
30
40

INTEGER NEQMAX, BNDMAX, EDFMAX, NEMMAX, NPEMAX, NODELM,
NDFNOD, NODE , ELMNUM, INOD , JNOD , IDOF ,
JDOF , ROW , COL , COL1 , EROW , ECOL

DOUBLE PRECISION ESTIF , EF , GSTIF , GF

DIMENSION NODE(NEMMAX,NPEMAX) , ESTIF(EDFMAX,EDFMAX), EF(EDFMAX),
GSTIF(NEQMAX,BNDMAX), GF(NEQMAX)

DO 40 INOD = 1, NODELM
ROW = NDFNOD * (NODE(ELMNUM,INOD) - 1)

DO 30 IDOF = 1, NDFNOD
ROW = ROW + 1
EROW = (INOD - 1) * NDFNOD + IDOF
GF(ROW) = GF(ROW) + EF(EROW)

DO 20 JNOD = 1, NODELM
COL1 = NDFNOD * (NODE(ELMNUM,JNOD) - 1)

DO 10 JDOF = 1, NDFNOD
COL = COL1 - ROW + JDOF + 1
ECOL = (JNOD - l) * NDFNOD + JDOF
IF (COL .GT. 0) THEN
GSTIF(ROW,COL) = GSTIF(ROW,COL) + ESTIF(EROW,ECOL)
END IF
CONTINUE
CONTINUE
CONTINUE
CONTINUE
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C

C
SUBROUTINE: CR_FI C
PURPOSE: EVALUATE FUNCTIONS FI AND ITS DERIVATIVES C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE CR_FI(NLAMAX, NELMSH, NEMMAX,
NLAYER, KZ, ZlLAY, ZZLAY,
AK, FI_D, D_FI_N, D_FI)

IMPLICIT LOGICAL (A-Z)

integer NLAMAX, NELMSH,NEMMAX,
NLAYER, I , K, L , N

double precision KZ, ZlLAY, Z2LAY, AK, FI_D,
H, R, D_FI_N, D_FI, SUM_AK

DIMENSION NLAYER(NEMMAX), KZ(nemmax,NLAMAX),
ZlLAY(nemmax,NLAMAX), ZZLAY(nemmax,NLAMAX),
AK(nemmax,NLAMAX), FI_D(nemmax), D_FI_N(nemmax,NLAMAX),
D_FI(nemmax,NLAMAX)

160

 

3O

4O

50

60

100

D0 100 I = 1, NELMSH

if (NLAYER(I) .gt. 1) then
SUM_AK = O
AK(I,1) = ( KZ(I,1)/KZ(I,2) - 1 )
SUM_AK = AK(I,1)

DO 30 K = 2, (NLAYER(I)—1)
AK(I,K) = ( KZ(I,K)/KZ(I,K+1) - 1 )*( 1 + SUM_AK)
SUM_AK = SUM_AK + AK(I,K)
CONTINUE
END IF
H = ( Z2LAY(I,NLAYER(I)) - ZlLAY(I,1) )
FI_D(I) = H
DO 40 R = 1, (NLAYER(I)-1)
FI_D(I) = FI_D(I) + ( H — (ZZLAY(I,R)- ZlLAY(I,1)) )*AK(I,R)
CONTINUE
D_FI_N(I,1) = 1
DO 50 L = 2, NLAYER(I)
D_FI_N(I,L) = D_FI_N(I,L-1) + AK(I,L—1)
CONTINUE
DO 60 N = 1, NLAYER(I)
D_FI(I,N) = D_FI_N(I,N) / FI_D(I)
CONTINUE
CONTINUE
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C
C
C
C

C

SUBROUTINE: FBNDRY C
PURPOSE: IMPOSE THE PRESCRIBED FORCE B.C.'S ON THE C
BANDED SYMMETRIC STIFFNESS MATRIX AND MODIFY C

THE FORCE VECTOR. C

C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE FBNDRY (NEQMAX, NBCMAX, NDFNOD, NSGF , ISGF , VSGF ,
GF )

IMPLICIT LOGICAL (A-Z)

INTEGER NEQMAX, NBCMAX, NDFNOD, NSGF , ISGF , IBC ,
IDOF

DOUBLE PRECISION VSGF , GF

DIMENSION ISGF(NBCMAX,2), VSGF(NBCMAX), GF(NEQMAX)

161

 

DO 10 IBC = 1, NSGF

IDOF = NDFNOD * (ISGF(IBC,1) - 1) + ISGF(IBC,2)

GF(IDOF) = GF(IDOF) + VSGF(IBC)

10 CONTINUE

RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C
C
C

C ___

10

20
C ___
C ___
C _-_
C __-

SUBROUTINE: GQUAD

PURPOSE: ASSIGN VALUES TO THE ARRAYS "GAUSPT" AND "GAUSWT"

USED IN GAUSSIAN QUADRATURE.

SUBROUTINE GQUAD (NGPMAX, GAUSPT, GAUSWT )

IMPLICIT LOGICAL (A-Z)

INTEGER NGPMAX, I , J
DOUBLE PRECISION GAUSPT, GAUSWT

DIMENSION GAUSPT(NGPMAX,NGPMAX),

INITIALIZE ARRAYS.

DO 20 I = l,NGPMAX
DO 10 J = l,NGPMAX
GAUSPT(I,J) = 0.000
GAUSWT(I,J) 0.000
CONTINUE
CONTINUE

GAUSSIAN QUADRATURE OF ORDER 1.

GAUSPT(1,1) = 0.000
GAUSWT(1,1) 2.000

GAUSSIAN QUADRATURE OF ORDER 2.

GAUSPT(1,2) - -1.000/SQRT(3.000)

GAUSPT(2,2) = -GAUSPT(1,2)
GAUSWT(1,2) = 1.000
GAUSWT(2,2) = GAUSWT(1,2)

GAUSSIAN QUADRATURE OF ORDER 3.

GAUSPT(1,3) = -SQRT(3.000/5.000)
GAUSPT(2,3) = 0.000

GAUSPT(3,3) = ~GAUSPT(1,3)
GAUSWT(1,3) = 5.000/9.000
GAUSWT(2,3) = 8.000/9.000
GAUSWT(3,3) = GAUSWT(1,3)

GAUSSIAN QUADRATURE OF ORDER 4.

GAUSPT(1,4) = -0.8611363ll6
GAUSPT(2,4) = -O.3399810436
GAUSPT(3,4) = -GAUSPT(2,4)

162

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

GAUSWT(NGPMAX,NGPMAX)

C
C
C
C
C
C

 

GAUSPT(4,4) = -GAUSPT(1,4)
GAUSWT(1,4) = 0.3478548451
GAUSWT(2,4) = 0.6521451549
GAUSWT(3,4) = GAUSWT(2,4)
GAUSWT(4,4) = GAUSWT(1,4)
RETURN

END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C

C
SUBROUTINE: INPELM C
PURPOSE: READ ELEMENT DATA. C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

1000
1010

SUBROUTINE INPELM (OUNIT , nodelm, ndfelm, ngpstf, ngpgrd,
ndfnod)

IMPLICIT LOGICAL (A-Z)

integer OUNIT , nodelm, ndfelm, ngpstf, ngpgrd, ndfnod
NDFNOD = l
NDFELM = 4
NGPSTF = 2
NGPGRD = 2
NODELM = 4

WRITE(OUNIT,1000)
WRITE(OUNIT,1010) NODELM, ndfnod, ndfelm, ngpstf, ngpgrd

FORMAT(' ',//3X,'ELEMENT DATA:')
FORMAT(' ',

//1OX,'NUMBER OF NODES PER ELEMENT, . ............. NODELM =',I4,
/10X,'NUMBER OF D.O.F. IN EACH NODE, ........... . NDFNOD =',I4,
/10X,'NUMBER OF D.O.F. IN EACH ELEMENT, ......... NDFELM =',I4,

//10X,'INTEGRATION RULE , ........................ NGPSTF =',I4,
/10X,'INTEGRATION POINTS , ...................... NGPGRD =',I4)

RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

000000

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
SUBROUTINE: INP_NL C
PURPOSE: READ OPTION PARAMETERS FOR NON LINEAR ANALYSIS C
FROM THE INPUT FILE. C
C
C
SUBROUTINE INP_NL (IUNIT , OUNIT , NLSTEP, DLAM, ITRMAX, EPSILN )
IMPLICIT LOGICAL (A-Z)

INTEGER IUNIT , OUNIT , NLSTEP, ITRMAX

DOUBLE PRECISION DLAM, EPSILN, NLS_TM

163

 

READ(IUNIT,*) NLSTEP, ITRMAX, EPSILN
NLS_TM = NLSTEP
DLAM = 1/NLS_TM

WRITE(OUNIT,1000)
WRITE(OUNIT,1010) NLSTEP, DLAM, ITRMAX, EPSILN

1000 FORMAT(' ',//3X,'NON LINEAR ANALYSIS DATA:')
1010 FORMAT(' ',

.//10X,'NUMBER OF LOAD STEP, .................... NLSTEP =',I6,
/10X,'LOAD INCREMENT SCALE, ................... DLAM =',E10.4,
/10X,'MAX. NUMBER OF ITERATIONS, .............. ITRMAX =',I6,
/10X,'TOLLERANCE, .... ......................... EPSILN =',E10.4)

RETURN

END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: INLAYR C
C PURPOSE: INPUT ELEMENT LAYERS C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE INLAYR(IUNIT , OUNIT , NELMSH, NLAYER, NLAMAX,NEMMAX,
KXO, KZO, ZlLAY, Z2LAY, LAMBDl, LAMBD2, PRTLAY,
NELM_X, NELM_Z, INPLAY)

IMPLICIT LOGICAL (A-Z)

INTEGER I, J, E, K, IUNIT, OUNIT, NELMSH, NLAYER, NLAMAX,
NEMMAX, PRTLAY, INPLAY, NELM_X, NELM_Z, ASSLAY,
NLATMP

DOUBLE PRECISION KXO, KZO, ZlLAY, ZZLAY, LAMBDl, LAMBD2,
KXOTMP, KZOTMP, LBlTMP, LB2TMP,
ZlLTMP, ZZLTMP

DIMENSION NLAYER(NEMMAX),
KXO (nemmax,NLAMAX) , KZO (nemmax,NLAMAX) ,
Z1LAY(nemmax,NLAMAX) , Z2LAY (nemmax,NLAMAX) ,
LAMBDl (nemmax,NLAMAX) , LAMBD2 (nemmax,NLAMAX)

IF (INPLAY .EQ. 1) THEN
E = 0
ASSLAY = 0
DO 80 I=1, NELM_Z
READ(IUNIT,*) NLATMP
DO 70 J = 1, NLATMP
E = 0
READ(IUNIT,*) KXOTMP, KZOTMP, LBlTMP, LBZTMP, ZlLTMP, zzLTMP
E = 1 + (I-1)*NELM_X
DO 60 K = 1, NELM_X
NLAYER(E) = NLATMP
KXO(E,J) = KXOTMP
KZO(E,J) = KZOTMP

164

LAMBD1(E,J) = LBlTMP
LAMBD2(E,J) = LBZTMP
Z1LAY(E,J) = ZlLTMP

ZZLAY(E,J) = Z2LTMP
E = E + 1
60 CONTINUE
70 CONTINUE

80 CONTINUE
ELSE

DO 180 I=1, NELMSH
READ(IUNIT,*) NLAYER(I)
DO 170 J=1, NLAYER(I)
READ(IUNIT,*) KXO(I,J), KZO(I,J), LAMBD1(I,J), LAMBD2(I,J),
. ZlLAY(I,J), Z2LAY(I,J)
170 CONTINUE
180 CONTINUE

END IF

IF (PRTLAY .EQ. 1) THEN
WRITE(OUNIT,1000)

DO 100 I=1, NELMSH

WRITE(OUNIT,1010)
WRITE(OUNIT,1020) I, NLAYER(I)
WRITE(OUNIT,1030)
DO 90 J=1, NLAYER(I)
WRITE(OUNIT,1040) J,KXO(I,J),KZO(I,J),ZlLAY(I,J),ZZLAY(I,J)

90 CONTINUE
100 CONTINUE
END IF
C --- FORMAT STATEMENTS.
1000 FORMAT(' ,///3X,'NUMBERS AND CHARACTERISTICS OF LAYERS: ')

1010 FORMAT('
1020 FORMAT('
1030 FORMAT('

,//6X,'ELEMENT',3X,'N LAYER',/)

,9X,13,7X,I3)

,/6X,'N LAYER',5X,'KXO',12X,'KZO',12X,'Z1LAY',9X,
. 'ZZLAY',5X,/)

1040 FORMAT(' ',9X,I3,4X,E10.4,4X,E10.4,4X,E10.4,4X,E10.4)

RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: INPMSH C
C PURPOSE: INPUT MESH PARAMETERS AND FORM THE MESH. C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE INPMSH (IUNIT , OUNIT , NEMMAX, NPEMAX, NNMMAX, nodelm,

ndfnod,DOM_X, DOM_Z, NELM_X, NELM_Z,
NELMSH, NDFMSH, nodmsh, node, BNDWTH, X, Z)

165

 

IMPLICIT LOGICAL (A-Z)

integer IUNIT , OUNIT , NEMMAX, NPEMAX, NNMMAX, nodelm,
ndfnod, NELM_X, NELM_Z, nodmsh, BNDWTH,
NDFMSH, node, IMESH, EORDER, NXINCR, NZINCR,
I , J , N, NELMSH, NW, NPE

DOUBLE PRECISION X , Z, DOM_X, DOM_Z, DX, DZ

DIMENSION NODE(NEMMAX,NPEMAX), X(NNMMAX), Z(NNMMAX),
DX(NNMMAX) , DZ(NNMMAX)

READ(IUNIT,*) IMESH, DOM_X, DOM_Z, NELM_X, NELM_Z

NODMSH = (NELM_X + l)*(NELM_Z + 1)
NELMSH = NELM_X * NELM_Z

IF (IMESH .EQ. 0) THEN
DO 10 I = 1, NELMSH
READ(IUNIT,*) (NODE(I,J), J=1,NODELM)
10 CONTINUE
DO 20 I = 1, NODMSH
READ(IUNIT,*) X(I), Z(I)
20 CONTINUE

ELSE

EORDER = 1

NXINCR = NELM_X * EORDER

NZINCR = NELM_Z * EORDER

READ(IUNIT,*) (DX(I), I = 1, NXINCR)

READ(IUNIT,*) (DZ(I), I = 1, NZINCR)

NPE = NODELM

CALL MESH (NPEMAX, NEMMAX, NNMMAX, EORDER, NODE , NPE ,

NODMSH, NELMSH, NELM_X, NELM_Z, X , Z ,
. DX , DZ)
END IF
NDFMSH = NDFNOD * NODMSH
BNDWTH = 0
DO 1 N = l,NELMSH
DO 2 I = l,NODELM
DO 3 J = l,NODELM
NW = ( ABS( NODE(N,I) - NODE(N,J) ) + 1 )

IF (BNDWTH .LT. NW) THEN
BNDWTH = NW

END IF
3 CONTINUE
2 CONTINUE
1 CONTINUE

WRITE(OUNIT,1000)
WRITE(OUNIT,1010) NELM_X, NELM_Z, NODMSH, NDFMSH, BNDWTH
WRITE(OUNIT,1020)
WRITE(OUNIT,1030)
DO 90 N = l,NELMSH
WRITE(OUNIT,1040) N, (NODE(N,I), I=1,NODELM)
90 CONTINUE

166

1000 FORMAT(' ',//3X,'MESH DATA:')
1010 FORMAT(' ',

//10X,'N. OF ELEMENTS IN THE MESH IN X DIRECTION, .NELMSH =',I4,
/10X,'N. OF ELEMENTS IN THE MESH IN Z DIRECTION, .NELMSH =',I4,
/10X,'NUMBER OF NODES IN THE MESH, .............. NODMSH =',I4,
/10X,'NUMBER OF D.O.F. IN THE MESH, ............. NDFMSH =',I4,
/10X,'HALF BAND WIDTH OF GLOBAL STIFFNESS MATRIX, BNDWTH =',I4)

1020 FORMAT(' ',//10X,'BOOLEAN (CONNECTIVITY) MATRIX: ')

1030 FORMAT(' ',/10X,'ELEMENT NO.',6X,'NODE(1,2,3,4)',/)
1040 FORMAT(' ',11X,I3,10X,9(I3,1X))

RETURN

END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C SUBROUTINE: MESH

C PURPOSE: GENERATE THE CONNECTIVITY ARRAY NODE(I,J), THE

C COORDINATES X(I), Z(I), AND MESH INFORMATION

C (NODMSH, NELMSH, NODELM) FOR RECTANGULAR DOMAINS.
C THE DOMAIN IS DIVIDED INTO RECTANGULAR ELEMENTS

C WITH FOUR, EIGHT, OR NINE NODES (NUMELX BY NUMELY
C NONUNIFORM MESH IN GENERAL).

C
C

000000000

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE MESH (NPEMAX, NEMMAX, NNMMAX, IEL , NOD , NPE,
NNM , NEM , NX , NZ , X , Z ,
DX , DZ )

IMPLICIT DOUBLE PRECISION (A-H,O-Z)

INTEGER NPEMAX, NEMMAX, NNMMAX, IEL , NOD , NPE ,
NNM , NEM , Nx , Nz

DOUBLE PRECISION x , z , DX , DZ

DIMENSION X(NNMMAX) , Z(NNMMAX) , NOD(NEMMAX,NPEMAX),

DX(NNMMAX), DZ(NNMMAX)

NEX1=NX+1

NEZl=NZ+1

NXX=IEL*NX

NZZ=IEL*NZ

NXX1=NXX+1

NZZl=NZZ+1

NEM=NX*NZ
NNM=NXX1*NZZ1-(IEL-1)*NX*NZ
IF(NPE.EQ.9)NNM=NXX1*NZZ1

C --- GENERATE THE CONNECTIVITY ARRAY NOD(I,J).

K0=0

IF(NPE.EQ.9)K0=1

NOD(1,1) = 1

NOD(1,2) = IEL+1

NOD(1,3) = NXX1+(IEL-1)*NEX1+IEL+1

167

200

210
220

230

240

250

260

270

280

290
300

310

IF (NPE .EQ. 9)NOD(1,3)=4*NX+5
NOD(1,4) = NOD(1,3) - IEL
IF(NPE .EQ. 4)GO TO 200

NOD(1,5) = 2
NOD(1,6) = NXXl + (NPE-6)
NOD(1,7) = NOD(1,3) - 1

NOD(1,8) = NXX1+1
IF (NPE .EQ. 9)NOD(1,9)=NXX1+2
IF(NZ .EQ. l)GOTO 230
M = 1
DO 220 N = 2,NZ
L = (N-1)*NX + 1
DO 210 I = l,NPE
NOD(L,I) = NOD(M,I)+NXX1+(IEL-1)*NEX1+KO*NX
M=L

IF(NX .EQ .1)GO TO 270

D0 260 NI = 2,NX
DO 240 I = l,NPE
K1 = IEL
IF(I .EQ. 6 .OR. I .EQ. 8)K1=1+K0
NOD(NI,I) NOD(NI-1,I)+K1
M = NI
DO 260 NJ = 2,NZ
L = (NJ-1)*NX+NI
DO 250 J = l,NPE
NOD(L,J) = NOD(M,J)+NXX1+(IEL-1)*NEX1+K0*NX
M = L

GENERATE THE COORDINATES X(I) AND Z(I)

DX(NXX1)=0.0

DZ(NZZl)=0.0

ZC=0.0

IF (NPE .EQ. 9) GOTO 310
D0 300 NI = 1, NEZl

I = (NXX1+(IEL-1)*NEX1)*(NI—1)+1
J = (NI-1)*IEL+1

X(I) = 0.0

Z(I) = ZC

DO 280 NJ = 1,NXX

I=I+1

X(I) = X(I-1)+DX(NJ)

Z(I) = ZC

IF(NI.GT.NZ .OR. IEL.EQ.1)GO TO 300
J = J+1

ZC = ZC+DZ(J—l)

I = I+1

X(I) = 0.0

Z(I) = ZC

DO 290 II = 1, NX

K = 2*II-1

I = I+1

X(I) = X(I-1)+DX(K)+DX(K+1)
Z(I) = ZC

ZC = ZC+DZ(J)

RETURN

DO 330 NI=1,NZZI
I=NXX1*(NI-l)

168

320
330

XC=0.0
DO 320 NJ=1,NXX1

I = I+1
X(I) = XC
Z(I) = ZC

XC = XC + DX(NJ)
ZC = ZC + DZ(NI)
RETURN

END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C

SUBROUTINE: INPNBC
PURPOSE: INPUT NODAL BOUNDARY CONDITIONS.

C
C
C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

10

20

30

SUBROUTINE INPNBC (IUNIT , OUNIT , NBCMAX, NNMMAX, NODMSH,
NSGD , NSGF , ISGD , ISGF , X, Z ,
VSGF, FLAG_D)

IMPLICIT LOGICAL (A-Z)

INTEGER IUNIT , OUNIT , NBCMAX, NNMMAX, NODMSH,
NSGD , NSGF , ISGD , ISGF , FLAG_D,
I , J , K , IBC , NOD ,

DOUBLE PRECISION X , Z, VSGD , VSGF

DIMENSION ISGD(NBCMAX,2), ISGF(NBCMAX,2), IBC(4),

ndfnod,
vsgd ,

ndfnod,

DOF

VSGD(NBCMAX) , VSGF(NBCMAX) , X(NNMMAX),

Z(NNMMAX)

READ(IUNIT,*) NSGD, NSGF

FLAG_D = 1

IF (NSGD .GT. 0) THEN
READ(IUNIT,*) ((ISGD(I,J),J=1,2),I=1,NSGD)
READ(IUNIT,*) (VSGD(I),I=1,NSGD)

END IF

IF (NSGF .GT. 0) THEN
FLAG_D = 0
READ(IUNIT,*) ((ISGF(I,J),J=1,2),I=1,NSGF)
READ(IUNIT,*) (VSGF(I),I=1,NSGF)

END IF

WRITE(OUNIT,1000)
WRITE(OUNIT,1010) NSGD
WRITE(OUNIT,1020)
DO 30 I = l,NODMSH
DO 10 K = l,NDFNOD
IBC(K) = 0
CONTINUE
DO 20 J = 1,NSGD
NOD = ISGD(J,1)
DOF = ISGD(J,2)
IF (NOD .EQ. I) IBC(DOF) = DOF
CONTINUE
WRITE(OUNIT,1030) I, X(I), Z(I), (IBC(K), K = 1, NDFNOD)
CONTINUE

169

4O

50

60

C ___

1000
1010
1020

1030

1035
1040
1050

1060
1070
2050
2060

WRITE(OUNIT,1035)
WRITE(OUNIT,1070) (VSGD(I), I = l, NSGD)

WRITE(OUNIT,2050)

IF (NSGF .NE. 0) THEN
WRITE(OUNIT,1040) NSGF
WRITE(OUNIT,1050)
DO 60 I = l,NODMSH
DO 40 K = l,NDFNOD
IBC(K) = 0
CONTINUE
DO 50 J = 1,NSGF
NOD = ISGF(J,1)
DOF = ISGF(J,2)
IF (NOD .EQ. I) IBC(DOF) = DOF
CONTINUE
WRITE(OUNIT,1030) I, X(I), Z(I), (IBC(K), K = 1, NDFNOD)
CONTINUE
WRITE(OUNIT,1060)
WRITE(OUNIT,1070) (VSGF(I), I = l, NSGF)
WRITE(OUNIT,2060)
END IF

FORMAT STATEMENTS.

FORMAT(' ',//3X,'NODAL BOUNDARY CONDITIONS: ',/)

FORMAT(' ',9X,'NUMBER OF SPECIFIED TEMPERATURE D.O.F =',I4)

FORMAT(' ',/10x,'NODE',4X,'X-COORD.',6X,'z-COORD.'5x,
'SPECIFIED TEMPERATURE D.O.F.',/)

FORMAT(' ',9X,I3,4X,E10.4,4X,E10.4,3X,12(I2,1X))

FORMAT(' ',/10X,'VALUES OF THE SPECIFIED TEMPERATURE D.O.F.:')

FORMAT(' ',lOX,'NUMBER OF SPECIFIED GENERALIZED FORCES =',I4)

FORMAT(' ',/10X,'NODE',4X,'X-COORD.',4X,'Z-COORD.',5X,
'SPECIFIED GENERALIZED FORCE',/)

FORMAT(' ',/10X,'VALUES OF THE SPECIFIED GENERALIZED FORCES:')
FORMAT(' ',/10X,6(E10.4,2X))

FORMAT(' ',/)

FORMAT(' ',/)

RETURN

END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C

SUBROUTINE: INPOPT
PURPOSE: READ OPTION PARAMETERS FROM THE INPUT FILE.

C
C
C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE INPOPT (IUNIT, OUNIT, PRTSTF, PRTGRD, PRTLAY, INPLAY,

PRTITR)

IMPLICIT LOGICAL (A-Z)

INTEGER IUNIT , OUNIT , PRTSTF, PRTGRD, PRTLAY, INPLAY, PRTITR

CHARACTER*72 TITLE

READ(IUNIT,1000) TITLE
READ(IUNIT,*) PRTSTF, PRTLAY, INPLAY, PRTITR

PRTGRD = 0

170

WRITE(OUNIT,1010)
WRITE(OUNIT,1020)

1000 FORMAT(A72)

TITLE

"THERMNL VERS.

1010 FORMAT(' ',//3X,
'RESULTS OF
1020 FORMAT(' 'I//3X,A72)
RETURN
END

1.0"

- WRITTEN BY ANTONIO PANTANO')

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: PRTGRD C
C PURPOSE: WRITE THE SOLUTION GRADIENTS TO THE OUTPUT FILE. C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

subroutine prtdu

IMPLICIT LOGICAL

integer

double precision

dimension

WRITE(OUNIT,1000)
WRITE(OUNIT,1010)

C ___
DO 50 ELMNUM = 1,
C --- DETERMINE THE LOCAL (ELEMENT)
C --- COORDINATES.
gdof = (node(elmnum,1) —

do 10 edof = 1, ndfelm
eu(edof) = gu(gdof+edof)
10 continue
do 20 enod = 1, nodelm
gnod = node(elmnum,enod)
ex(enod) = x(gnod)
20 continue

(ounit , npemax,
ngpmax, nodelm,
node ,
gaUSPtI 9“ I
gdsf , a0 ,
(A—Z)
ounit , npemax,
ngpmax, nodelm,
node ,
elmnum, idof ,
gdof
gaUSPtI gu I
gdsf , a0 ,
xc , dudx ,
gu(neqmax) ,
sf(npemax) ,
a0(nemmax) ,
x(nnmmax) ,

neqmax, nnmmax,
ndfnod, ndfelm,
an I ex I
al I a2 I
neqmax, nnmmax,
ndfnod, ndfelm,
igp , enod ,
eu 1 ex I
a1 , a2 ,
adudx , a ,
eu(edfmax) I
gdsf(npemax) ,
a1(nemmax) ,

gauspt(ngpmax,ngpmax),

NELMSH

1)

* ndfnod

171

BEGIN LOOP OVER NUMBER OF ELEMENTS IN THE MESH.

edfmax, nemmax,
nelmsh, ngpgrd,
sf ,
x )
edfmax, nemmax,
nelmsh, ngpgrd,
gnod , edof ,
sf ,
X I
xi , elenth
ex(npemax) ,
a2(nemmax) ,

SOLUTION VECTOR AND NODAL

node(nemmax,npemax)

elenth = ex(nodelm) - ex(l)
C —-- BEGIN LOOP ON NUMBER OF GAUSS POINTS (NGPGRD).

DO 40 IGP = 1, NGPGRD
XI = GAUSPT(IGP,NGPGRD)

C --- COMPUTE THE SOLUTION GRADIENTS AT THE GAUSS POINTS.
CA CALL SHPL1D (NPEMAX, NODELM, XI , ELENTH, SF , DSF ,
CA . GDSF , JACOBN)

DUDX = 0.000
XC = 0.000
DO 30 IDOF = 1, NODELM
DUDX = DUDX + EU(IDOF) * GDSF(IDOF)
XC = XC + EX(IDOF) * SF(IDOF)
30 CONTINUE

A
ADUDX

A0(elmnum) + Al(elmnum)*XC + A2(e1mnum)*XC*XC
A*DUDX

WRITE(OUNIT,1030) ELMNUM, XC, DUDX, ADUDX

40 CONTINUE
50 CONTINUE

C --- FORMAT STATEMENTS.
1000 FORMAT(' ',/3X,'SOLUTION GRADIENTS: ',/)
1010 FORMAT(' ',/1OX,'ELEMENT NO. X-COORD. DU/DX A*DU/DX',/)
1030 FORMAT(' ', 11X,13,6X,4(E10.4,2X))
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: PRTGU C
C PURPOSE: WRITE THE SOLUTION VECTOR TO THE OUTPUT FILE. C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

subroutine prtgu (OUNIT , NEQMAX, NNMMAX, NDFNOD, NODMSH, GU ,
X, Z )

IMPLICIT LOGICAL (A-Z)

INTEGER OUNIT , NEQMAX, NNMMAX, NDFNOD, NODMSH,
U]. I 02 f I I J

DOUBLE PRECISION GU , X, Z

DIMENSION GU(NEQMAX) , X(NNMMAX), Z(NNMMAX)

WRITE(OUNIT,1000)
C --- OUTPUT OF SOLUTION VECTOR.

WRITE(OUNIT,1010)
U2 = 0

172

1000
1010
1020

DO 20 I = 1, NODMSH

U1 = 02 + 1

U2 = U2 + ndiOD

WRITE(OUNIT,1020) I, X(I), Z(I), (GU(J), J = U1, 02)
CONTINUE

FORMAT STATEMENTS.

FORMAT(' ',///3x,'FINAL SOLUTION OF THE NON LINEAR ANALYSIS: ',/)
FORMAT(' ',lOX,'NODE',7X,'X-COORD.',5X,'Z—COORD. '6X, 'T',/)
FORMAT(' ', 10x,I3,6x,E10.4,2x,E10.4,2x,E11.6)

RETURN

END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C
C
C
C
C

C ___

C _._—

10

C

SUBROUTINE: SHPLlD_T C
PURPOSE: EVALUATE THE l-D LAGRANGIAN INTERPOLATION FUNCTIONS C
AND THEIR GLOBAL DERIVATIVES AT THE GAUSS POINTS. C

C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE SHPLlD_T (NPEMAX, XI , ELENTH_X, DSF,
GDSF , JACOBN, SF2, IGP)

IMPLICIT LOGICAL (A-Z)

INTEGER NPEMAX, I, IGP

DOUBLE PRECISION XI , ELENTH_X, DSF , GDSF , JACOBN,
SF2

DIMENSION DSF(NPEMAX), GDSF(NPEMAX),
SF2(2,2)

LINEAR LAGRANGE INTERPOLATION FUNCTIONS FOR 2-NODED ELEMENTS.

SF2(1,IGP) = 0.5000 * (1.000 - XI)
SF2(2,IGP) = 0.5000 * (1.000 + XI)
DSF(I) = -0.5000

DSF(2) = 0.5000
COMPUTE THE GLOBAL DERIVATIVES OF SF(I).
JACOBN = ELENTH_X * 0.5000

DO 10 I = 1, 2
GDSF(I) = DSF(I) / JACOBN
CONTINUE
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C
C

C
SOLVE: C

173

THIS PROGRAM SOLVES A BANDED SYMMETRIC SYSTEM OF EQUATIONS. C
THE BANDED MATRIX IS INPUT THROUGH BAND(NEQNS,NBW), AND C
RHS IS THE RIGHT HAND SIDE (FORCE VECTOR) OF THE EQUATION. C
NEQNS IS THE NO. OF EQUATIONS AND NBW IS THE HALF BAND WIDTH. C
IN RESOLVING, IRES .GT. 0, LHS ELIMINATION IS SKIPPED. C
C
C

0000000

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE SOLVE (NRM , NCM , NEQNS , NBW , BAND, RHS, IRES)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)

INTEGER NRM , NCM , NEQNS , NBW , IRES
DOUBLE PRECISION BAND , RHS

DIMENSION BAND(NRM,NCM), RHS(NRM)

MEQNS=NEQNS-1

IF (IRES.GT.0) GO TO 40

DO 30 NPIV=1,MEQNS

NPIVOT=NPIV+1

LSTSUB=NPIV+NBW-l

IF (LSTSUB.GT.NEQNS) LSTSUB=NEQNS
DO 20 NROW=NPIVOT,LSTSUB

C --- INVERT ROWS AND COLUMNS FOR ROW FACTOR.
NCOL=NROW-NPIV+1

FACTOR=BAND(NPIV,NCOL)/BAND(NPIV,1)
DO 10 NCOL=NROW,LSTSUB
ICOL=NCOL-NROW+1
JCOL=NCOL-NPIV+1
10 BAND(NROW,ICOL)=BAND(NROW,ICOL)-FACTOR*BAND(NPIV,JCOL)
20 RHS(NROW)=RHS(NROW)-FACTOR*RHS(NPIV)
30 CONTINUE
GO TO 70
40 DO 60 NPIV=1,MEQNS
NPIVOT=NPIV+1
LSTSUB=NPIV+NBW-1
IF (LSTSUB.GT.NEQNS) LSTSUB=NEQNS
DO 50 NROW=NPIVOT,LSTSUB
NCOL=NROW-NPIV+1
FACTOR=BAND(NPIV,NCOL)/BAND(NPIV,1)
50 RHS(NROW)=RHS(NROW)-FACTOR*RHS(NPIV)
6O CONTINUE

C --- BACK SUBSTITUTION.

70 DO 90 IJK=2,NEQNS
NPIV=NEQNS-IJK+2
RHS(NPIV)=RHS(NPIV)/BAND(NPIV,1)
LSTSUB=NPIV—NBW+1
IF (LSTSUB.LT.1) LSTSUB=1
NPIVOT=NPIV-1
DO 80 JKI=LSTSUB,NPIVOT
NROW=NPIVOT-JKI+LSTSUB
NCOL=NPIV-NROW+1
FACTOR=BAND(NROW,NCOL)

80 RHS(NROW)=RHS(NROW)-FACTOR*RHS(NPIV)

174

90 CONTINUE
RHS(1)=RHS(l)/BAND(1,1)
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE: STIFF_T

C
C
C PURPOSE: COMPUTE THE ELEMENTAL STIFFNESS AND FORCE ARRAYS
C
C
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE STIFF_T(EDFMAX, NPEMAX, NGPMAX, ndfelm,
NGPstf, ELMNUM,
DSF, GDSF, EX, ESTIF, EF, GAUSPT,

GAUSWT,

EZ, KX , KZ, ZlLAY, Z2LAY, NEMMAX,

NLAMAX, NLAYER,AK, FI_D, D_FI, SF2)

IMPLICIT LOGICAL (A-Z)

integer EDFMAX, NPEMAX, NGPMAX, ndfelm,
NGPstf, I , J I igp.
NLAMAX, NLAYER, NEMMAX, ELMNUM, K

double precision DSF , GDSF, EX, ESTIF, EF,
GAUSPT, GAUSWT,
xi , const , jacobn,
ELENTH_X, ELENTH_Z,
EZ, KX , KZ, ZlLAY, Z2LAY,
AK, FI_D, D_FI, SF2, K_LAY,
FI, D_FI_F, FI_N, H_LAYER, FXI, DXI

dimension dsf(npemax), gdsf(npemax),
EX(NPEMAX), ef(edfmax),

GAUSPT(NGPMAX,NGPMAX), GAUSWT(NGPMAX,NGPMAX)

I

ESTIF(EDFMAX,EDFMAX), EZ(NPEMAX), NLAYER(NEMMAX),

KX(nemmax,NLAMAX), KZ(nemmax,NLAMAX),
ZlLAY(nemmax,NLAMAX), ZZLAY(nemmax,NLAMAX),
AK(nemmax,NLAMAX), FI_D(nemmax),

D FI(nemmax,NLAMAX), SF2(2,2), K_LAY(4,4,NLAMAX),

FI(2,3), D_FI_F(2), FXI(2,3), DXI(2)

ELENTH_X = EX(2) - EX(l)
ELENTH_Z = EZ(4) - EZ(l)

C --- INITIALIZE THE ARRAYS.

DO 2 I= 1, 4
DO 3 J= 1, 4
DO 4 K= 1, NLAYER(ELMNUM)
K_LAY(I,J,K) = 0 000

4 CONTINUE
3 CONTINUE
2 CONTINUE

DO 20 I = 1, NDFELM

175

C
C
C
C
C
C

EF(I) = 0.000
DO 10 J = 1, NDFELM
ESTIF(I,J) = 0.000
10 CONTINUE
20 CONTINUE

C --- BEGIN GAUSS QUADRATURE LOOP FOR FULL INTEGRATION.
DO 21 IGP = 1, NGPstf

XI = GAUSPT(IGP,NGPstf)

CALL SHPLlD_T (NPEMAX, XI, ELENTH_X, DSF,
GDSF , JACOBN, SF2, IGP)

CONST = JACOBN*GAUSWT(IGP,NGPstf)
21 CONTINUE

DO 22 I= 1, 4
DO 23 J= l, 4
DO 24 K= 1, NLAYER(ELMNUM)

CALL FI_AT_N(I, J, K, ELMNUM, ZlLAY, z2LAY,AK,
NEMMAX, NLAMAX, npemax, FI_D, D_FI,
D_FI_F, FI, FI_N, FXI, DXI, SF2, GDSF)

H_LAYER = ( Z2LAY(ELMNUM,K) - ZlLAY(ELMNUM,K) ) / 2

K_LAY(I,J,K)= H_LAYER/3 * (
JACOBN * KX(ELMNUM,K) * ( 2*DXI(1)*DXI(2) ) *
( FI(1,1)*FI(2,1) + 4*FI(1,2)*FI(2,2) + FI(1,3)*FI(2,3) )+
6 * JACOBN * KZ(ELMNUM,K) *
( FXI(1,1)*FXI(2,1) + FXI(1,2)*FXI(2,2) )*
( D_FI_F(l) * D_FI_F(2) ) )

ESTIF(I,J) = ESTIF(I,J) + K_LAY(I,J,K)
24 CONTINUE
23 CONTINUE

22 CONTINUE

50 CONTINUE

RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: UBNDRY C
C PURPOSE: IMPOSE THE PRESCRIBED DISPLACEMENT B.C.'S ON THE C
C BANDED SYMMETRIC STIFFNESS MATRIX AND MODIFY C
C THE FORCE VECTOR. C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE UBNDRY (NEQMAX, BNDMAX, NBCMAX, NDFNOD, NDFMSH, BNDWTH,
NSGD , ISGD , VSGD , GSTIF , GF )

176

IMPLICIT LOGICAL (A-Z)

INTEGER NEQMAX, BNDMAX, NBCMAX, NDFNOD, NDFMSH, BNDWTH,
NSGD , ISGD , IBC , IDOF , BWMl , ROW ,
COL , K

DOUBLE PRECISION VSGD , GSTIF , GF , UVALUE

DIMENSION ISGD(NBCMAX,2), VSGD(NBCMAX), GSTIF(NEQMAX,BNDMAX),
GF(NEQMAX)

DO 40 IBC = 1, NSGD

IDOF = NDFNOD * (ISGD(IBC,1) - 1) + ISGD(IBC,2)
UVALUE = VSGD(IBC)

BWMl = BNDWTH - 1

ROW = IDOF - BNDWTH

DO 20 K = l, BWMl

ROW = ROW + 1

IF (ROW .GT. 0) THEN
COL = IDOF - ROW + l
GF(ROW) = GF(ROW) - GSTIF(ROW,COL) * UVALUE
GSTIF(ROW,COL) = 0.000

END IF

20 CONTINUE

GSTIF(IDOF,1) = 1.000
GF(IDOF) = UVALUE
ROW = IDOF

DO 30 K = 2, BNDWTH
ROW = ROW + 1
IF (ROW .LE. NDFMSH) THEN
GF(ROW) = GF(ROW) - GSTIF(IDOF,K) * UVALUE
GSTIF(IDOF,K) = 0.000
END IF
30 CONTINUE
40 CONTINUE
RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: FI_AT_N C
c PURPOSE: CHOOSE THE FUNCTION FI AND XI TO BE USED IN THE C
C EVALUATION OF THE ACTUAL K_LAY C
c
C

C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE FI_AT_N (I,J,K, ELMNUM,ZlLAY,ZZLAY,AK,
NEMMAX, NLAMAX, npemax, FI_D, D_FI,
D_FI_F, FI, FI_N, FXI, DXI, SF2, gdsf)
IMPLICIT LOGICAL (A-Z)

integer I, J, K, R, Q,
NLAMAX, NEMMAX, ELMNUM, npemax

177

double precision ZlLAY, ZZLAY,
AK, FI_D, D_FI,
FI, D_FI_F, FI_N, ZZ, FXI, DXI,
gdsf,SF2

dimension ZlLAY(nemmax,NLAMAX), ZZLAY(nemmax,NLAMAX),
AK(nemmax,NLAMAX), FI_D(nemmax),
D_FI(nemmax,NLAMAX), gdsf(npemax), SF2(2,2),
FI(2,3), D_FI_F(2), FXI(2,3), DXI(2)

DO 90 R = 1,2

IF (I .EQ. 1 .OR. I .EQ. 4) THEN

DXI(l) = gdsf(l)
FXI(1,R) = SF2(1,R)
ELSE

DXI(l) = gdsf(2)
FXI(1,R) = SF2(2,R)
END IF

IF (J .EQ. 1 .OR. J .EQ. 4) THEN
DXI(2) = gdsf(l)
FXI(2,R) = SF2(1,R)
ELSE
DXI(2) = gdsf(2)
FXI(2,R) = SF2(2,R)
END IF

90 CONTINUE
DO 110 R = 1,3

FI_N = o 000
IF (R .EQ. 1) THEN
zz= ZlLAY(ELMNUM,K)-Z1LAY(ELMNUM,1)
ELSE IF (R .EQ. 2) THEN
zz=( (ZZLAY(ELMNUM,K)-ZlLAY(ELMNUM,1)) +

(ZlLAY(ELMNUM,K)-ZlLAY(ELMNUM,1)) )/2

ELSE IF (R .EQ. 3) THEN
zz= Z2LAY(ELMNUM,K)-ZlLAY(ELMNUM,l)

END IF
FI_N = zz
DO 100 Q = 2, K
FI_N = FI_N +(zz - (ZZLAY(ELMNUM,Q-1)-Z1LAY(ELMNUM,l)))
*AK(ELMNUM,Q-1)
100 CONTINUE

IF (I .LT. 3) THEN

D_FI_F(I) = —D_FI(ELMNUM,K)
FI(1,R)= (1 - FI_N/FI_D(ELMNUM) )
ELSE

D_FI_F(I) = D_FI(ELMNUM,K)

FI(1,R) = FI_N/FI_D(ELMNUM)
END IF

IF (J .LT. 3) THEN
D_FI_F(2) = —D_FI(ELMNUM,K)
FI(2,R) = (1 - FI_N/FI_D(ELMNUM) )

178

ELSE

D_FI_F(2) = D_FI(ELMNUM,K)
FI(2,R) = FI_N/FI_D(ELMNUM)
END IF

110 CONTINUE

RETURN
END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C C
C SUBROUTINE: PR_L_T C
C PURPOSE: PRINT LAYER SOLUTION C
C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

SUBROUTINE PR_L_T(OUNIT, NPEMAX, NLAMAX, NELMSH, NEMMAX,
NLAYER, NEQMAX, ZlLAY, Z2LAY,
AK, FI_D, NODE, GUL, GF, T_LAY)

IMPLICIT LOGICAL (A-Z)

integer OUNIT ,NPEMAX, NLAMAX, NELMSH,NEMMAX,
. NEQMAX, NLAYER, I , J, K1, K2,NODE

double precision ZlLAY, Z2LAY, AK, FI_D,
GUL, GF, T_B, T_T, FI_N, FI_T, Z_CENT, T_LAY

DIMENSION NLAYER(NEMMAX),NODE(NEMMAX,NPEMAX),
ZlLAY(nemmax,NLAMAX), Z2LAY(nemmax,NLAMAX),
AK(nemmax,NLAMAX), FI_D(nemmax), GUL(nemmax,2,NLAMAX),
GF(NEQMAX), T_LAY(nemmax,NLAMAX)

DO 10 I = 1, NELMSH
DO 20 J = 1, 2

IF (J .EQ. 1) THEN
T_B = GF(NODE(I,1))
T_T = GF(NODE(I,4))
FI_T = (ZZLAY(I,1)- ZlLAY(I,1)) / FI_D(I)
GUL(I,1,1) = T_B * (1 — FI_T) + T_T * FI_T
FI_T = 0.0000

DO 30 K1 = 2, NLAYER(I)

FI_N = ZZLAY(I,K1) - ZlLAY(I,1)
DO 40 K2 - 2, K1
FI_N = FI_N +(ZZLAY(I,K1) - ZZLAY(I,K2-1))*AK(I,K2-1)

40 CONTINUE

FI_T = FI_N / FI_D(I)

GUL(I,1,K1) = T_B * (1 — FI_T) + T_T * FI_T
30 CONTINUE

ELSE

T_B = GF(NODE(I,2))

T_T = GF(NODE(I,3))

FI_T = (22LAY(I,1)- ZlLAY(I,1))/FI_D(I)

179

GUL(I,2,1) = T_B * (1 — FI_T) + T_T * FI_T
FI_T a 0.0000
DO 50 K1 = 2, NLAYER(I)
FI_N = ZZLAY(I,K1) - ZlLAY(I,1)
DO 60 K2 = 2, K1
FI_N = FI_N +(ZZLAY(I,K1) — ZZLAY(I,K2-1))*AK(I,K2-l)
6O CONTINUE
FI_T = FI_N / FI_D(I)
GUL(I,2,K1) = T_B * (1 — FI_T) + T_T * FI_T
50 CONTINUE

END IF
20 CONTINUE
10 CONTINUE

WRITE(OUNIT,1000)
DO 100 I=1, NELMSH

WRITE(OUNIT,1010)
WRITE(OUNIT,1020) I, NLAYER(I)
WRITE(OUNIT,1030)
DO 90 J=1, NLAYER(I)
Z_CENT = ZlLAY(I,J) + (Z2LAY(I,J) - ZlLAY(I,J))/2
WRITE(OUNIT,1040) J, z2LAY(I,J), GUL(I,1,J), GUL(I,2,J),
Z_CENT, T_LAY(I,J)

90 CONTINUE
100 CONTINUE

C --- FORMAT STATEMENTS.

1000 FORMAT('
1010 FORMAT('
1020 FORMAT('
1030 FORMAT('

,///3X,'SOLUTION AT EACH LAYER: ')
,//6X,'ELEMENT',3X,'LAYER',/)

,9X,I3,7X,I3)

,/6X,'LAYER N.',5X,'Z LAYER',6X,'SIDE 1-4',6X,

. 'SIDE 2-3',7X,'Z(CENTER)',7X,'T(CENTER)'/)

1040 FORMAT(' ',9X,I3,4X, E10.4,4X, E11.6,4X, E11.6,3X,E10.4,4X, E14.6)

RETURN
END

180

REFERENCES

181

[1]

[2]

[3]

[4]

[5]

l6]

[7]

[8]

REFERENCES

P. C. Yang, C. H. Norris and Y. 1966. Stavsky. “Elastic Wave Propagation in
Heterogeneous Plates”, Int. J. Sol. Structures, Vol. 2:665-684.

M. Ganapathi, M. Touratier. 1997. “A Study on Thermal Postbuckling Behavior
of Laminated Composite Plates Using a Shear-Flexible Finite Element”, Finite
Element in Analysis and Design, Vol. 28:115-135.

J.N. Reddy. 1990. “A Review of Refined Theories of Laminated Composite
Plate”, Shock Vibr. Dig, Vol. 22:3-17.

T. Kant, R. K. Khare. 1994. “Finite Element Thermal Stress Analysis of
Composite Laminates Using a Higher-Order Theory”, Journal of Thermal
Stresses, 17:229-255.

K. N. Cho, A. G. Striz, C. W. Bert. 1989. “Thermal Stress Analysis of Laminate
Using High-Order Theory in Each Layer”, Journal of Thermal Stresses, Vol.
12:321-332.

R. Rolfes. 1991. “Higher Order Theory and Finite Element for Heat Conduction
in Composites”, Numerical Methods in Thermal Problems, Proceedings of the
Seventh International Conference, Vol. VII, Part 2:880-889. R.W. Lewis, J.H.
Chin, (1 GM. Homsy.

A.K. Noor and L.H. Tenek. 1992. “Steady-State Nonlinear Heat Transfer in
Multi-layered Composite Panels”, ASCE Journal of Engineering Mechanics,
118:1661-1679.

R. Rolfes. 1990. “Efficient Thermal Analysis of Anisotropic Composite Plates
Using New Finite Elements”, Developments in the Science and Technology of
Composite Materials, Proceedings of the Fourth European Conference on

182

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[13]

Composite Materials, 743-748. J. Fuller, G. Gruninger, K. Schulte, A.R. Bunsell,
and A. Massiah.

N. Mukherjee, K. Sinha. 1993. “A Finite Elements Analysis of Inplane
Thennostructural Behavior of Composite Plates”, Journal of Reinforced Plastic
and Composites, Vol. 12:1026-1042.

K. Matsumoto, M. Zako. 1995. “Study on Three Dimensional Thermal Stress
Analysis for Laminated Composite Materials”, ASME I995, PVP-Vol. 302:75-80.

J. Argyris, L. Tenek, F.Oberg. 1995. “A Multilayer Composite Triangular
Element for Steady-State Conduction/Convection/Radiation Heat Transfer in
Complex Shells”, ASCE Computer Method in Applied Mechanics and
Engineering, 120:271-301.

K. Tamma, A. Yurko. 1989. “Finite-Element Thermal Modeling/Analysis
Formulations for Layered Composite Materials”, Numerical Heat Transfer, Part B
vol. 15:73-97.

R. Rolfes, A. K. Noor, H. Sparr. 1997. “Evaluation of Transverse Thermal
Stresses in Composite Plates Based on First-Order Shear Deformation Theory”,
ASME 1997, PVP-Vol. 369:167-184.

A. K. Noor, Y. H. Kim, J. M. Peterson. 1994. “Transverse Shear Stresses and
Their Sensitivity in Multilayered Composite Plates”, AIAA Journal, Vol. 32 No.
6:1259-1269.

J .N. Reddy. 1987. “Generalization of Two-Dirnensional Theories of Laminated
Composite Plates”, Comm. Appl. Num. Meths., Vol. 3:173-180.

A. Toledano, H. Murakami. 1987. “A Composite Plates Theory for Arbitrary
Laminate Configuration”, J. App. Mech., Vol. 54:181-189.

C.Y. Lee, D. Liu. 1992. “An Interlaminar Shear Stress Continuity Theory for
Laminated Composite Analysis”, Computers & Structures, Vol. 42 No. 1:69-78.

X. Lu, D. Liu. 1992. “An Interlaminar Shear Stress Continuity Theory for Both
Thin and Thick Composite Laminate”, J. App. Mech., Vol. 59:502-509.

183

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

128]

[29]

N. Watanabe, K. Teranishi. 1995. “Thermal Stresses Analysis for Al Honeycomb
Sandwich Plates With Very Thin CFRP F aces”, AIAA, 1334-CP: 1955-1963.

A. K. Pandey, P. Dechaumphait, Allan R. Wieting. 1992. Thermal-Structural
Finite Element Analysis Using Linear Flux Formulation”, Journal of
Ihermophysic and Heat Transfer, Vol. 6 No. 2:341348.

M. DiSciuva. 1987. “An Improved Shear-Deformation Theory for Moderately
Thick Multilayered Anisotropic Shells and Plates”, Journal of Applied
Mechanics, 54:589-596.

M. DiSciuva. 1985. “Development of an Anisotropic Multilayered Shear
Deformable Rectangular Plate Element”, Computers & Structures, Vol. 21 :789-
796.

R.C. Averill. 1994. “Static and Dynamic Response of Moderately Thick
Laminated Beams with Damage”, Composites Engineering, 4:3 81-3 95.

RC. Averill, Y.C, Yip. 1996. “Development of Simple, Robust Finite Elements
Based on Refined Theories for Thick Laminated Beams”, Computers &
Structures, 59:529-546.

Y.C. Yip, R.C. Averill. 1996. “Thick Beam Theory and Finite Element Model
with Zig-Zag Sublaminate Approximations”, AIAA Journal, 34: 1627-1632.

Y.C. Yip, R.C. Averill. 1997. “A 3D Laminated Plate Finite Element Model with
Zig-Zag Sublaminate Approximations”, submitted for publication.

Y.B, Cho, R.C. Averill. 1997. “An Improved Theory and Finite Element Model
for Laminated Beams Using First Order Zig-Zag Sublaminate Approximations”,
Composite Structures, 37:281-298.

M. Cho and RR. Parrnerter. 1993. “Efficient Higher Order Composite Plate
Theory for General Lamination Configurations”, AIAA Journal, 31 : 1299-1306.

M. DiSciuva. 1993. “A General Quadrilateral Multilayered Plate Element with
Continuous Interlaminar Stresses”, Computers & Structures, Vol. 47:91-105.

184

[30]

[31]

[32]

[33]

[34]

[35]

[36]

[37]

133]

[39]

[40]

Murakami. 1986. “Laminated Composite Theory with improved In-plane
Response”, J. App. Mech., Vol. 53:661-666.

Y.B, Cho, R.C. Averill. 1999. “First Order Zig-Zag Sublaminate Plate Theory and
Finite Element Model for Laminated Composite and Sandwich Panels”, accepted
for publication.

M. Savoia and J .N. Reddy. 1995. “Three-dimensional Thermal Analysis of
Laminated Composite Plates”, Int. J. Solid Structures, Vol. 32 No. 5:593-608.

A. Pantano, R.C. Averill. 1999. “An Improved Thermal Lamination Theory and
Finite Element Model for Analysis of Heat Transfer in Composite Structures”,
Submitted to the Journal of Composite Materials in May 1999.

G. Prathap, b. P. Naganarayana. 1995 . “Consistent Thermal Stress Evaluation in
Finite Elements”, Computers & Structures, Vol. 54 No. 3:415-426.

G. Prathap, P. Naganarayana. 1997. “Accurate Thermal Stress Predictions Using
Co-continuous High- Order Shear Deformable Elements”, Computers Meth. Appl.
Mech. Eng, Vol. 144: 61 -.75

Maenghyo Cho, Heung Suk Park. 1997. “Comments on Three-dimensional
Thermal Analysis of Laminated Composite Plates”, Int. J. Solid Structures, Vol.
34 No. 35-36z4633-4651.

A. Pantano, R.C. Averill. 1999. “A First Order Zig-Zag Sublaminate Model for
Analysis of Thermal Stresses in Laminated Composite and Sandwich Panels”,
Submitted to the Journal of Sandwich Structures and Materials, in October 1999.

ANSYS Inc. . 1996. User’s Manual, edition 5.3.
R. M. Jones. 1975. “Mechanics of Composite Materials”, Taylor & Francis.

E. Reissner. 1986. “On a Mixed Variational Theorem and Shear Deformation
Rectangular Plate Theory”, Int. J. Num. Meth. Eng, Vol. 23:93-98.

185

[41]

[42]

[43]

[44]

[45]

[46]

[47]

[48]

[49]

[50]

[51]

A. Tessler, T. J. R. Huges. 1983. “An Improved Treatment of Transverse Shear in
Mindlin-Type F our-Node Quadrilateral Element”, Computer Method In Applied
Mechanics and Eng., Vol. 39:311-335.

G. Prathap, B. R. Somashekar. 1988. “Field and Edge Consistency Synthesis of a
4-Noded Quadrilateral Plate Bending Element”, Int. J. Num. Meth. Eng., Vol.
26:1693-1708.

J. Barlow. 1976. “Optimal Stress Locations in Finite Elements Models”, Int. J.
Num. Meth. Eng., Vol. 10:243-251.

Dan Zenkert. 1997. “Introduction to Sandwich Construction”, EMAS Publishing.

L. C. Davis, B. E. Atrz. 1995. “Thermal Conductivity of Metal-Matrix
Composites”, .1. App. Phys, Vol. 77:4954-4960.

M. M. Yovanovich. 1970. “On the Temperature Distribution and Constriction
Resistance in Layered Media”, Composite Material, Vol. 4:567-570.

P. S. Ramensh, M. H. Lean. 1991. “An indirect Boundary Element Method for
Heat Transfer in Composites Structures”, Proc. 0f the Seventh Int. Conference
held in Stanford “Numerical Method in Thermal Problem ”, Vol. 7 part. 2:890-
897.

L. Xingde. 1995. “Temperature Calculation in the Thickness Direction for
Composite of Arbitrary Layers”, Computers and Structures, Vol. 62 N. 4:763-
770.

L. Wang, C. Wang. 1993. “Analysis of Thermal Stresses in Composite Laminates
by Assumed Stress Hybrid Multilayer Element”, Journal of Thermal Stresses,
Vol. 13:103-112.

K. D. Cole. 1998. “Thermal Properties of Honeycomb Panels — Analysis of
Transient Experiments with One-Dimensional Models”, Procedings of the ASME
Heat Transfer Division, Vol. 5:83-93.

V. Vihak, M. Yuzvyak, A. Yasinskij. 1998. “The Solution of the Plane
Thermoelasticity Problem for a Rectangular Domain”, Journal of Thermal
Stresses, Vol. 21 :545-561.

186

[52]

[53]

[54]

[55]

[56]

[57]

[58]

[59]

N. H. Huang, T. R. M. Tauchert. 1991. “Thermal Stresses in Doubly-Curved
Cross-ply Laminates”, Journal of Solid Structures, Vol. 29 No. 8991-1000.

N. L. Breivik, M. W. Heyer. 1998. “Buckling and Postbuckling Behavior of
Curved Composite Panels Due and Mechanical Loading”, J. of Reinforced
Plastics and Composites, Vol. 17 No. 14:1292-1306.

Z. Guiping, T. Limin. 1993. “A Semi-Analytical Solution for Thermal Stress
Analysis of Laminated Composite Plate in the Hamiltonian System”, Computers
and Structures, Vol. 55 N. 1:113-118.

T. Hause, L. Librescu, T. F. Johnson. 1998. “Thermomechanical Load-Carring
Capacity of Sandwich Flat Panels”, Journal of Thermal Stresses, Vol. 21:627-
653.

J. A. Eberhardt, J. E. Locke. 1994. “Thermal-Stress Analysis Using of
Anisotropic Plates with Temperature Dependent Material Properties and Free
Edge Boundary Conditions”, AIAA, 1376-CPz552-558.

T. E. Tay. 1992. “Finite Element Analysis of Thennoelastic Coupling in
Composite”, Computers and Structures, Vol. 43 No. 1:107-112.

J.N. Reddy. 1993. “An Introduction to the Finite Element Method”, McGraw-Hill.

J. G. Crose, R. L.Holrnan, N. J. Pagano. 1987. “Validation of Advanced
Composite Thermal Stress Analysis Method”, Transaction of the ASME, Vol.
109:40-46.

187