Assessment of unreinforced concrete and masonry structures against potentially disastrous load cases such as seismic or settlement loading remains a challenge. A finite element based technique, Sequentially Linear Analysis, avoids many of the issues commonly associated with the modelling of brittle materials. One such issue is the failure of non-linear iterative finite element methods to converge on a solution when modelling brittle materials exhibiting a global ‘snap-back’ behaviour. Sequentially Linear Analysis solves such issues by linking together a series of linear-elastic analyses with a unit of material cracking applied between each stage.
While Sequentially Linear Analysis has shown promise in modelling concrete structures failing in rocking and sliding, it employs a simple isotropic uniaxial failure criterion and as a result falls short of adequately modelling anisotropic materials such as masonry and consistently overestimates the capacity of structures where shear is a predominant failure mechanism. In addition, Sequentially Linear Analysis takes no account of stress redistribution occurring between analysis stages. This thesis resolves these issues by incorporating an anisotropic biaxial failure criterion and new stress redistribution criteria into Sequentially Linear Analysis. The validity of the solutions is presented by comparison of Sequentially Linear Analysis simulations to destructive laboratory tests on masonry structures where shear is a predominant failure mechanism.
The new model is found to provide much more accurate prediction of the capacities of the structures simulated, however, would tend to overestimate the global stiffness of the structure.
Unreinforced brick and mortar type structures exist throughout the world and need to be assessed for potentially disastrous loads, e.g. seismic and settlement loads. To accomplish this, researchers have experimentally tested large masonry panels, or even full scale structures. Testing can be prohibitively expensive for the majority of structural assessments, as such these tests need to be interpreted for use with computer modelling, finite element analysis is one approach. Finite element approaches range from micromechanical modelling of individual brick and mortar elements to large scale continuum approximations to masonry; however, the softening behaviour of brittle materials still provides a challenge to computational methods. Advanced iterative methods have been developed to handle this behaviour but still fail to converge when the model exhibits a large 'snap-back' global response.
Due to these challenges, simple hand calculations provide a first-order alternative assessment approach. One may consider failure at first cracking, however this can understate the ultimate strength of the structure by a significant margin. Additionally, modelling methods which assume drystone blockwork failure mechanisms do not consider the tensile strength at all and often do not allow for fracture of bricks. Thus, despite the challenge of finite element modelling, it remains the predominant approach.
Various methods of computational modelling of masonry structures have been developed. Early work by Lourenço1 involved finite element techniques using individual brick and mortar elements, however the very large mesh sizes required for reasonably sized structures make such an approach computationally difficult outside of research use. Homogenised masonry finite elements – approximating masonry as a continuum – exhibiting a non-linear strain-softening behaviour greatly reduce the number of required elements. Lourenço further developed a homogenised anisotropic continuum model2 for analysis of masonry. Existing commercial finite element packages can reproduce this homogenised behaviour using iterative methods, however they are often slow to converge and convergence can fail due to the effects of the negative tangent stiffness associated with material softening.
Rots3 proposed an alternative method of analysing structures, Sequentially Linear Analysis, which overcomes iteration issues by performing a sequence of linear analyses and applying a unit of cracking – or 'damage' – to the model between each step. This method has been shown to be effective at modelling concrete structures, including reinforced concrete structures. Sequentially Linear Analysis has also been used to model masonry structures but still retains some significant limitations, namely:
In general, the first of these two of these issues are particularly important for modelling masonry, as they result in the overestimation of the shear capacity of masonry panels. The third issue needs to be addressed for Sequentially Linear Analysis in general.
The main aims of this thesis are in extending Sequentially Linear Analysis to address the issues listed in the previous section, these aims are:
In this report, chapter 2 reviews Sequentially Linear Analysis (SLA). The details of the implementation of extending SLA to include a biaxial failure criterion for both concrete and masonry are then presented in chapter 3. Additionally, chapter 3 discusses the details for improving the load redistribution criteria between damage increments. In chapter 4, laboratory tests on two masonry structures are reviewed: a simple masonry panel, followed by a full scale building façade. Each test is simulated using the revision of Sequentially Linear Analysis from before and after the objectives listed in this section were achieved. Finally, comparisons between the two are made with one another and with the experimental results.
This chapter contains the prerequisite information for proceeding with the achieving the objectives listed in chapter 1. Section 2.1 describes some general terminology for masonry structures. Section 2.2 reviews the latest developments in Sequentially Linear Analysis in order to lay the foundations on which this thesis is based.
Three main failure mechanisms of masonry associated with large horizontal loads are discussed in this report: rocking, sliding and shear. Rocking and sliding appear similar in that both manifest in horizontal cracks (cracking parallel to the bedding joints) whereas shear failure manifests in diagonal cracks. Rocking and sliding differ from one another in the type of rigid body motion that occurs once a large enough crack has formed for the structure to break in two. If the forces applied to the free section are large enough for it to overturn, i.e. if there is an unbalanced moment applied, rocking failure occurs. If the forces applied are large enough to overcome the friction between the two sections of the structure, sliding failure occurs. As rocking failure requires as large a moment as possible, it is typical of rocking failure cracks to form at the very base of a structure in order to maximise the lever arm of the overturning moment, whereas sliding failure can occur at any height. Shear failure occurs as a result of large shear stresses and forms diagonal cracks. Panels or piers of masonry structures with particularly low aspect ratios (height:width ratio of around 1.0 to 1.5) under large compressive forces are particularly susceptible.
As mentioned in chapter 1, Sequentially Linear Analysis uses a sequence of linear elastic analyses to model the highly non-linear response of brittle structures. For each stage in the analysis, the structure is loaded up until one point fails. At the point of failure, a unit of damage is applied and the process is repeated. Figure 1 illustrates the process of Sequentially Linear Analysis in the abstract sense.
Sequentially Linear Analysis is implemented in FORTRAN as a user defined material model for the finite element package, DIANA. The tasks of applying and removing loads, and performing the linear analysis are features included in DIANA. The steps “Scale up the load until one element fails” and “Modify the model to account for cracking at the point of failure” are performed within the user defined material model. Sequentially Linear Analysis has been refined to consider mesh objectivity4, to model reinforced concrete structures5, to consider non-proportional loading6 and to perform 3D analysis using shell elements7.
“Modify the model to account for cracking at the point of failure” – the damage step – is already well implemented and no changes are proposed in this report.
“Scale up the load until one element fails” – the scaling step – does not consider biaxial effects on the strength of the structure. The strength of structures in combined shear and compression will always be overestimated.
Sequentially Linear Analysis uses a smeared crack damage model: the material properties at the integration points within an element are reduced to represent the amount of cracking contained within the element. Many schemes have been proposed for performing this material damage and these generally fall into two categories: fixed crack models and rotating crack models.
Fixed crack models allow cracking to initiate in any direction from an integration point, but then must continue to crack in the same direction regardless of how the stresses subsequently change. Fixed crack models are physically more realistic but will often overestimate strengths as stresses will tend to reorientate to a direction in which cracking is forbidden; this is referred to as 'stress locking'. Stress locking can be alleviated by allowing more than one crack to form in different directions from the same integration point. Additionally, fixed crack models must allow for the transfer of shear stresses across the crack due to friction or aggregate interlock effects. The amount of shear transfer is controlled by a shear retention factor. The shear retention factor is described later in this section.
Rotating crack models also may initiate cracking in any direction, but will reorientate the crack according to the principal stress directions for each step of the procedure. In this case, the transfer of shear stresses across the crack becomes irrelevant (in the principal stress orientation, the shear stress will always be zero). Rotating crack models are physically unrealistic and will usually underestimate strength, however they remain popular for assessment purposes due to their conservative nature.
Sequentially Linear Analysis uses a fixed-orthogonal crack model: cracking is fixed in the direction it initiates in, and a second crack is allowed to propagate perpendicularly to the first. Empirically, this has been found to offer the best compromise solution. However, the shear stiffness needs to be adjusted to limit stress locking.
The stiffness at an integration point can be built up by considering two Young's moduli, \(E_{1}\) and \(E_{2}\), representing the stiffness in the two crack directions; two associated Poisson's ratios, \(\nu_{12}\) and \(\nu_{21}\); and a reduced shear modulus, \(G_{red}\), allowing for the transfer shear across the cracks. A simple linear elastic constitutive relationship can then be built up.
\[ \begin{bmatrix}\sigma_{11} \\ \sigma_{22} \\ \tau_{12}\end{bmatrix} = \begin{bmatrix} \displaystyle \frac{E_{1}}{1-ν_{12} ν_{21}} & \displaystyle \frac{ν_{12} E_{1}}{1-ν_{12} ν_{21}} & 0 \\ \displaystyle \frac{ν_{21} E_{2}}{1-ν_{12} ν_{21}} & \displaystyle \frac{E_2}{1-ν_{12} ν_{21}} & 0 \\ 0 & 0 & G_{red} \\ \end{bmatrix} \begin{bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \gamma_{12} \end{bmatrix} \label{eq2_1}\tag{2.1} \]
where \(\sigma_{11}\), \(\sigma_{22}\) and \(\tau_{12}\) represent the stress state in the orientation of the crack and \(\epsilon_{11}\), \(\epsilon_{22}\) and \(\gamma_{11}\) represent the strain state in the same orientation.
In order to calculate the reduced Young’s moduli, a complete stress-strain response must be selected for the material. A simple linear softening curve is shown in Figure 2a, this is a first order approximation of the typical response seen in brittle materials such as concrete and masonry. The shape of the curve is defined by a tensile strength, \(f_{t}\), a Young’s modulus, \(E\), a fracture energy, \(G_{f}\), and a crack bandwidth, \(h\).
The fracture energy, \(G_{f}\), can be considered a material property and is defined as the energy required per unit area to form a crack. The crack bandwidth, \(h\), is the characteristic size of the computational finite element (units of length) and exists to solve the problem of mesh dependence. Rots suggested empirical rules for calculating the crack bandwidth for various element shapes8.
Figure 2b shows how the stress-strain curve is discretised into a number of finite damage increments. While the shape of this ‘sawtooth’ curve is chosen to preserve its area (and therefore the fracture energy) compared with the original softening curve, the model allows the user to set the number of sawteeth in order to trade off computational time with accuracy. As damage is applied to an integration point, the model steps down from tooth to tooth. For each stage of damage, the new tensile strength becomes the peak of the tooth, and the new Young's modulus becomes the slope. The Poisson’s ratio is scaled pro-rata with the Young’s modulus. The reduced shear modulus is calculated using the material properties of the larger of the two cracks and a shear retention factor, \(\beta\). For example, if \(E_{1} < E_{2}\).
\[ G_{red} = \beta \frac{E_{1}}{2(1 + \nu_{12})} \label{eq2_2}\tag{2.2} \]
The shear retention factor, \(0.0 \leq \beta \leq 1.0\), is a measure of the degree of interlock at crack faces. In modelling concrete structures, \(\beta\) tends to depend on aggregate size. In masonry structures, it can usually be taken to equal 1.0.
In order to determine the point to apply damage, the loads must be scaled up until one point fails.
Sequentially Linear Analysis accepts two load cases: load case A, an optional fixed load, typically used to represent dead loads; and load case B, the load case to be scaled, typically used to represent the live loads, e.g. seismic or settlement loading.
Elastic analyses for load case A and load case B are run separately and the stress states for each integration point evaluated, \(\sigma_{A}\) and \(\sigma_{B}\). A load factor, \(\lambda\), is then calculated for each point such that the combined stress state, \(\sigma_{A} + \lambda \sigma_{B}\), is on the point of failure. The existing model uses a simple uniaxial failure criterion: for crack initiation this is when either principal stress exceeds a set tensile strength; for this case, DeJong et al.6 showed that \(\lambda\) is given by:
\[ \lambda = \frac{-b \pm \sqrt{b^{2} - 4ac}}{2a} \] where
\[\begin{align} a & = \sigma_{xy,B}^{2} - \sigma_{xx,B}\sigma_{yy,B} \\ b & = f_{t}(\sigma_{xx,B} + \sigma_{yy,B}) + 2 \tau_{xy,A} \tau_{xy,B} - \sigma_{xx,B}\sigma_{yy,A} \\ c & = \sigma_{xy,A}^{2} - \sigma_{xx,A} \sigma_{yy,A} + f_{t}(\sigma_{xx,A} + \sigma_{yy,A}) - f_{t}^{2} \label{eq2_3}\tag{2.3} \end{align}\]
where \(f_{t}\) is the material tensile strength and \(x\), \(y\) and \(z\) are the global co-ordinate axes.
If stresses are being considered for a partially damaged integration point with one or two existing cracks, further damage is controlled only by the stress normal to a crack. The stresses are resolved in the orientation of the primary crack and the following formulae are used:
\[ \lambda = \frac{f_{t} - \sigma_{nn,A}}{\sigma_{nn,B}} \textrm{ or } \lambda = \frac{f_{t} - \sigma_{pp,A}}{\sigma_{pp,B}}\label{eq2_4}\tag{2.4} \]
where \(n\) and \(p\) are the axes normal and parallel to the primary crack – the fixed-orthogonal crack model dictates that the secondary crack will always form normal to the primary crack.
Note that if only one of the two cracks allowed by the crack model has begun to form, the unformed crack is treated as a crack of size zero and the above formulae are still used with the full material tensile strength.
Once all values of \(\lambda\) have been calculated for every integration point, the critical value is chosen, and damage applied.
DeJong et al.6 explain in detail the procedure for selecting the critical load factor, \(\lambda\). A brief overview is given here. The calculated load factors are split into two categories, \(\lambda_{t}\) and \(\lambda_{c}\). The \(\lambda_{t}\) values cover the majority of elements: where load case A is well within the failure surface, and an increasing load case B is required to fail the element. The \(\lambda_{c}\) values arise where load case A has already exceeded the strength of the material, and load case B is required to restore it back to safety.
In most cases all the \(\lambda_{t}\) values will be greater than all (if any) the \(\lambda_{c}\) values and the minimum \(\lambda_{t}\) is selected as the critical \(\lambda\); however, when a \(\lambda_{c}\) value exceeds a \(\lambda_{t}\), no equilibrium solutions exist within or on the failure surface so the element corresponding to the maximum \(\lambda_{c}\) is selected to apply damage at an arbitrary load factor. Extra care must be taken in the event of negative load factors.
The failure criterion described in the previous chapter is now extended. In section 3.1, the failure criterion is extended to take account of isotropic biaxial effects. This section looks at unreinforced concrete, a brittle, isotropic material. Biaxial tests on concrete panels are reviewed and, from which, a numerical failure surface is derived and incorporated into Sequentially Linear Analysis. Section 3.2 investigates biaxial tests on masonry panels and uses the framework introduced in section 3.1 to incorporate an anisotropic biaxial failure surface appropriate to masonry into SLA. Sections 3.3 and 3.4 contain further notes on implementation of the biaxial framework. Section 3.5 looks at introducing the stress redistribution criteria.
Unreinforced concrete exhibits an isotropic behaviour that is otherwise similar to masonry. Kupfer et al.9 studied the in-plane behaviour of various unreinforced concrete panels and their results are shown in Figure 3a. More recent tests by Lee et al.10 offer similar results. Note that in the compression-tension region, the tensile strength tapers off with increasing compression.
Figure 3b shows an idealised version of the failure surface in terms of the principal stresses; Figure 3c shows the same idealised failure surface in terms of the axial and shear stresses in the global co-ordinate system as this representation was found to greatly simplify the maths and provide scope for later including anisotropic effects.
The three cones that make up the failure surface in Figure 3c are defined using the formula:
\[ A\sigma_{yy}^2+B\sigma_{xx}^2+C\tau_{xy}^2+D\sigma_{xx} \sigma_{yy}+E\sigma_{yy}+F\sigma_{xx}+1=0\label{eq3_1}\tag{3.1} \]
where the constants are shown in Table 1 (normalised by the tensile strength).
| Cone | A | B | C | D | E | F |
|---|---|---|---|---|---|---|
| T-T | 0.00 | 0.00 | -1.00 | 1.00 | -1.00 | -1.00 |
| C-T | -0.10 | -0.10 | -1.21 | 1.01 | -0.90 | -0.90 |
| C-C | 0.00 | 0.00 | -0.01 | 0.01 | 0.10 | 0.10 |
Table 1: Isotropic biaxial failure surface cone formula (equation 3.1) constants
The strength of this approach comes from its flexibility. A wide variety of possible failure surfaces can be represented by varying the coefficients and number of cones used. Further details on the use of this formula can be found in section 3.3. The T-T cone corresponds exactly to the tension-tension region of the failure surface, and the C-T cone covers the two compression-tension regions. The C-C (compression-compression) cone has been listed for completeness, however the current implementation of Sequentially Linear Analysis for non-proportional loading does not consider compression failure.
Using the new cone model, the formula for determining the load factor, \(\lambda\), becomes:
\[ \lambda = \frac{-b \pm \sqrt{b^{2} - 4ac}}{2a} \] where
\[\begin{align} a&=A\sigma_{yy,B}^2+B\sigma_{xx,B}^2+C\tau_{xy,B}^2+D\sigma_{xx,B} \sigma_{yy,B} \\ b&=2A\sigma_{yy,A} \sigma_{yy,B}+2B\sigma_{xx,A} \sigma_{xx,B}+2C\tau_{xy,A} \tau_{xy,B}+D(\sigma_{yy,A} \sigma_{xx,B}+\sigma_{yy,B} \sigma_{xx,A})+E\sigma_{yy,B}+F\sigma_{xx,B} \\ c&=A\sigma_{yy,A}^2+B\sigma_{xx,A}^2+C\tau_{xy,A}^2+D\sigma_{xx,A} \sigma_{yy,A}+E\sigma_{yy,A}+F\sigma_{xx,A}+1) \label{eq3_2}\tag{3.2} \end{align}\]
For crack initiation, this formula can be used as is, however for continued propagation, further damage is controlled only by stresses normal and perpendicular to the existing crack (the shear stress in this orientation has no effect).
In order to determine the new failure surface for the cracked element, a slice through the original failure surface should be taken at the angle of the newly formed crack and used as a two dimensional failure surface. Rather than take this approach directly, it is simpler to project the element stresses onto the plane that this slice forms in the global axes. The new stresses can then be used with the cone formula as above. The describing formula for the plane formed by the slice is as follows:
\[ \tan \theta = \frac{2\tau_{xy}}{\sigma_{xx}-\sigma_{yy}} \label{eq3_3}\tag{3.3} \]
where \(\theta\) is the angle of the crack with respect to the global x-axis.
The following formulae are used to transform the stresses:
\[\begin{align} \sigma_{xx}'&=(\sin^{4} \theta+cos^4 \theta)\sigma_{xx}+2\sin^{2} \theta\cos^{2} \theta\sigma_{yy}-2\sin\theta\cos\theta(\sin^{2} \theta-\cos^{2} \theta)\tau_{xy} \\ \sigma_{yy}'&=2\sin^{2} \theta\cos^{2} \theta\sigma_{xx}+(\sin^{4} \theta+cos^4 \theta)\sigma_{yy}+2\sin\theta\cos\theta(\sin^{2} \theta-\cos^{2} \theta)\tau_{xy} \\ \tau_{xy}'&=-(sin\theta\cos\theta)(\sin^{2} \theta-\cos^{2} \theta)\sigma_{xx}+(sin\theta\cos\theta)(\sin^{2} \theta-\cos^{2} \theta)\sigma_{yy}+4\sin^{2} \theta\cos^{2} \theta\tau_{xy} \label{eq3_4}\tag{3.4} \end{align}\]
These equations perform the following task:
This has the effect of locking the principal stresses in the orientation of the crack and setting them equal to the stresses normal and parallel to the crack. This method has the advantage that any future changes to the failure surface for crack initiation will not require additional changes to the propagation phase.
Dhanasekar et al.11 combined the results of a number of earlier tests on masonry panels and suggested an anisotropic failure surface using the constants in Table 2 with the cone formula (equation 3.1).
| Cone | A | B | C | D | E | F |
|---|---|---|---|---|---|---|
| 1 | 0.00 | 0.00 | -1.47 | 1.38 | -1.33 | -1.04 |
| 2 | 0.01 | -0.01 | -0.73 | 0.48 | -0.81 | -0.58 |
| 3 | 0.00 | 0.00 | -0.01 | 0.00 | 0.06 | 0.05 |
Table 2: Masonry constants (normalised by the tensile strength parallel to the bedding joints)
The three cones have been numbered, as they do not match exactly the three regions of the failure surface as with the isotropic cones. Investigation into these constants found them to be a poor fit of the data collected. Furthermore, it is rather user-unfriendly of the model to expect these 18 figures to define each new failure surface – none of the constants correspond to the results of any standard material tests.
If the general shape of the failure surface is assumed to remain roughly the same between different masonry types, simple relationships can be built up between the results of standard material tests and the cone constants required by the model. An obvious choice for the material tests on which to base these relationships would be the four uniaxial strengths: tensile strength normal and parallel to the bedding joints (\(f_{t,n}\) and \(f_{t,p}\)) and compressive strength normal and parallel to the bedding joints (\(f_{t,n}\) and \(f_{t,p}\)). The failure surface proposed by Dhanasekar et al. was redefined in terms of the four uniaxial strengths and is shown in Figure 4.
This failure surface can be constructed using two cones. The formulae used to generate the required cone constants for this failure surface are as follows:
\[\begin{align} A &= 0.0 \\ B &= 0.0 \\ C &= -\bigg(\frac{4}{9} \frac{f_{t,n}^{2}}{f_{t,p}^{2}} + 1 \bigg)^{-1} \\ D &= \frac{f_{t,p}}{f_{t,n}} \\ E &= -\frac{f_{t,p}}{f_{t,n}} \\ F & = -1.00\label{eq3_5}\tag{3.5} \end{align}\]
\[\begin{align} A &= -\frac{2f_{t,p}^{2}}{3f_{t,n}f_{c,n}} \\ B &= -\frac{f_{t,p}}{2f_{c,p}} \\ C &= -\bigg(\frac{16}{9} \frac{f_{t,n}^{2}}{f_{t,p}^{2}} + 1 \bigg)^{-1} \\ D &= \frac{f_{t,p}}{3f_{t,n}} + \frac{f_{t,p}^{2}}{f_{t,n}f_{c,n}} \\ E &= \frac{f_{t,p}}{f_{c,n}} - \frac{2f_{t,p}}{3f_{t,n}} \\ F &= \frac{f_{t,p}}{f_{c,p}} - \frac{1}{2}\label{eq3_6}\tag{3.6} \end{align}\]
These formulae have been chosen such that the stresses are normalised by the uniaxial tensile strength parallel to the bed joints.
3.3 NOTES ON IMPLEMENTATION OF CONE FORMULA
Consider the cone formula in the general form:
\[ A_{xx}x^{2} + 2A_{xy}xy + A_{yy}y^{2} + A_{zz}z^{2} + 2B_{x}x + 2B_{y}y + C = 0\label{eq3_7}\tag{3.7} \]
Firstly, select \(C = 1\) in order that inserting any set of coordinates into the left hand side of the equation will return a positive number if it is inside the cone, and a negative number if it is outside the cone. Every load factor calculated must be categorised dependent on whether the scalable load is bringing the stresses into the failure surface or taking the stresses out of the failure surface (i.e. whether to categorise the load factor, \(\lambda\), as a \(\lambda_{t}\) load factor or a \(\lambda_{c}\) load factor – see section 2.2.3). This can be determined by increasing and then decreasing the calculated \(\lambda\) by a very small amount to see if the value returned by the cone formula changes from positive to negative or vice versa. Note that in the event of a negative \(\lambda\), the effect is reversed.
Next, consider the slice of the cone in the x-y plane, i.e. where \(z = 0\).
\[ A_{xx}x^{2} + 2A_{xy}xy + A_{yy}y^{2} + 2B_{x}x + 2B_{y}y + 1 = 0\label{eq3_8}\tag{3.8} \]
This is the formula for a hyperbola in Cartesian co-ordinates. A special case of the hyperbola – the degenerate hyperbola – consists of the two required intersecting lines. For this equation to form a degenerate hyperbola, the following must be satisfied:
\[ \begin{vmatrix} A_{xx} & A_{xy} & B_{x} \\ A_{xy} & A_{yy} & B_{y} \\ B_{x} & B_{y} & 1 \\ \end{vmatrix} = 0 \label{eq3_9}\tag{3.9} \]
As the cone is formed by rotating two intersecting lines about their axis of symmetry, the cone will not terminate at its tip, but also form a second mirror-image cone. When calculating load factors, solutions on the second mirror-image cone should be ignored.
As the failure surface can be built up from a number of intersecting cones, it is necessary to ignore any solutions found that do not lie on the inner surface. Suppose the following solution is found on a particular cone:
\[\begin{align} x &= p \\ y &= q \\ z &= r \label{eq3_10}\tag{3.10} \end{align}\]
The calculated values for \(x\) and \(y\) are then inserted into the cone formula for every other cone and solved for the positive value of \(z\).
\[ z = \sqrt{\frac{-1}{Az}(A_{xx}p^{2} + 2A_{xy}pq + 2A_{yy}q^{2} + 2B_{x}p + 2B_{y}q + 1)}\label{eq3_11}\tag{3.11} \]
If a value of \(z\) is found such that \(z < \left| r \right|\) or \(\textrm{Im}(z) \neq 0\) then that the solution \((p, q, r)\) is outside the failure surface and should be ignored. Note that this relies on the assumption that the failure surface is convex.
As two crack directions are considered for each integration point, it is sensible to consider two separate failure surfaces for each integration point. Each failure surface can shrink independently of the other as the crack it represents grows. As a result of this, each crack only makes use of half of its failure surface and solutions in the opposite half should be discarded.
3.4 NOTES ON IMPLEMENTATION OF THE QUADRATIC FORMULA
There a few caveats to watch out for when solving the quadratic formula using standard floating point arithmetic. Consider the equation:
\[ \lambda = \frac{-b \pm \sqrt{b^{2} - 4ac}}{2a}\label{eq3_12}\tag{3.12} \]
Firstly, when \(a = 0\), a singularity will occur, and instead the following should be used:
\[ \lambda = \frac{-c}{b}\label{eq3_13}\tag{3.13} \]
Secondly, when \(4ac \ll b\) a rounding error will occur for the positive solution. Consider the equation rewritten in the following form:
\[ \lambda = \frac{-b}{2a}(1 - \sqrt{1 - \epsilon}) \]
where
\[\epsilon = \frac{4ac}{b^{2}}\label{eq3_14}\tag{3.14}\]
Taking the first few terms of the binomial expansion:
\[ \lambda = \frac{-b}{2a}(1 - 1 + \frac{1}{2}\epsilon - \frac{1}{8}\epsilon^{2} + O(\epsilon^{3}))\label{eq3_15}\tag{3.15} \]
The source of the rounding error, the \(1 - 1\), can be eliminated, and the binomial expansion taken to a sufficient number of terms.
\[ \lambda = \frac{-b}{2a}(\frac{1}{2}\epsilon - \frac{1}{8}\epsilon^{2} + O(\epsilon^{3}))\label{eq3_16}\tag{3.16} \]
Section 2.2.3 described the scheme for selecting the critical load factor, \(\lambda\). As it stands, no regard is given to the stress history of the structure when selecting the critical load factor. In reality, after cracking occurs at a single integration point, the stress will be redistributed to surrounding region. However, Sequentially Linear Analysis does not consider whether this is feasible before starting the next analysis as it will restart from a state of zero stress. The criteria for selecting \(\lambda\) must be extended to also take account of the critical \(\lambda\) from the previous stage of the analysis.
Figures 5 to 10 list all the possible cases for various arrangements of the critical \(\lambda_{t}\) (the minimum load factor required to fail an integration point inside the failure surface) and the critical \(\lambda_{c}\) (the minimum load factor required to bring all integration points outside the failure surface back inside). For each case, the rule for selecting the critical load factor, \(\lambda_{i}\) and the integration point to crack, \(e_{i}\), is shown dependent on whether the previous load factor, \(\lambda_{i-1}\), is smaller than the two critical load factors, greater than the two critical load factors, or lies somewhere between the two.
For each case, if no stress redistribution is taking place, the new load factor, \(\lambda_{i}\), will become either \(\lambda_{t}\) or \(\lambda_{c}\) as before. Otherwise, if stress redistribution is taking place, the new load factor will roll over from the previous stage and remain \(\lambda_{i-1}\). The integration point to crack, \(e_{i}\), can only be the point corresponding to one of the two critical load factors, \(e_{t}\) or \(e_{c}\).
Consider Case 1 (Figure 5), this case represents the arrangement that will cover the majority of analysis steps. Without the stress redistribution, \(\lambda_{t}\) would be selected as the critical load factor every time. This is acceptable if the previous load factor was greater than the calculated \(\lambda_{c}\), otherwise the problem of stress redistribution arises: the stress state cannot rise from \(\lambda_{i-1}\) to \(\lambda_{t}\) without passing through a region where the \(e_{c}\) integration point is outside of the failure surface. The new criteria distinguish between these two scenarios in order to produce more accurate results.
In this chapter, two laboratory tests of masonry structures are investigated. Both of these tests have been selected as shear failure plays a significant role in determining their capacity. In section 4.1 a destructive test of a simple wall panel under combined compression and shear is reviewed and the results of simulations using the various Sequentially Linear Analysis failure criterions described in this thesis are compared with one another and with the experimental data. In section 4.2 a laboratory test of a full scale masonry façade under a pseudo-seismic load is reviewed and the experimental data are compared with the old uniaxial model and the new biaxial anisotropic model. In section 4.3 the effects of the new stress redistribution criteria on the façade simulation are investigated. Finally section 4.4 discusses the results and proposes further improvements.
Various destructive tests have been carried out on masonry which provide benchmarks for computational modelling strategies. Anthoine et al. performed laboratory tests on simple 1.0m wide rectangular masonry walls12. One such test involved a 1.35m tall wall under a uniformly distributed 150kN compressive load. The top of the wall was fitted to a horizontal actuator to cyclically load the structure in shear until failure. The experimental setup is shown in Figure 11a. This simple example neatly illustrates the differences between various approaches.
Anthoine et al.12 also performed a non-linear iterative finite element analysis of this structure and the results are plotted along with the cyclic test envelope in Figure 11b. First note, the analysis fails as the structure begins to soften, as is common with non-linear iterative methods. Second, note that the ultimate load is overestimated. The uniaxial stress-strain response of the elements have failed to take account of biaxial effects – a reduced shear strength under compression.
Sequentially Linear Analysis, without any of the improvements detailed in this report, was applied to Anthoine’s shear wall test12 using the material properties in Table 3. The Young's modulus, \(E\), tensile strength, \(f_{t}\), and fracture energy, \(G_{f}\), were estimated by scaling from a diagram of the uniaxial stress-strain response used by Anthoine et al.12 for the non-linear iterative finite element analysis. The value of Poisson's ratio and shear retention factor were not stated so typical values were taken. The crack bandwidth, \(h\), for a one Gaussian integration point triangular finite element is found empirically8 to be the average length of the sides of the element. The current implementation of Sequentially Linear Analysis makes assigning unique crack bandwidths to every element very time consuming, as such an average value across all elements of 102 mm was chosen.
| \(E\) | \(\nu\) | \(f_{t}\) | \(G_{f}\) | \(\beta\) | \(h\) |
|---|---|---|---|---|---|
| 1760 MPa | 0.2 | 0.15 MPa | 0.2 J mm-2 | 1.00 | 102 mm |
Table 3: Material Properties for Uniaxial SLA Test
A fixed uniformly distributed vertical load of 150kN was applied to the top of the structure along with a scalable uniformly distributed horizontal shear force. The model was built-in at the base, and the top of the structure was assumed to remain planar, with a constant horizontal displacement along its length. The mesh and loading arrangement are shown in Figure 12. The load displacement diagram and simulated crack patterns are shown in Figure 13. Although Anthoine et al.12 did not provide crack patterns from the test for comparison, they remain in line with expected shear crack patterns. Note rocking failure begins to occur, however shear failure dominates.
Figure 13 also shows that Sequentially Linear Analysis has no problem modelling the post-critical behaviour of the structure and effectively captures the snap-back behaviour where the non-linear iterative analysis failed, however, as with the non-linear method, SLA still overestimates the strength of the structure when biaxial effects are not considered.
Anthoine's non-linear analysis12 predicted a slightly higher peak load than the Sequentially Linear Analysis simulation. The material properties used by Anthoine et al.12 for the analysis were poorly specified (described only on a stress-strain diagram) and as such are likely to differ from the material properties used in the SLA simulation. Additionally, Anthoine et al.12 used a square mesh to model the structure as opposed to the Delauney mesh used with Sequentially Linear Analysis. Square meshes are particularly bad at modelling diagonal crack propagation. The cracking will tend to align itself along the grid pattern, therefore a choice of mesh that doesn't allow smooth crack propagation in the required direction will usually overestimate capacity.
The effects of using the isotropic biaxial failure surface (defined in section 3.1) derived from the unreinforced concrete material tests was then investigated. The same material properties were used as with the uniaxial simulation, the new failure surface selects a uniaxial compressive strength of 10 times the tensile strength. The results are shown in Figure 14.
The biaxial model underestimates shear strength by neglecting the bedding direction of masonry – the shear strength of concrete tapers off much more quickly with increasing compression than for masonry. Note from the crack pattern in Figure 14b that failure is almost entirely dominated by shear and very little rocking failure occurs – in contrast with the uniaxial crack pattern where a large amount of rocking failure was predicted. Modelling masonry as behaving similarly to concrete of comparable strength is therefore a poor approximation.
The shear wall example is now investigated using the anisotropic failure surface (defined in section 3.2) with the same material properties and the following uniaxial strengths:
| \(f_{t,p}\) | \(f_{t,n}\) | \(f_{c,p}\) | \(f_{c,n}\) |
|---|---|---|---|
| 0.15 MPa | 0.10 MPa | 2.49 MPa | 2.96 MPa |
The value for \(f_{t,p}\) is the same as the tensile strength used for the uniaxial test. However, no laboratory test data were available on the other axial strengths, these have been scaled as per typical ratios of the axial strengths suggested by the Dhanasekar et al.11 panel test review. The results are shown in Figure 15.
The anisotropic model does an excellent job of estimating the strength of the structure, peaking at 84 kN, slightly higher than the laboratory test. An overestimate is to be expected due to the monotonic loading used in the model compared with the cyclic test. The simulation also predicts the peak to occur at a smaller displacement than the test, this could be due to the cyclic test reducing the overall stiffness of the structure. Alternatively, it could be due to the elastic constitutive model not consider the directional nature of Young’s modulus in the anisotropic material. More of the failure in the anisotropic model is concentrated in rocking rather than shear compared with the isotropic model. This reflects the increased shear strength and reduced normal tensile strength of the anisotropic model.
A full scale two-story unreinforced masonry building was destructively tested under cyclic loading by Magenes et al.13 This experiment provides a benchmark for the computational modelling of masonry structures.
The geometry and loading arrangement are detailed in Figure 16a. The structure consists of four walls marked A to D. Walls A, B and C were connected by interlocking bricks whereas wall D was isolated from the rest of the structure. For simplicity, this report only investigates the changes effected by the new failure surface on wall D.
The mesh is shown in Figure 16b. Hatched elements were forbidden from cracking due to lintels and experimental apparatus in these portions of the building. A dead load of 248.4 kN was applied to the first floor and a dead load of 236.8 kN was applied to the second floor. Loads were applied via a series of slender steel beams that connected wall B and wall D. The beams had negligible horizontal stiffness so no horizontal forces could be applied through the floor. Two equal horizontal live loads were applied at the two floor levels to simulate a seismic load monotonically. The material properties in Table 4 were utilised for the uniaxial simulation.
| \(E\) | \(\nu\) | \(f_{t}\) | \(G_{f}\) | \(\beta\) | \(h\) |
|---|---|---|---|---|---|
| 1410 MPa | 0.2 | 0.10 MPa | 0.05 J mm-2 | 1.00 | 331 mm |
Table 4: Material properties for façade simulation.
| \(f_{t,p}\) | \(f_{t,n}\) | \(f_{c,p}\) | \(f_{c,n}\) |
|---|---|---|---|
| 0.100 MPa | 0.075 MPa | 2.100 MPa | 2.500 MPa |
Table 5: Axial strengths for anisotropic façade simulation.
For the anisotropic model, the four tensile strengths are shown in Table 5. The tensile strength parallel to the bedding joints, \(f_{t,p}\), is the same as the uniaxial tensile strength from the uniaxial model. Magenes et al.13 noted the normal compressive strength, fc,n, was between 2 and 3 MPa, as such 2.5 MPa was taken as an average. The other two strengths were taken as similar ratios to 0.1 MPa and 2.5 MPa as those found by Dhanasekar et al.11 in the panel test review.
For further details of the finite element model, the reader is referred to DeJong et al.7
The results of the simulation are shown in Figure 17. In order to mimic the effects of cyclic loading, the results have been copied, reversed and superimposed. The anisotropic model does not appear to reduce the strength of the structure by a particularly large amount – both the isotropic and the anisotropic models are still overestimating the base shear capacity by a significant margin. This could be due to rocking or sliding failure dominating and the overestimation arising from other factors, but it is most likely due to the square mesh having a detrimental effect on shear crack propagation.
In order to investigate this, each element was subdivided into two right-angled triangular elements each with one Gaussian integration point, the crack bandwidth for each element adjusted was to 267mm (the average length of the sides). The elements were divided in such a way that the expected shear crack pattern could propagate easily from element to element. The new mesh is shown in Figure 18 and the results using this mesh are shown in Figure 19. Various crack patterns for different displacements are shown in Figure 20.
The new mesh greatly improves the prediction of the structure's peak and residual capacity. Although, as with the Anthoine panel test (section 4.1), the simulation predicts these loads at much smaller displacements than was found in the test. This is illustrated by the considerably larger amount of cracking seen in the simulations than the laboratory test at comparable displacements (Figure 20). This could be due to the cyclic loading reducing the overall stiffness of the structure. Alternatively, although the model has been extended to include the anisotropy in the failure surface, the anisotropy in Young's modulus and Poisson's ratio are not taken into account. The stiffness of the material in the normal direction should be reduced slightly.
The effect of the stress redistribution described in section 3.5 is now investigated. The façade simulation was run with and without the effects on stress redistribution enabled; the load-displacement responses are shown in Figure 21.
The response is largely unchanged aside from some small effects in the extreme post critical region: with the stress redistribution enabled, one of the snap-backs occurs at a slightly lower displacement, also, the simulation continues to provide results at very high displacements.
The lower displacement for the third snap-back is good: it implies the new stress redistribution method applied fewer damage increments to the structure in order to reach the peak from the previous snap-back, which in turn implies damage was applied in a less haphazard manner.
The results at very high displacements are of questionable value as sufficient damage has been applied to the structure by this point for sliding failure to occur. The model only continues to return results as an artifact of the finite element representation of fully cracked elements. The finite element software does not allow elements of zero stiffness so fully cracked elements are given a small but non-zero stiffness. As such, regardless of the amount of cracking in a model, elastic solutions will continue to be found even after a mechanism should have formed.
The two examples in this report were selected as shear played an important role in their failure. Shear failure occurs in masonry panels under large compressive forces with small aspect ratios (e.g. a height:width ratio of around 1.0 to 1.5). Both tests involved large compressive forces, the Anthoine et al.12 panel test had an aspect ratio of 1.35 and the central pier of the Magenes et al. façade test13 had an aspect ratio of 1.29. In both cases the new biaxial anisotropic model predicted much more accurately the peak loads when compared with the previous uniaxial model, but only when using an appropriate mesh. However, the simulation consistently underestimated the displacements of the structure.
One explanation for the underestimation of displacements could be that the material properties other than strength should also be extended to be anisotropic. The Young's modulus and Poisson's ratio chosen in the simulations do not have a directional component: the reduced stiffness in the normal direction is not taken into account in the simulation. As with the tensile strengths, stiffness properties could be specified in the two axial directions in order to better predict the load-displacement response of the structure. Additionally, the shape of the linear softening curve is simply scaled pro-rata with the reduction of the tensile strength. Specifying directional fracture energies would allow for better prediction of load-displacement response, although it should be noted that this effect will be small: the pro-rata scaling provides a reasonable approximation and the response of the structure is very insensitive to the value of the fracture energy14.
Another explanation for the discrepancy in the displacements was that a cyclic load was modelled as a monotonic load. Cyclic loading effectively doubles the amount of cracking in the structure, although most cracks formed in one direction will have closed up for the opposite direction. Additional cracking will also occur from the hysteretic effects that arise from cyclic loading, e.g. additional damage from the opening and closing of cracks. Cyclically loading a structure using Sequentially Linear Analysis would prove difficult without losing the benefits SLA has to offer. While it would be possible to specify a sequence of load factors for the simulation to work its way through, the linearisation is a poor approximation for elements that are frequently alternating between compression and tension. A cracked masonry element may have a particular reduced stiffness in the tensile region, but if it were to return to the compression region, the crack would close up and the full stiffness regained. Specifying separate stiffnesses for tension and compression would introduce a discontinuity to the stress-strain curve and the advantage of the linear approximation is lost. In monotonic loading it is reasonable to assume that the majority of elements will go into tension or compression initially and remain there; however this is not the case for cyclic loading.
The primary contribution herein is a new framework for specifying biaxial and anisotropic failure surfaces for use with Sequentially Linear Analysis. Furthermore the effects of stress redistribution have been addressed. Comparison of the model with laboratory testing has led to the following conclusions:
Lourenço PB. Computational strategies for masonry structures. Ph.D. dissertation. Delft University of Techonology; 1996.↩︎
Lourenço PB, Rots JG. Continuum model for masonry: Parameter estimation and validation. ASCE J Struct Eng 1998;124(6):642–52.↩︎
Rots, J.G.: Sequentially linear continuum model for concrete fracture. In: Fracture Mechanics of Concrete Structures, eds. R. de Borst, J. Mazars, G. Pijaudier-Cabot, J.G.M. van Mier, A.A. Balkema: Lisse, The Netherlands 2001, 831-839.↩︎
Rots, J.G. & Invernizzi, S.: Regularized sequentially linear saw-tooth softening model. Int. J. for Numerical and Analytical Methods in Geomechanics, 28 (2004), 821-856.↩︎
Rots, J.G.; Belletti, B. & Invernizzi, S.: Robust modeling of RC structures with an “event-by-event” strategy. Engineering Fracture Mechanics, 75 (2008), 590-614.↩︎
DeJong, M.J.; Hendriks, M.A.N. & Rots, J.G.: Sequentially linear analysis of fracture under non-proportional loading. Engineering Fracture Mechanics, 75 (2008), 5042-5056.↩︎
DeJong, M.J.; Belletti, B.; Hendriks, M.A.N. & Rots, J.G.: Shell elements for sequentially linear analysis: Failure of masonry structures under lateral loading. Engineering Structures, 31 (2009), 1382-1392.↩︎
Rots, J.G.: Computational modeling of concrete fracture. PhD-Thesis. Delft University of Technology, 1988. Delft University Press: Delft 1988.↩︎
Kupfer H, Hilsdorf H.K., Rusch H. Behaviour of Concrete Under Biaxial Stresses. ACI Journal; August 1969; 656-666.↩︎
Lee S.K., Song Y.C., Han S.H., Biaxial behaviour of plain concrete of nuclear containment building. Nuclear Engineering and Design 227; 2004; 143-153.↩︎
Dhanasekar M., Page A.W., Kleeman P.W., The failure of brick masonry under biaxial stresses. Proc. Instn Civ. Engrs, Part 2, 1985, 79, June, 295-313.↩︎
Anthoine A, Magonette G, Magenes G. Shear-compression testing and analysis of brick masonry walls. 10th European Conference on Earthquake Engineering; 1995; 1657-1662.↩︎
Magenes G, Kingsley GR, Calvi GM. Seismic testing of a full-scale, two-story masonry building: Test procedure and measured experimental response. Experimental and numerical investigation on a brick masonry building prototype. GNDT Rep. 3.0. Pavia (Italy); 1995.↩︎
Lourenço, P. B. Sensitivity analysis of masonry structures, Proc. 8th Canadian Masonry Symp., Jasper, Canada, 563-574 (1998)↩︎