5.1 Introduction

Finite element analysis (FEA) is an approximate numerical method that solves for the distribution of forces, displacements, or strains in a system. To complete an FEA, it is necessary to build an analytical model, which idealizes the variables of the structure under study: geometry, material properties, and boundary conditions, as well as existing or possible loadings. The idealization is based on judgment of the analyst about how the system behaves and what features are important to represent while solving the structure.

A full knowledge of the theoretical foundations behind computer programs is of primary importance to develop an accurate mathematical model. However, unlike EMA, FEA is a well-known technique used widely in engineering, and it is also presented extensively in various references [Cook et al. 2002, Bathe 1982]. Thus, an overall introduction to the theory of FEA is excluded from the scope of this thesis. Instead, the treatment of the subject is confined to the analysis of masonry domes in particular.

In this chapter, we will discuss the application of the finite element method to the solution of the dynamic parameters and the static forces in the system. The basic procedure is initiated by the development of an FE model. The next phase involves finetuning the FE model by comparisons of analytical and experimental frequencies and mode shapes. The procedure is finalized with the determination of the static state of stress, strain and displacements.

5.2 FE Model Development Applied to Domes

It is prudent to prove the assumptions in order to increase confidence in the success and accuracy of the FEA. In this study, several assumptions are accepted while adopting the FEA; however, the later analysis will show the validity of some of these assumptions.

• The deflections are small so that geometric nonlinearity does not occur.
• The materials present linearly elastic behavior.
• The material is isotropic.
• The material is homogeneous.
• The structure is composed of axisymmetric segments, and minor construction imperfections have little influence.
• The thickness of the shell is uniform.
• Basic material properties of the system (inertia and elastic forces) are independent from each other.

Modeling a real system may be problematical when all the interactions of parameters are taken into account. By the assumptions underlying most commercial FE software, all the inertia, elastic and dissipative forces are assumed independent of each other. However, the separation of the basic parameters into discrete items followed by superposition in the mathematical model is generally considered to give satisfactory results (Allemang 1993).

The commercially available FE program ANSYS 9.0 is employed throughout the analysis. In the following sections, the phases of analytical model creation are briefly discussed.

5.2.1 Geometric Model Creation

Creation of the geometric model is straightforward as long as the correct geometric representation of the actual physical structure is provided. The physical dimensions of the structural model are established either by on-site measurements or by available drawings. In cases where the domes are configured in sections of a sphere, they can satisfactorily be modeled with a very limited set of measurements of span, rise, web thickness, and surcharge height.

If a spherical dome is bounded by a square plan by cutting out four orthogonal vertical planes equidistant from the center, the edges are curved in a circular arc. The spandrel portions are called pendentives, and the dome is called a pendentive-dome [Figure 5-1]. A dome and its pendentives can simply be modeled by creating a spherical surface in the computer program. For pendentive domes, the span, radius of curvature, web thickness and plan dimensions are sufficient to build the model. Ribbed domes require additional surveying for the dimensions of the rib cross sections.

The volume formed between the adjacent elements (e.g. walls, adjacent pendentives) and the pendentives, from supports up to haunches, is usually filled with a solid mass, called surcharge volume, in order to increase the stability. The surcharge material and dimensions also need to be surveyed.

For cases when the dome webbing thickness is not definite in drawings or by visual inspection on site, some auxiliary non-destructive techniques, such as impact-echo, need to be employed.

In FE modeling, the first decision to be made is the level of simplification in depiction of geometric characteristics. As Cook et al. (2002) explained, “The modeler should seek to exclude superfluous details but include all essential features, so that analysis of the model is not unnecessarily complicated, yet provides results that describe the actual problem with sufficient accuracy.” Often masonry buildings are highly decorated, and they involve complex molding profiles, which are practically impossible to model accurately in an FE model. Attempts to represent every profile detail would increase computational cost and would not provide a justifiable increase in the accuracy [Figure 5-2].

Although a higher level of detail may be desirable in some instances, the analysis of load-bearing masonry structures does not require an extraordinary level of detail in the determination of the physical dimensions of the structure. Simplifying the complicated rib cross sections into rectangular forms of the same mechanical properties (area, moment of inertia) is a widely recommended practice.

Some specific considerations need to be observed while building an FE model for dynamic analysis. Although domical structures have an axisymmetrical arrangement, which enables the analyst to apply static analysis on one repetitive module, the use of symmetry in dynamic analysis excludes the unsymmetrical modes. Thus, particular care must be exercised while utilizing symmetry conditions; it is preferable to model the whole dome with its supporting elements for dynamic analysis. Another possible option is to use the feature of modal cyclic symmetry, which is available in ANSYS, for cyclically symmetric systems. This feature allows the review the mode shapes of a structure by modeling just a sector of it (ANSYS, 2003).

Additionally, unless the entire section responding to the force input is included in the FE model, the dynamic behavior of the structure cannot be estimated accurately. It is necessary to determine which components of the structure are involved in the dynamic response. This task is often challenging for masonry structures, and engineering intuition may not be sufficient. If possible, employing an experimental technique for this identification is expedient.

Additional mass, which is normally represented as an external loading in static analysis, is an important part of the system dynamics. Therefore, at the model creation stage, existing nonstructural mass must be denoted as rigid mass particles and be made a part of the structure for a correct dynamic representation.

5.2.2 Element Type Selection and Meshing

The discretized geometric model is composed of nodal points and connecting elements through which interelement compatibility is maintained. The field variation of state variable within the element is approximated from the nodal data: the term “finite element” refers to each individual segment over which the interpolation is applied. The designation of nodal degrees of freedom (DOF) limits the outcome and governs the spatial variation of the variable over the finite element.

Today, FE software provides numerous element types and options. Each element type has its own advantages and drawbacks; selection mainly depends on the DOF they provide. Shell elements are particularly well-suited to the modeling of vault webbing, as they are economical compared to solid elements and include all of the load effects that characterize vault webbing, including the interactions between membrane forces and bending moments.

Shell elements are developed to model singly curved or doubly curved surfaces of varying thickness. In general, they are available in both quadrilateral (four node) and triangular (three node) form elements. Shell elements have six DOF at each node: translations in the nodal x, y, and z directions and rotations about the nodal x, y, and zaxes. They have bending, membrane and shear stiffnesses. In this study, two shell elements offered by ANSYS are utilized: SHELL63 and SHELL93.

Three nodes (four nodes for a quadrilateral) with a thickness value for each node and orthotropic material properties are required to define a SHELL63 triangular-shaped element [Figure 5-3]. The thickness as well as the gradient of the variable is assumed to vary linearly throughout the element.

SHELL63, which provides limited warping ability, has both bending and membrane capabilities with a user defined “bending stiffness to membrane stiffness” ratio option. It allows application of in-plane and normal loads. SHELL63 also provides a reduced stress stiffness matrix that improves mode shape calculations (ANSYS 8.0, 2000).

SHELL93 has better ability to fit curved surfaces by the addition of midside nodes. Describing six nodes and the three thickness values for corner nodes is necessary along with orthotropic material properties to form the triangular-shaped element [Figure 5-4]. For quadrilateral forms, the definitions of eight nodes and four corner thickness values are required. The variable is interpolated in quadratic fashion within the element (ANSYS 8.0, 2000).

For homogenous materials, the shell is usually represented by the middle surface, also called reference surface, which is located in the middle of the thickness. In ANSYS, the user has a choice of where the reference surface is located within the thickness of the shell, which can be very convenient for modeling the contact of the dome web with the rest of the structure.

Once the element type selection is completed, the solid model is then subdivided into selected finite elements through a computerized systematic procedure known as meshing. The topology and size of the individual elements are influential on the analysis results, and they are controlled by different features in FE analysis. The degree of degeneracy of element shapes is controlled by aspect ratios (the ratios of the key dimensions of elements), the corner angle deviations, or by the parallelism of opposite edges. Although the elements with deformed shapes might be necessary in some cases, it is important to assure certain limitations for these shape factors by means of the mesh controls available in FE programs.

The mesh size evaluation is based on the level of approximation acceptable for the investigation. In theory, infinite mesh refinements result in a true representation of the model and give exact solutions to the mathematical problem, which is hypothetically possible but actually impractical. Therefore, it is necessary to decide the level of refinement in mesh size, remembering that too a coarse mesh causes very poor analytical results; while too fine a mesh results in computational problems.

Since the field solution is based on nodal displacement, and strains and stresses are calculated from this outcome, less accurate results for stress are found while displacements may be substantially correct. Therefore, the mesh should be further refined for regions where the accuracy of the stress distribution is important. It is therefore imperative to establish the refinement before initiating the FE calibration process.

As in all other features of FE modeling, the final determination of the acceptable level for element shapes and sizes is based on goals of the analysis and remains the analyst’s responsibility.

5.2.3 Material Properties

The real behavior of masonry is non-linear, in the sense that the stress-strain law, even for elastic compression behavior in compression, may show softening. The behavior in tension at low stress levels is inelastic due to the influence of cracking. The behavior in compression at higher stress levels is also inelastic due to irreversible softening effects. Due to the orientation of the mortar joints, the material is anisotropic and inhomogeneous.

Because it is very complex and computationally expensive to model all the material characteristics of masonry, this study proceeds on the assumption of linearly elastic, homogeneous, isotropic material properties for masonry and tile assembly. Although this is a crude approximation, the validity of linearly elastic behavior assumption is checked through experimental studies. Through material sample tests, the isotropic characteristics of the materials are investigated. Then, the obtained values are homogenized in order to achieve a reasonable isotropy-based approximation. As long as the structures under study are free from severe cracking, substantial progress in structurala assessment can be accomplished through FE analysis based on the simple form of the constitutive law.

Another challenge while defining the material properties of existing masonry structures is the fact that the documentation of the construction phases and materials are overlooked or no longer preserved. For a study on Guastavino vaulting, the situation is even more difficult since the available literature on the material properties of Guastavino vaults is very slim.

Saliklis (2003), perhaps the only scholar who studied the material characteristics of Guastavino tile, has conducted a survey on numerous small tile samples and found that the terracotta tile has significant anisotropic behavior:

The most notable result was consistent orthotropic properties. On average, the dynamic elastic modulus E_x was 21520 MPa in the X-direction, but E_y was only 12000 MPa in the Y-direction. For all tiles, E_x exceeded E_y by a factor of about 1.8. Examining the statistical variation in both directions, it is clear that the orthotropic properties are significant throughout all of the tested tiles. However, despite the presence of orthotropic properties, it is necessary to use an average value for E of about 16548 MPa in structural models, because it is assumed that
the tiles were randomly oriented in the structure.

In his further studies on the FE modeling of the Guastavino tile arches, he promotes the use of isotropic material properties for the tile. Based on this assumption, he investigates the effects of the variation in mortar modulus of elasticity, from 0.69 GPa to 16.5 GPa. He concludes that the traditional mortar used in Guastavino vaulting does not have a substantial effect on the performance of the Guastavino system. His paper suggests 0.1 for the Poisson’s ratio of the tiles.

Despite the existence of such a study on Young’s modulus and Possion’s ratio for the tile itself, quantifiable material data regarding the mortar or combinations of tile and mortar are not found in literature prior to the present study. In cases where the mortar and tile properties are known, a simple formula can be used to homogenize the two materials. However, no specifications regarding the material properties exist on mortar used in the Guastavino system. In addition to this, the increase in mortar strength due to aging is a source of additional inaccuracy.

A refined study on the density of either material has not been found in literature. Moreover, the extent to which generalization is possible for the tile and mortar of the Guastavino system is doubtful–since the process underwent a continuous experimentation within the active life of the Guastavino Fireproof Construction Company.

Establishing the material properties–Young’s modulus, Poisson’s ratio and density–for Guastavino’s assembly based the construction system is accomplished in this study by laboratory experimentation on the tile and mortar samples obtained on site. This experimentation, not intending to fill the void of literature on material properties of Guastavino vaulting, provides a modest knowledge to establish the usable range of elastic properties. This experimentation is presented in Section 5.3.

5.2.4 Boundary Conditions and Loads

The replication of existing boundary conditions is particularly challenging in solid modeling of masonry structures since the joint restraints are dependent on the configuration of members as well as the physical properties of masonry. None of the theoretically fixed or free boundary conditions truly exist in a real structure. At all times, their use introduces approximation, and the choice of one or the other  depends on the analyst’s understanding and insight as well as the technology and capacity of the FE software.

Because boundary conditions are applicable to the degrees of freedom at the nodes, rather than to the elements, the prescription of boundary conditions is related to the selection of element types. In most commercially available FE software, physical constraints are invoked by zero displacement and/or rotations at the user-defined nodes. For situations of partial restraints, an elastic foundation, which is usually simulated by a series of springs, can be utilized.

For the two case studies in this thesis, the verification of selected boundary conditions is done by comparison of dynamic parameter estimates with those of experimental results and thin elastic shell theory. A further verification is done by the visual comparisons of tension zones in the FEA results with the existing cracks on dome shells for the CCB.

The arrangement of loading conditions depends on the analysis type. In modal analysis, the results are independent from force input, and the only possible loading condition is a zero displacement constraint that defines boundary conditions. However, transient and harmonic analysis options require a general user-defined time-dependent load in addition to the provided boundary conditions. The details regarding the characteristics of the analysis options and underlying reasons for their preferences is described in following sections.

5.2.5 Analysis Type

ANSYS software provides a variety of analysis options, which have their own features. The analysis methods adopted in this study are listed below and provided with brief technical descriptions.

5.2.5.1 Dynamic Analysis

Dynamic analysis examines how the external loads are balanced with inertial, dissipative and elastic resistance of the system within the discretized time domain. It is important to note that in addition to discretization in space, dynamic analysis also requires the discretization in time. Two different alternatives of dynamic analysis are adopted in this study and briefly mentioned below.

5.2.5.2 Modal Analysis

Modal analysis is a useful feature in FE programs enabling the comparisons of analytical estimates with experimental data. FEA associates the static stiffness matrix with the mass matrix and solves for eigenvalues and eigenvectors, which represent the natural frequencies and associated mode shapes.

The basic equation solved in a typical undamped modal analysis is the classical eigenvalue problem [Equation 5-1] (ANSYS 8.0, 2000):

(Equation 5-1)

In ANSYS, modal analysis proceeds on the linearly elastic material assumption, and thus the influence of nonlinearities on the modal parameters is assumed insignificant prior to the analysis. FEA uses consistent matrices by default. There are several different modal-extraction methods available. The Block Lanczos method uses a sparse matrix solver and is preferred when the structure under consideration is large (40,000+ DOF) and symmetric. The program solves for a user-defined number of natural frequencies and mode shapes within the user-defined frequency range. The subspace method is preferred for large-scale symmetric eigenvalue problems, while the PowerDynamics method is especially preferred for obtaining the first few modal parameters of very large models (100,000+ DOF). The reduced method works with master DOF and then expands the computed modes to a full set. The accuracy depends on the master DOF selection. This method is recommended for models of less then 10,000 DOF. Certain mode-extraction techniques, also available in ANSYS, allow the analyst to introduce proportional damping, and they yield complex eigenvalues and eigenvectors (ANSYS 8.0, 2000).

Assuming the existence of reasonable meshing, doubling the number of elements reduces the frequency errors in a beam by a factor of 1/16. Therefore, obtaining a proper mesh has high importance for modal analysis results in FEA. Upon providing a correct representation of the system geometry and material characteristics with an appropriate mesh, the natural frequency estimations are upper bounds for the exact frequencies of the mathematical model (Cook et al. 1974).

5.2.5.3 Transient Analysis

Transient analysis is commonly employed to determine the time-varying response of the structure in terms of displacements, strains, stresses, or reaction forces under a user-defined loading. Essentially, the software performs a series of static analyses for each time step while considering the dissipative and inertial forces. The time increment as well as the loading conditions are also user-defined.

Transient analysis is applied to the preliminary FE model while deciding the excitation level and excitation locations prior to the field experiment. The possible time response and FRF estimates are obtained and used for the selection of exciter size and capacity.

5.2.5.4 Static Analysis

The main difference between static and dynamic analysis is that static analysis constrains all the forces to be applied slowly so that inertia or dissipative resistance are not activated. Since architecture-based systems are mostly subjected to static loads, obtaining static analysis results in order to evaluate the static state behavior is one of the ultimate goals of this study. Static analysis under gravity loading is applied to the validated FE models for both structural systems by which load distribution within the doubly curved domain and static-state behavior is studied.

5.2.6 Calibration and Validation of a Model

FE validation is essentially assessing the accuracy and reliability of the mathematical model. In this study, the reliability of the FE model is defined in relation to the experimental data and the closed form solutions of elastic thin shell theory. Assuming that the experimental results are reflections of the actual behavior of the system and remembering that many scholars have already verified the solution delivered by thin shell theory, the calibration of the FE model with respect to the two established solutions presents a tool that can provide accurate estimates for static state behavior.

It is important for the reader to remember that when testing on existing masonry structures the level of desired accuracy is not as high as the tests conducted under fully controlled conditions in labarotories. The assumptions established on the material properties of masonry prior to the study inhibit a certain level of approximation.

There are many computerized model-updating programs available, which iteratively calibrate the FE model by modifying the mass and stiffness matrices. These programs yield an FE model where dynamic parameter estimates correlate with the experimental data. In this study, however, a certain level of compromise is involved in the desired accuracy. Since the experimentally obtained results are fairly accurate, the modification of inertia and elastic resistance of the system with respect to the experimental values is not applicable.

Although the correlation of the natural frequencies obtained in this study is not as accurate as the results that might be delivered by advanced FE updating software, this study only seeks for an overall correlation in order to reach the most realistic boundary conditions and material properties. Therefore, the results are satisfactory for our purposes.

Information is gained from the experimentally acquired mode shapes about the admissible boundary conditions, and an iterative calibration process is initiated. Within each step, boundary conditions are tuned, and the mathematical model is solved for natural frequencies and mode shapes until an overall agreement between experimentally acquired modal parameters and FE estimates is achieved.

There are always problems in the correlation between FEA and EMA. The EMA data, having relatively fewer nodal points than the FEA, delivers a limited number of mode shapes. The mode shapes estimated in the FEA results that are not acquired in the experiments due to the test limitations are omitted during the validation process.

5.3 Material Sample Tests

A tile sample of 15.8 cm by 20 cm along with a 10 cm by 14 cm mortar lump is obtained from the site while executing the dynamic tests [Figure 5-5]. The tiles are cut into a square of 15.0 cm by 15.0 cm in order to eliminate the influence of differing size in perpendicular directions. The edges are ground to ensure parallel alignment of two loading edges and perpendicular angles between surfaces. Plaster capping is provided to prevent flaking of the tile edges. Since the mortar specimen has an irregular shape, the two sides are ground and also capped with plaster in order to obtain an even distribution of compressive force.

An epoxy-based compound is applied to both the tile and mortar surfaces to obtain a smooth and even face before mounting the gauges. Three electrical resistance strain gauges along each perpendicular direction, six in total, are mounted on the tile in order to average the compressive strain readings and to obtain Poisson’s ratio. Two additional strain gauges are mounted in the center on the other surface of the tile. On the mortar specimen, two strain gauges are mounted on each side–four in total.

A hydraulic testing machine is utilized to apply two-point uniform compression stress on the specimens in two orthogonal directions separately [Figure 5-6]. The tests on the tile give a stress-strain curve of linearly elastic behavior whose slope reveals the Young’s modulus of the tile [Figure 5-7]. Unlike Saliklis’s (2003) findings, similar properties for the two axes are observed. A Young’s modulus of 13.2 GPa along the longitudinal direction and 15.4 GPa along the transverse direction are found.

In two orthogonal directions, Poisson’s ratio does not exhibit a significant variation. Along the longitudinal axis, v is found to be approximately 0.182 while the value on the transverse axis is 0.194. Based on the assumption that the tiles are oriented equally in both directions, these two v values are linearly averaged to obtain the Poisson’s ratio of tile.

A similar compression test is completed on mortar sample. Since the form of the sample is amorphous, the mortar is only loaded in a longitudinal direction. Therefore, it is assumed that the values obtained through these tests are indicative for both directions. Figure 5-8 presents the experimentally obtained stress-strain graph of the mortar where it can be seen that there is a certain deviation from linear behavior. With this experiment, the approximate Young’s Modulus of 2.97 GPa and the Poisson’s ratio of 0.32 is obtained for mortar. According to the values published in McNary and Abrams for elastic constants in mortar, this is approximately twice the value of E expected for a type O mortar, but significantly lower than the expected E for a type M mortar (McNary and Abrams 1985).

The density of the tiles is obtained from a simple computation of weight and volume of the tile and mortar samples. The existing literature on the Guastavino construction system claims an equal contribution of the tile and mortar to the volume of the whole assembly. Therefore, the densities of tile and mortar are linearly averaged. The findings of this test are presented in Table 5-1.

Table 5-1: The experimentally obtained material properties for tile.

5.4 Applications to Structures in Study

The FE models for two Guastavino domes are developed by using ANSYS v. 8.0. The analytical model creation process, execution of the analysis, and inferences made from the obtained results are discussed in the subsequent paragraphs.

5.4.1 State Education Building (SEB)

5.4.1.1 Geometric Model Creation

The structural system of study in the State Education Building is arranged with twelve side-by side domes filling a three by four grid. Each dome is configured from a segment of a sphere over a square plan form supported by pendentives [Figure 5-9]. For the SEB, the architectural drawings of the masonry assembly are available to determine the physical dimensions of the dome. Survey data obtained during the site visit are also used to verify the geometric accuracy of the model. Ignoring construction imperfections, an exact spherical segmental dome is modeled in ANSYS [Figure 5-10].

There are, however no specifications of the webbing thickness on the available drawings. Further studies to experimentally determine the shell thickness reveal a thickness of 13.5 cm, and this value is approximated as a uniform and constant quantity in the FE solid model. The discussion on these test methods and results is presented previously in section 3.5.1.4.

Initial observations inside the attic space reveal that the unreinforced masonry domes are not load bearing. The steel trusses, which touch the dome shell periphery every 45° [Figure 5-9], carry the upper floor loads to the columns. At these contact locations, small tile blocks of 40 cm by 150 cm plan dimension connect the steel frame to the dome shell. These buttresses introduce restraint to the dome structures.

Since the arches and ribs are slender, the interaction between the side-by-side domes under dynamic loading is not clear [Figure 5-11]. Moreover the complex composition of the different structural members -steel trusses, tile domes, tile buttresses, and rubble-stone surcharge–presents challenges for the assessment of the dynamic characteristics of the system [Figure 5-9]. In Chapter 3, the experimental techniques are described in order to identify the integration of structural members with the center tile dome, which is selected to be the focus of the study. As illustrated in Figure 3-19 in Chapter 3, the steel truss members do not exhibit a noticeable dynamic contribution in the frequency range of interest. Adjacent domes are similarly found to be structurally independent, based on discussion in Chapter 3.

Additionally, it is considered that the volume and material of the surcharge have a limited effect on the dome behavior, and thus the surcharge volume is excluded from the solid model. In this formulation, the analysis of the whole vault system is reduced to a single dome model, which is determined adequate for the FE calibration and validation procedure. Consequently, the influence of the adjacent structural members–surcharge, steel frame, buttresses, adjacent domes–is represented by boundary condition definitions.

5.4.1.2 Element Type and Selection

The Guastavino tile domes studied in this thesis fall into the definition of “thin shell theory.” Thus, the two previously discussed shell elements, SHELL63 and SHELL93, are utilized separately in ANSYS models.

5.4.1.3 Material Properties

When attempting to analyze a tile masonry structure, the fundamental material stiffness property is an effective modulus of elasticity representing the combined effect of tile masonry units, mortar, joints and voids. The tile units were oriented in both directions in the construction of the SEB [Figure 5-12]. As noted by Saliklis (2003), a homogenization study is necessary to find a Young’s modulus value for the combined effect of the mortar and tile units that compensates anisotropic characteristics in two directions.

The tile units have different dimensions in two directions (2 cm×15 cm×30 cm) while the typical mortar joint is 2 cm in all directions. Since the present study proceeds with the isotropic material assumption, the composite material needs to be homogenized in three directions [Figure 5-13].

Discounting the slight anisotropic behavior of the laboratory results, the derived values of Young’s modulus in x and y directions are first averaged to obtain an isotropybased property. This approach assumes that the orientation of the tiles in the x and y directions have equal percentage throughout the dome shell. Subsequently, the modulus of elasticity value for the tile is homogenized with the elastic properties of the mortar based on the constant stress approximation [Equation 5-2].

For the entry of Lt and Lm in the formulae, the ratios of mortar and tile length for these two directions are averaged.  is found to be 0.08 while  is 0.92 based on the mortar joint and tile dimensions.

(Equation 5-2)

Equation 5-3 results in the homogenized Young’s modulus value of 11×10⁹ Pa along the x-y plane. The homogenization of mechanical properties of tile and mortar is completed with utilizing the rule of mixtures along the z direction [Equation 5-5]. According to the findings of the impact-echo testing, the dome is composed of four courses of tiles and hence three mortar layers. The Equation 5-4 results in 7.56×10⁹ Pa isotropic E value for the tile and mortar assembly.

(Equation 5-5)

Recalling the equal contribution of mortar and tile to the entire volume of the assembly, the Poisson’s ratio of tile and mortar is linearly averaged based on Equation 5- 6 [Table 5-2].

Table 5-3: The material properties of tile and mortar assembly input for the preliminary FE model.

These preliminary material properties of the assembly are then fine-tuned based  on the experimentally obtained natural frequencies, which are assumed to represent the actual behavior of the structure. Since only one type of material is modeled in this analysis, the natural frequencies from the FE model are dependent on the material properties. The material properties obtained through the fine-tuning process are presented in Table 5-4.

Table 5-4: The material properties of tile and mortar assembly derived from FE model calibration process.

5.4.1.4 Boundary Conditions

In order to reach accurate results, the boundary conditions are of particular interest. Because boundary conditions are applicable to the degrees of freedom at the nodal points rather than to the elements, at the points where a constraint is desired a keypoint must be defined during the geometric model creation.

Since the physical configuration of masonry assemblies places further challenges on the identification of boundary conditions, it is crucial to adjust them through a calibration procedure. An extensive experimental data set is available for the SEB, which gives information on the possible boundary conditions to work with during the calibration process. It is useful to start with physically reasonable boundary conditions in the initial FE model.

The peripheral edges of the dome are restrained from a horizontal movement perpendicular to the plane of symmetry. At the base line where the dome webbing meets the supports, translation in all three axes is restricted. The tile buttresses between the adjacent domes provide vertical restraint while the diagonal buttresses allow vertical movements. The different behavior of the buttresses is determined via the experimental studies and discussed in Chapter 3. The calibration of boundary conditions reveals that the steel I girders, contacting the dome surface along the pendentive edges at four corners, restrict rotation [Figure 5-14]. The boundary conditions obtained through the fine-tuning process are presented in Figure 5-15.

5.4.1.5 FE Model Validation

The iterative calibration process is completed when an overall agreement between experimental and analytical data is reached. The mode shapes used for comparison with the experimental data are listed in Figure 5-16.

Since the dome apices are not accessible, an impact that purely excites the axisymmetric modes is not achieved. The natural frequencies of axisymmetric mode shapes, desired by the thin shell theory, exhibit a weak peak in the FRF’s. However, it is so dominated by the adjacent bending and torsional peaks that the identification of the mode shape is not possible. For this particular building, the validation of the axisymmetric modes will only be achieved through the estimates of the thin elastic shell theory [Table 5-5].

Table 5-5: The natural frequencies obtained through three methodologies, SEB, NY.

5.4.1.6 Sensitivity Analysis

The sensitivity analysis is initiated by investigating the influence of changes in material properties. Since the solid model is reduced to the tile dome, only one type of material is modeled in the FE model. Therefore, the natural frequencies obtained through FE analysis are completely dependent on the E/ρ ratio entered into the program. As long as this ratio is kept constant, the natural frequencies and the mode shapes remain identical. In other words, the fine-tuning of material properties in relation to experimental data only reveals information on the E/ρ ratio, not the values of each. Therefore, a further study to determine the individual values for E and ρ is necessary before initiating the static state analysis.

For the SEB, material samples form the site are available. The values obtained through the homogenization of laboratory test results are then calibrated based on the dynamic testing data. The results of the sensitivity analysis of 1% and 10% are presented in Table 5-6. It is important to note the linear variation of the natural frequencies due to a change in the material properties.

Table 5-6: The change in natural frequency due to 1% and 10% increase in density.

In order to determine the influence of the mesh size, type and refinement on the computed results, a further sensitivity analysis is completed. A model meshed with the SHELL63 element results in 2.5% higher natural frequency for the first mode, while a closer correlation is observed for higher order modes. Although the bending and torsion based mode shapes exhibit reasonable correlation, the axisymmetric modes are not estimated at all when the FE model is meshed with SHELL63 elements. This is probably due to the insufficient warping capacity of these elements. For further studies, SHELL93 elements are recommended.

Both quadrilateral and triangular elements are used in meshing as illustrated in Figure 5-17. The solutions of both for natural frequencies show a reasonable agreement within a range of 1%.

Table 5-7: The variation of natural frequencies (Hz) due to automatic mesh refinement.

The influence of the mesh size refinement is also investigated in three steps by the size refinement feature in the FE program [Figure 5-18]. Through a sensitivity study, it is found that the FE model converges when a mesh size of approximately 20 cm is provided for Guastavino tile masonry structures at the scale of 5-10 m radius [Table 5-7]. A further mesh refinement is not necessary, as for the first mode the difference in natural frequencies from the refined mesh (5%) is the same as the uncertainty in material properties. Average of values of E from 7.4 to 7.8 Gpa was found, giving a variation of 5%.

5.4.1.7 Making Inferences from the FE Model

A careful visual site survey on the tile domes at the SEB revealed that the domes do not exhibit any apparent cracks. Hence, the maximum tensile or compressive stress can satisfactorily be computed via linearly elastic FE analysis.

The element solutions for principal stress distributions are calculated. The most critical tensile stresses can be found from an evaluation of the first principal stress plot. The first principal stress contour plot reveals tensile stress concentrations along the edges [Figure 5-19]. The magnitudes of these stresses are around 0.17 MPa, which is approximately 8.7% of the allowable tensile stress coefficient presented by Guastavino (Guastavino 1892). The compressive stresses in this plot are as low as 5.4 x10^-5 MPa.

The second and third principal stresses present the highest magnitude compressive stress concentrations around the tile buttresses, with a compressive stress magnitude of approximately 1.0 MPa or 7% of the compressive capacity specified by Guastavino (Guastavino 1892). While the second principal stress displays tensile stress values of approximately 0.069 MPa, the third principal stress contour plot reveals almost no significant tensile stress [Figure 5-20].

In all cases, none of the principal stresses is high enough to induce inelastic or nonlinear behavior of the tile. Similarly, the stress distribution in global x, y and z-axes exhibit parallel results [Figure 5-21]. In all cases, the stress level is less than 10% of the allowable tensile stresses coefficient as stated by Guastavino (Guastavino 1892).

Subsequently, reaction forces at boundary locations–illustrated in Figure 5-22 – are obtained under gravity loading. The results show that the steel girder contacting the dome shell at the tile buttresses between adjacent domes contributes significantly to the vertical load transfer of the total dome weight of 10.6 tons. Approximately 15.7 kN is carried by the steel truss at the buttress locations between the domes, while approximately 10.9 kN is transferred through the surcharge volume.

The horizontal reaction force counterbalanced by the adjacent domes is calculated to be 20.7 kN in total along one side of the pendentive dome. On the other hand, in the x and y directions, the horizontal reaction forces where the dome web meets the surcharge is observed to be limited to 3.4 kN.

The buttress locations between the domes present very little reaction force, 0.007 kN, in the direction parallel to the symmetry plane. Although the proceeding analysis maintains that support restraint, the necessity of such a restraint is doubtful.

5.4.2 City County Building (CCB)

A brief description of this building is provided in section 1.4. The highly decorated domes have varying thicknesses ranging between three and seven tiles [Figure 5-23]. Based on the discussion on the simplification of the superfluous details in solid modeling, a constant and uniform thickness of four tile layers (13.5 cm) is utilized.

5.4.2.1 Geometric Model Creation

All the construction drawings of the building are available, and the physical model is developed according to the construction dimensions found in these documents. Although the measurements taken on site showed small differences for the side lengths of the square, the construction imperfections are neglected, and the dome is modeled as a symmetric spherical segment with a square plan form.

The three side-by-side domes are located above the vestibule entrance. Massive brick arches at four sides support these pendentive domes. During the field tests, the dynamic interaction between the adjacent domes as well as between the dome webbing and supportive arches is investigated. The dynamic test results–discussed in Chapter 3– reveal that an excitation on the dome webbing is not transferred to either the adjacent domes or the brick arches. This can be explained with the fact that the massive nature of the arches is providing a fixed end condition along the peripheries of the dome shell. This observation leads to the modeling of a single dome, and representation of the physical configuration of the adjacent members as boundary conditions.

Additionally, tile buttresses, extending between the dome webbing and the steel columns at four diagonal pendentive corners at the height of the top of the arches, support the dome as seen in Figure 5-24. The influence of these buttresses is also reproduced with boundary conditions.

5.4.2.2 Element Type and Meshing

As in the SEB, the dome shell is meshed with SHELL93 elements. Based on the sensitivity analysis completed for the SEB, a maximum mesh size of 20 cm is adopted [Figure 5-25].

5.4.2.3 Material Properties

The contract specifications of the CCB building reveals extensive information on the material quality Guastavino has been asked to provide. The contract specifications calls for “sound, hard-burned, semi porous terra cotta tiles of 7/8 “ thick by 6” wide and 12” long” (CCB Contract and Specifications 1915). The mortar is specified to be composed of one part American Portland cement to two parts sand (CCB Contract and Specifications 1915).

Although the contract descriptions for “timbrel arching” provides valuable information regarding the scope and conditions of the work Guastavino was required to do, no quantifiable information regarding the material properties is provided. Additionally, due to the lack of samples or coupons for the CCB, the validated material properties of the SEB are used for the preliminary analysis. Through a calibration process, these material properties are fine-tuned, and the final values, presented in Table 5-8, are obtained.

Table 5-8: The material properties of tile and mortar assembly, CCB, PA.

5.4.2.4 Boundary Conditions

The massive and stiff peripheral arches are considered to restrain the dome edge as a fixed end condition. The four tile buttresses in the middle of the pendentive corners are modeled as horizontal displacement restraints in both the x and y directions. Although not measured, the surcharge volume is observed to have a low rise, and its influence on the dynamic behavior is therefore excluded from the analysis.

The tile domes contact the balcony floor at their apex. Since the balcony floor has closely spaced beams (two on each side of the apex), the domes are presumed nonbearing. However, through the calibration process, the influence of the contact surface between the dome and the upper floor is found to restrain the dome webbing horizontally in both the x and y directions [Figure 5-26].

5.4.2.5 FE Model Validation

The fine-tuning of boundary conditions and material properties is completed when an overall correlation between experimental and analytical data is achieved. Since the experimental data in the CCB presents limited knowledge on the mode shapes, the calibration process is progressed using comparisons of the sign of the displacements of available data points rather than global mode shapes. As the experimental data are of lower quality compared to the SEB, a higher degree of compromise in the correlation is established. Still, an overall agreement with the mode shapes and with the natural frequencies is achieved. Figure 5-27 presents the mode shapes, which are matched to the experimental data. As explained previously, the modes, which are estimated by the FE model, but not acquired by EMA, are not taken into account in this study.

For a second time, it is observed that only the first two axisymmetric modes are observable in the pendentive domes of the CCB: the higher order axisymmetric domes are not observed [Table 5-9].

Table 5-9: The natural frequencies obtained through three methodologies, CCB, PA.

5.4.2.6 Making Inferences from the FE Model

The condition survey of the back of the tile domes revealed that approximately 2 mm cracks, in line with the pendentives and directed toward the apex, have developed [Figure 5-28].

Unlike the SEB, the intrados of vestibule domes of the CCB are exposed to the exterior fluctuation in temperature while the extrados within the building are less susceptible to temperature changes. Although the initial cause of these cracks may be the structural loads, the cyclic temperature loading on this structure certainly influences the development of the cracks.

An interesting point to note about this building is that the original construction drawings do not indicate the tile buttresses. They may have been built to prevent further crack development after completion of the construction or may simply be a revision on the design of Guastavino.

The principal stresses are obtained on completion of the static analysis on the validated FE model [Figure 5-29]. Similar to the SEB, the first principal stresses revealed compressive stress concentrations along the edges where the arches and dome web meet, however, with insignificant amplitudes. No crack development or separation is observed between these two elements of the system.

In the plot of the first principal stress, at the exact crack locations, side by side tensile and compressive stress concentrations are also observed. The tensile stresses are of approximately 2.2 MPa magnitude. Assuming that Guastavino’s tests and coefficients are credible, the tensile stress concentrations exceed the tensile capacity of tile and mortar assembly in those regions.

The intensity of the compressive stresses around the buttresses can be identified on the second principal stress plot [Figure 5-30]. The intensity of the stresses is 0.5 MPa, considerably lower than the compressive capacity coefficient of Guastavino. However being adjacent to a high tensile zone, the compressive stresses possibly contribute to the cracks observed at those regions. Similar to the results of the SEB, third principal stress is purely compressive, with a maximum value of 0.23 MPa.

The first principal stress reveals tensile stress concentrations high enough to cause the material to exhibit nonlinear behavior in addition to the nonlinearities introduced due to opening and closing of the existing cracks. This divergence from the linearly elastic behavior is also confirmed by the reciprocity and linearity checks on modal testing data– discussed in Chapter 3. However, although the experimental results are displayed nonlinearity, in this study the level of non-linearity in this system of domes is considered to have an insignificant affect on the load path analysis. The following paragraphs describe the load transfer between the elements for domes of the CCB.

The static analysis proceeds with the evaluation of the reaction forces. The results obtained from the final FE model shows that the buttresses account for a limited amount of horizontal forces, while most of the horizontal thrust is transferred to the massive arches (54.7 kN). The tile arches carry the entire vertical load (72.6 kN), and ultimately transfer it to the piers. Almost no vertical load is resisted at the base where the vault web meets the surcharge volume.

5.5 Concluding Remarks

Using the computerized tools of FEA, it becomes possible to solve for the governing equations of domed structures. The objective of this chapter was to present the steps necessary to follow while utilizing the FE analysis tools particularly for large-scale masonry structures.

The complete process of the concept is exemplified in two structures, starting from the analytical model creation and ending with the inferences from the results. Within this procedure, the supportive data on the material properties and dome web thickness are used to achieve a better representation of the actual physical structure. The calibration and validation of the FE models precede based on the dynamic parameter results obtained from EMA and the thin elastic shell theory.

The sensitivity analysis, which gives valuable information on the significance of the input entry, is completed. A 1% change in the material properties causes an approximately 0.5% change in the natural frequencies. The precision required for the material property input is dependent on the desired accuracy of the solution. The results indicate that shell elements with midside nodes give more accurate results compared to those with only corner nodes. It is seen that the choice of three or four node elements does not make a significant difference on the results. A mesh element size approximately around the size of the tiles is observed to give accurate results.

As a final step, the stress distribution and reaction forces are obtained through the static analysis under gravity loading. The results show low stresses generally for the SEB and critical tensile stress concentration around the buttresses for the CCB. The findings of the static analysis are then used in the discussion on Guastavino domes in Chapter 6.