Selected Problems of Contemporary Thermomechanics

The electrochemical performance of solid oxide fuel cell (SOFC) is significantly influenced by three-phase boundary (TPB) zones in the microstructure. TPB zones are locations where all three phases comprising the microstructure such as the two solid phases and the pore phase are present. Electrochemical reactions such as oxygen reduction occur near TPBs, and TPB density is believed to affect the polarization resistance of the cathode. In this regard, the effect of interface degradation under repeated thermal loading on the mechanical integrity and electrochemical performance of solid oxide fuel cell (SOFC) electrodes is studied through finite element simulations. Image-based 3-D models are used in this study, with additional interface zones at the boundaries between dissimilar solid phases. These interface zones are composed of 3-D cohesive elements of small thickness. The effect of interface degradation on mechanical integrity is studied by subjecting 50:50 LSM:YSZ wt.% cathode models to increasing levels of thermal load from room temperature (20°C) up to operating temperature (820°C). Energy quantities (e.g., strain energy and damage dissipation) for cathode models with and without cohesive interface zones are obtained through finite element analysis (FEA). These quantities are compared using energy balance concepts from fracture mechanics to gain insight into the effects of interface degradation on mechanical integrity.


Introduction
Thermomechanics is a scientific discipline which investigates the behavior of bodies (solid, liquid, and gas) under the action forces and heat input. Thermomechanical phenomena commonly occur in the human environment, from the action of solar radiation to the technological processes. The analysis of these phenomena often requires extensive interdisciplinary knowledge, e.g., thermodynamics, continuum mechanics (solid and fluid), soil mechanics, biomechanics, metallurgy, hydraulics, civil engineering, and materials science and even anatomy, chemistry, meteorology, or hydrology. The wide range of thermomechanics applications depends on the field of science and the areas of knowledge in which phenomena are considered. The description of these phenomena requires not only knowledge of the laws of physics but the use of advanced mathematical apparatus, tensor algebra, and methods for solving differential and integral equations. Thermomechanical phenomena are analyzed using analytical and numerical methods. The analytical solution offers a quicker assessment of the searched values and its dependence on the various parameters, but for more complex problems, they are difficult or even impossible to apply. Some problems can be solved only with numerical methods, of which the finite element method is commonly used, but also methods of boundary elements, finite differences and elementary balances. In addition to the mentioned above methods, one needs to know how to solve complex equation systems (in case of the author's original software) or to possess the ability to handle professional engineer's packages.
Thermomechanics therefore describes a broad category of phenomena. It is a generalization of classical mechanical theory and thermodynamic theory. Currently, thermomechanical coupling is a fully formed issue. Basic dependences and differential equations have been formulated based on mechanical and thermodynamic laws. Numerous methods and algorithms for solving differential equations of thermomechanical coupling have been developed, including the finite element method.
Looking at the development of thermomechanics, we cannot omit scientists who laid the foundation for this area of science. First and foremost, Isaak Newton, the author of the three principles of dynamics [1], an outstanding physicist and mathematician, parallels with G.W. Leibniz who developed the theory of differential and integral calculus. In turn, the development of thermomechanics (and not only) was contributed by Fourier, the creator of the Fourier transform and Fourier series theories, which he used in his fundamental work on the theory of heat conduction [2]. One should also mention eminent scientists, the creator of the law of thermal radiation, Kirchhoff [3] and Maxwell [4,5]. Over the past half-century, a number of
We perform finite element analysis (FEA) of thermal stresses induced in reconstructed SOFC cathode microstructures under thermal cycling. 3-D finite element (FE) models of SOFC cathode microstructures are generated from a stack of 2-D microstructure images. Anode (50:50 wt.% NiO:YSZ) and cathode (50:50 wt.% LSM:YSZ) microstructures have previously been analyzed and validated for thermal stress using finite elements by the authors [17]. This study extends the work by investigating the effect of interface degradation under repeated thermal loading on the mechanical integrity and electrochemical performance of SOFC electrodes through finite element simulations. The effect of interface degradation on mechanical integrity is studied by subjecting 50:50 LSM:YSZ wt.% cathode models to increasing levels of thermal load from room temperature (20°C) up to operating temperatures (820°C). Energy quantities (e.g., strain energy and damage dissipation) for cathode models with and without cohesive interface zones are obtained through FEA. These quantities are compared using energy balance concepts from fracture mechanics to gain insight into the effects of interface degradation on mechanical integrity. The electrochemical performance of SOFCs is significantly influenced by three-phase boundary (TPB) zones in the microstructure. TPB zones are locations where all three phases comprising the microstructure-the two solid phases and the pore phase-are present. Electrochemical reactions such as oxygen reduction occur near TPBs, and TPB density is believed to affect the polarization resistance of the cathode [5,18]. In this study, we hypothesize that degradation of weak interfaces under thermal cycling has an adverse effect on TPB zones in the microstructure of the SOFC, leading to a reduction in electrochemical performance over time. Interface degradation under thermal cycling is implemented in the FE electrode models through a simplified scheme. The scheme consists of five successive monotonic, steady-state heating operations from room temperature (20°C) up to operating temperature (820°C), combined with interface strength and fracture energy degradation in each heating operation. Each thermal loading operation represents the heating phase of a normalized thermal cycle, where one normalized cycle represents 1000 actual thermal cycles. SOFCs have been known to perform with great reliability, showing no reduction in performance, for more than 2 years at a stretch [1]. Thus, it may be reasonable to assume that one normalized heating cycle in this study represents 1000 actual thermal cycles in the operating life of the cell. The interface degradation scheme is explained in detail later.

Image-based finite element microstructure models
3-D FE microstructure models were reconstructed from 41 2-D crosssectional images of cathode microstructures [9,10], using focused ion beamscanning electron microscopy (FIB-SEM). Examples of the 2-D images are shown in Figure 1. These images are of the real cathode microstructures having 50:50 LSM:YSZ composition.
Thermomechanics of Solid Oxide Fuel Cell Electrode Microstructures Using Finite Element… DOI: http://dx.doi.org /10.5772/intechopen.76118 In the SEM cathode images, white (pixel value = 255) represents LSM, gray (pixel value = 127) represents YSZ, and black (pixel value = 0) represents the pores. The original cathode image has 217 pixels (width) × 147 pixels (height). The x-y (in-plane) spatial resolution between pixels is 40.8 nm and the z-spacing between images is 53.3 nm. The cathode images are thus of size 8.85 μm (width) × 6.00 μm (height). Finite element modeling was carried out using the commercial FE software Abaqus [19] for which MATLAB ® [20] was to create input files. The 3-D FE model was reconstructed from the 2-D images using 3-D finite elements (eight-node brick elements). In order to increase computational efficiency for 3-D analysis using finite element method, we reduced the full model to a representative model in which we sacrifice some details of phase geometry but microstructural skeleton that is crucial to stress analysis remains almost unchanged. To provide validation data for this present approach, we used 2-D cathode microstructure images and calculated the three-phase boundary (TPB) density and phase surface area density for the original, full-size 50:50 wt.% LSM:YSZ cathode, as well as the reduced 50:50 cathode. The data agreed reasonably well with those reported in the literature [10,12,21].
The original images of the cathode were simplified by sampling pixels at regular intervals to reduce the image resolution while retaining a detailed microstructure. The simplified cathode images were 3.51 μm (width) × 3.02 μm (height) in size. The depth of the voxel was 53.3 nm. A stack of all the 2-D images was created in the z-direction by using a cell array construct to arrange the images consecutively. An initially blank "buffer" plane was then introduced between consecutive images. The gaps between consecutive images were filled by assigning one eight-node brick element to each voxel. Thus, the geometry of the cathode microstructure was reconstructed in the 3-D FE models, with a step variation in material regions between consecutive images. The volume fractions of phases in each 3-D model were calculated by counting the number of voxels corresponding to each phase (based on pixel value) and dividing by the total number of voxels in the model. Information concerning the material properties, boundary conditions, initial temperature, temperature field and required outputs was also specified in the input file.

Thermomechanical material properties
Finite element analyses of the LSM/YSZ (50:50 wt.%) cathode microstructure models were carried out to investigate thermal stresses due to various temperature fields. The effects of varying phase volume fractions and temperature-dependent material properties on thermal stresses and probability of failure were investigated. The FE model was subjected to fixed boundary conditions. The behavior of the Selected Problems of Contemporary Thermomechanics 8 model with increasing temperature loads was investigated by subjecting the model to eight different spatially uniform predefined temperature fields of magnitude 120, 220, 320, …, 820°C. In each analysis, the initial temperature was specified as 20°C (room temperature), so that the model was subjected to eight different magnitudes of temperature change (ΔT = 100, 200, 300, …, 800°C). Table 1 lists the room temperature material properties used for YSZ and LSM [3,22]. Figure 2 shows the variation of the CTE of YSZ with temperature [22] and the variation of Young's modulus of LSM and YSZ with temperature [4]. The CTE of LSM was assumed constant over the temperature range considered. The room temperature value of the CTE of LSM was used in the FE analyses.

Finite element analysis using cohesive interface model
An interface damage initiation and evolution model is implemented in the 3-D finite element (FE) SOFC electrode microstructures. The thin interface zones between dissimilar solid phases consist of 3-D cohesive elements. The 3-D cohesive elements, which are of negligible thickness compared with the neighboring solid elements, are used to simulate debonding between different solid phases at their interface. The 3-D cohesive elements are inserted between eight-node 3-D linear brick elements belonging to dissimilar solid phases, as shown in Figure 3. Based on the in-plane (x-y) dimensions of the solid elements in the anode (14.0 nm) and cathode (40.8 nm), the thickness of the thin interface cohesive elements is chosen as 1 nm. Interfacial (fracture) energy is used as a reference value for bond strength between LSM and YSZ (cathode). Perfect bonding is assumed between elements belonging to the same solid phase. An algorithm based on the boundary pixel identification scheme is used to identify elements lying on solid-phase boundaries.  [3,22].

Figure 2.
Variation of CTE of YSZ with temperature [16]; variation of Young's modulus of LSM and YSZ with temperature [12].

Traction-separation behavior of cohesive elements
The behavior of the interface cohesive zone elements is described using a traction-separation relation. Such a cohesive zone model (CZM) describes local material separation behavior by relating tractions (t) on the interface surfaces to material separation (δ). The CZM simulates progressive damage in the cohesive zone. The CZM is used to model imperfect bonding in highly porous electrodes of SOFCs. In three dimensions, the traction-separation model for cohesive elements consists of three deformation modes: one normal mode and two in-plane shear modes. In the normal mode, the deformation is purely in a direction perpendicular to the plane of the interface. In the two in-plane shear modes, the deformation is parallel to the plane of the interface. Thus, for 3-D cohesive elements, the traction vector (t) has three components: one normal component (i.e., perpendicular to the plane of the interface), t n , and two in-plane shear components, t s and t t , parallel to the plane of the interface. The separation vector (δ) has three corresponding components, δ n , δ s , and δ t . The nominal strain components can then be written as follows: Here, T 0 is the initial thickness of the interface. If we set T 0 = 1, the nominal strain components become numerically equal to the corresponding separation components. This rescaling can be achieved by defining the interface stiffness (K i ) as follows: E i is the Young's modulus of the interface material and T 0 is the initial thickness of the interface. In this study, the weak interfaces between dissimilar solid phases are modeled using cohesive elements with linear elastic traction-separation behavior up to damage initiation, followed by linear stiffness degradation, to simulate damage evolution. A linear stiffness degradation model has been used by Xiao et al. [23] to study the fracture behavior of ultra-high strength concrete. Thus, the linear damage evolution model may be appropriate for describing the degradation of SOFC electrode interfaces between dissimilar solid (ceramic) phases that exhibit brittle behavior. The normal and in-plane modes of the traction-separation model used in this study are assumed to be uncoupled. The uncoupled linear elastic constitutive behavior of the cohesive elements can be expressed as follows, in terms of tractions t and separations δ (which are numerically equal to strains ε): The uncoupled behavior of the cohesive elements is indicated by the fact that all off-diagonal terms are zero in the interface stiffness (K i ) matrix in the constitutive traction-separation (t-δ) relation given earlier. Figure 4 is a schematic representation of the traction-separation curve for LSM/YSZ interface cohesive elements used in this study. In Figure 4, δ 0 is the critical separation at which damage is initiated in the cohesive element. Similarly, δ f is the final separation at the ultimate failure of the cohesive element.
Based on the approach adopted by Nguyen et al. [24], a quadratic stress interaction equation is used to describe the damage initiation condition: Here, t n 0 , t s 0 , and t t 0 denote the maximum nominal stress when the deformation is either purely in a direction perpendicular to the plane of the interface or in one of the two orthogonal directions lying in the plane of the interface. The Macaulay brackets < > enclosing t n indicate that only tensile normal stress causes damage initiation in the cohesive elements. Purely compressive normal stress causes no damage in the cohesive layer [25]. This equation is used to describe the damage initiation criterion for the interface, accounting for the interaction between the normal (tensile) and in-plane (shear) stresses. Due to lack of detailed experimental data on the strength and stiffness properties of the LSM/YSZ interfaces, the approach used by Nguyen et al. [24] is adopted for the cathode: the normal tensile strength and in-plane shear strengths of the LSM/YSZ interface are all assumed to be equal to the Weibull strength of the weaker phase (LSM), that is, 52 MPa [26,27]. The stiffness of the LSM/YSZ interface is assumed to be characterized by the Young's modulus of   [4]. Stiffness degradation does not occur for the cohesive elements under pure normal compressive stress.
The definitions of damage evolution used in this study are based on the fracture energies of the LSM/YSZ interface (cathode). The fracture energy (per unit interface area) is the area under the traction-separation curve for the interface. The critical energy release rate (G i ) for the LSM/YSZ interface is taken as 7.80 J m −2 [11]. Based on numerical simulation results reported in the literature [11], the interface is assumed to fail in pure Mode I (normal mode) loading conditions. This leads to a mode-mix ratio ψ = 0° based on traction components. The Young's modulus and original thickness of the LSM/YSZ interface are used to calculate its normalized stiffness K i (40 x 10 9 N mm −3 ). The stiffness K i and interface strength are used to calculate the critical separation δ 0 (0.0013 nm). The area under the tractionseparation curve is equal to the fracture energy per unit interface area, 7.80 J m −2 . This is used to calculate the final separation at ultimate failure, δ f (0.3 μm).

Interface degradation scheme
Degradation of the interface and TPB zones between dissimilar solid phases under repeated thermal cycling may be one of the major reasons for degradation of SOFC performance over time. In this study, interface degradation with thermal cycling is simulated using a simplified scheme. Only the heating phase of the thermal cycle is simulated. Steady-state conditions are assumed and transient effects are neglected. The models are subjected to five successive steady-state heating operations. Each analysis utilizes a spatially uniform temperature field to simulate steady-state heating of the model from room temperature (20°C) up to operating temperature (820°C). Progressive interface degradation is simulated by decreasing the stress at which damage begins in the cohesive interface layers in each successive analysis. The critical and final separations for the cohesive elements are assumed to remain unchanged over the five heating cycles, leading to a progressive reduction of both interface stiffness and critical energy release rate, as shown schematically in Figure 5.
Critical energy release rate (i.e., fracture energy per unit interface area) is a measure of the resistance of the interface material to damage [28,29], while stiffness is indicative of the amount of damage (i.e., stiffness decreases as damage progresses). Therefore, both stiffness and fracture energy are assumed to progressively decrease with thermal cycling. Thus, this scheme implements interface degradation in a  manner. Mechanical degradation of the interface is considered, without considering interface degradation due to electrochemical or redox processes. The earlier scheme is only a first approximation to the complex degradation processes that occur at the actual interfaces in a real SOFC electrode. Multiphysics simulations considering both mechanical and electrochemical degradation processes may be an interesting topic for future research.

Study design
Two sets of electrode microstructure models are generated. One set has perfect bonding throughout the entire model and the other has cohesive zones between dissimilar solid phases. Image-based 3-D models similar to those described in previous sections are used here, with one important difference. The models used are derived through image simplification by analyzing (2 × 2) pixel squares in the original images. This scheme increases the effective pixel size but preserves both the overall size of the image and the depth of the voxel, leading to a decrease in the reconstructed volume due to the use of a limited number of images. However, this scheme leads to a reasonably good approximation of the original microstructure. In the present study, a similar approach is used to simplify the model. An algorithm is devised to analyze consecutive pixel squares of (n × n) pixels in each original image. For each such square, the number of pixels belonging to each phase is counted. The phase that contributes the maximum number of pixels to the n × n pixel square is assigned as the final phase of the single pixel at that location in the simplified image. Additionally, in the present approach, preservation of the original volume of the electrode model in the simplified version is implemented by causing an increase in both effective pixel size and effective voxel depth. The volumes of the simplified models used here are 112.18 μm 3 for the anode and 670.35 μm 3 for the cathode. The cathode volume may be compared against the sampled volume of 178 μm 3 for the simplified cathode model. The values of n used in this study are as follows: for the cathode, n = 6, and for the anode, n = 10. These levels of simplification are necessary to achieve reasonable computational time for each FE analysis. A detailed validation study for this new approach based on volume fractions, phase surface area densities, and TPB densities may be an interesting topic for future research. In the present study, the plausibility of this new approach is demonstrated by showing that TPB density values obtained for cathode models lie within physically reasonable limits of experimentally determined ranges for these values.
Increasing levels of thermal load are applied through spatially uniform temperature fields of magnitude 120°C, 220°C,…, 820°C, to simulate steady-state heating from room temperature (20°C) up to operating temperature. Energy quantities (e.g., strain energy and damage dissipation) are calculated from FEA for each ΔT = 100°C, 200°C,…, 800°C. The energy quantities for models with and without cohesive zones are compared to gain insight into the effect of weak interface zones between dissimilar solid phases on the mechanical integrity of the cathode. Energy concepts from fracture mechanics are used to interpret the results obtained. This procedure is performed using the cathode microstructure model, considering the original composition (50:50 wt.% LSM:YSZ). In the second part of the study, threephase boundary (TPB) evolution under repeated, steady-state, monotonic thermal loads is studied using 50:50 wt.% electrode microstructure models. Electrode models with cohesive interface zones are subjected to five successive monotonic, steady-state, heating operations from room temperature up to operating temperature (ΔT = +800°C) using a spatially uniform temperature field. During each successive analysis, the strength of the cohesive interface between dissimilar solid phases is decreased to simulate the effect of interface degradation under thermal cycling. Progressive interface degradation under thermal cycling is implemented using the simplified scheme explained earlier in this chapter. The strains induced in the cohesive layers are calculated using FEA and are used to quantify TPB zone damage using a strain-based criterion, that is, separation in the cohesive layers, as explained later.

Results and discussion
Cohesive elements give rise to numerical convergence and stability issues that are well documented in the literature. Unstable behavior of cohesive elements due to negative stiffness (i.e., softening) in the damage evolution phase causes local instability in the electrode models used in this study. The local stabilization algorithm available in Abaqus is used along with viscous regularization for the cohesive elements to address convergence and stability issues. Small values of the viscosity parameter (μ = 0.001-0.01) are used to ensure that the energy quantities dissipated in regularization and stabilization are small compared with the strain energy in the model. Fixed boundary conditions (BCs) are prescribed by fixing all degrees of freedom (DOFs) at each node on each face of the models. The 3-D cathode model with cohesive zones is subjected to steady-state temperature change from room temperature (20°C) up to operating temperature. Spatially uniform temperature fields of increasing magnitude (120°C, 220°C,…, 820°C) are used to apply increasing levels of thermal loads (ΔT = 100°C, 200°C,…, 800°C). Similar analyses are performed for the 3-D cathode model without cohesive zones. Energy quantities (e.g., strain energy and damage dissipation) are computed from each FE analysis. Figure 6 also shows the interface zones between dissimilar solid phases that are composed of cohesive elements. Figure 7 shows energy quantities obtained from the FE analyses of cathode models with and without cohesive interface zones. The energy quantities are plotted as functions of temperature, which show a comparison between the strain energies of the cohesive model and the non-cohesive model. The energy dissipated due to damage (i.e., the damage dissipation) in the interface zones of the cohesive cathode model is plotted as a function of temperature in Figure 8. It is observed from Figure 7 that the strain energy of the cohesive and non-cohesive cathode models increases with increasing temperature, as expected from the larger thermal stresses induced at higher temperatures, with fixed boundary conditions. Similarly, from Figure 8, it is seen that the energy dissipated due to interfacial damage in the cohesive model also increases with increasing temperature, since more damage occurs at higher temperatures. A significant feature of Figure 7 is that the strain energy of the non-cohesive model is slightly higher than that of the cohesive model at all temperatures. The physical explanation of the earlier observations is as follows. From the first law of thermodynamics and the Griffith energy balance statement for fracture (i.e., damage), we know that a damage process (e.g., cracking) can occur in a system only if this process causes the total energy of the system to either decrease or remain constant. The limiting condition is attained when damage occurs at equilibrium conditions, with the total energy remaining constant [30]: where E is the total energy of the system, Π is the potential energy of the system = U-F, U is the strain energy of the system, F is the work done by external   forces, W s is the work required for creating new surfaces, and A is the damage variable (e.g., crack length).
In the cathode models analyzed, fixed boundary conditions are prescribed. This means that all translational and rotational DOFs at each node on each face of the model are set to zero. Thus, this model can be regarded as a displacement controlled model [30], with the displacement of the faces fixed at zero. For a displacement controlled system, we know that the external work is zero (F = 0). This is confirmed by the numerical values of the external work obtained from FEA, which are on the order of 10 −40 . Thus, for a displacement controlled system, Π = U. Finally, from Irwin's energy release rate analysis [30], we know that for displacement control the strain energy of the system decreases as damage (e.g., cracking) occurs. This explains the lower strain energy of the cohesive model as compared to the noncohesive model at a given temperature.
The second objective of this study is to approximately simulate the effects of thermal cycling on electrode integrity and performance using monotonic, steadystate, thermal loads coupled with an interface strength degradation scheme. In order to achieve this, the electrode models (cathode: 50:50 wt.% LSM:YSZ) with cohesive zones are subjected to five successive analyses. Each analysis represents the heating phase of a normalized thermal cycle, with one normalized cycle being taken equal to 1000 cycles. Only the heating phase of the thermal cycle is simulated in each analysis. Steady-state conditions are assumed and transient effects are neglected. Starting from a uniform initial room temperature of 20°C, the models are subjected to steady-state temperature change up to an operating temperature of 820°C using a spatially uniform temperature field. During each successive analysis, progressive interface degradation is simulated by decreasing the interface strength and fracture energy of the cohesive layers. Two separate interface strength degradation schemes are studied. In the first scheme (Scheme 1), the initial interface strength is assumed to decrease by 5 MPa in each successive cycle. In the second scheme (Scheme 2), the interface strength is assumed to decrease by 10 MPa in each successive cycle. The constitutive thickness of each cohesive element is set equal to 1.0, so that the nominal strain components (ε) calculated using FE analysis are numerically equal to the respective separation components (δ). Figure 9 shows the interface strength reduction schemes employed to study TPB evolution in the cathode under repeated thermal loading, that is, Scheme 1 (Δ t i 0 = − 5 MPa per cycle, where t i 0 = interface strength) and Scheme 2 (Δ t i 0 = − 10 MPa per cycle). The numerical values of interface strength and fracture energy used in Scheme 1 (Δ t i 0 = − 5 MPa per cycle) and Scheme 2 (Δ t i 0 = − 10 MPa per cycle) are given in Table 2. Pitakthapanaphong and Busso [11] have performed numerical simulations to determine the critical energy release rates (G i ) for several material interfaces encountered in SOFC systems. Numerically, these values vary from 7.80 J m −2 for LSM/YSZ interface to 38.70 J m −2 for LSM-LSCoO interface. From these data, it can be concluded that the interface fracture energy values used in this study ( Table 2) are physically reasonable. Cathode models with fixed boundary conditions (BCs) are used in this study.
The models are subjected to repeated, monotonic, steady-state heating as explained earlier, and the thermal stresses and strains are calculated through FEA. The strains induced in each cohesive element are written to the output data file associated with each analysis, along with strain energy and damage dissipation values. The strains are used to calculate the evolution of three-phase boundaries (TPBs) with thermal cycling, as explained in the next subsection. The energy quantities are used to study the mechanical degradation of the overall model. The TPB length at each temperature is calculated as follows. First, an electrode model without cohesive zones is created by image stacking (using a cell array data structure in MATLAB). The original 2-D images used to reconstruct the 3-D microstructure are stored as pixel value arrays in the cell array. In each 2-D image, boundary pixels are identified by pixel value comparison of each pixel with its eight adjacent in-plane neighbors. Boundary pixels are stored as black pixels (pixel value = 0) at corresponding locations in an initially empty "boundary point" cell array of identical size as the image cell array. Interior (i.e., non-boundary) locations are stored as white pixels at corresponding locations in the boundary point cell array. Using the boundary point cell array, locations of boundary pixels in the image array are identified as black pixels in the boundary point array, and the pixel value of the boundary pixel at the corresponding location in the image array is checked. The pixel values of the eight pixels adjacent to that boundary pixel are also checked to test whether all three phases are present at that boundary location. If so, the location is stored as a three-phase boundary (TPB) point (black pixel) in another initially empty "TPB" cell array. All other non-TPB locations are stored as white pixels in the TPB array. After TPB locations have been identified, the total number of TPB points is counted, and multiplied by the z-distance between images, to obtain the total undamaged TPB length (μm) in the electrode models. The total undamaged TPB length (μm) is divided by the volume (μm 3 ) of the 3-D electrode model to obtain the TPB density (μm −2 ) in the model. Damage accumulation over five normalized cycles is calculated by determining the failed TPB cohesive elements in each cycle that had not failed in the previous cycle. The difference between the number of failed TPB elements in the (n + 1) th cycle and in the n th cycle is used to quantify the additional TPB loss in each  successive cycle, and the cumulative TPB loss over N such cycles is thus determined by successively adding the new damage in each cycle to the cumulative total damage over the previous (N − 1) cycles. The model used in this study for evaluating cumulative TPB damage with thermal cycling may be regarded as a simple model for calculating cumulative fatigue damage. A detailed mathematical model for predicting the lifetime of planar SOFCs subjected to thermal cycling has been developed by Liu et al. [31]. That model uses Paris' law and cracks nucleation concepts to derive expressions for damage distribution in the interfacial layers of planar SOFCs under thermal cycling. It predicts that the number of cycles required for failure will decrease with increase in electrolyte thickness and electrode porosity [31]. Figure 10 depicts the change in the strain energy content of the cathode with thermal cycling. Figure 10 shows that for both interface degradation schemes, the strain energy of the model progressively decreases over five normalized thermal cycles. This is expected from the progressive mechanical degradation of the interfaces within the model with thermal cycling.  Such progressive interface degradation leads to a cumulative dissipation of energy due to damage and hence to a progressive decrease in strain energy content of the model. This is further illustrated in Figure 11, which shows the cumulative damage dissipation with thermal cycling. As expected, Scheme 2 (Δ t i 0 = − 10 MPa per cycle) leads to lower strain energy content, and larger cumulative damage dissipation, than Scheme 1 (Δ t i 0 = − 5 MPa per cycle). Figure 12 shows the evolution of the TPB density with thermal cycling for the anode and cathode models considered in this study. Figure 12 simultaneously compares the TPB density evolution for cathode versus anode and for Scheme 1 (Δ t i 0 = − 5 MPa per cycle) versus Scheme 2 (Δ t i 0 = − 10 MPa per cycle). The TPB density of cathode model decreases over the five successive normalized thermal cycles (∆T = +800°C), as expected from the progressive degradation of the interface with thermal cycling. As expected for the cathode model, Scheme 2 (10 MPa degradation per cycle) leads to a lower final TPB density than Scheme 1 (5 MPa degradation per cycle), after five normalized heating cycles. The TPB density of the cathode undergoes a large reduction (5.36% for Schemes 1, 11.18% for Scheme 2). Wilson et al. [10] have cited TPB density values lying in the range of 1.7-6.5 μm −2 for cathodes and have reported cathode TPB density values in the range 6-9 μm −2 [12]. From these experimental data, it may be concluded that the TPB density values obtained in this study for the cathode (2.138-2.407 μm −2 ) are physically reasonable.

Conclusion
2-D images of solid oxide fuel cell (SOFC) cathode microstructures (50:50 wt.% LSM:YSZ) are used to perform 3-D finite element models. The effect of LSM/YSZ interface degradation under repeated thermal loading on mechanical integrity and electrochemical performance of SOFC electrodes is simulated by implementing a simplified damage scheme in cathode FE models with cohesive interface zones. The cathode model is first subjected to increasing levels of thermal load using spatially uniform temperature fields. Energy quantities for models with and without cohesive interface zones are obtained through FEA. These quantities are compared using energy balance concepts from fracture mechanics to gain insight into the effects of interface degradation on mechanical integrity. The effect of interface degradation with increasing temperature on the mechanical integrity of a displacement controlled cathode microstructure model is clearly seen as a reduction in strain energy   of the model due to damage dissipation. The evolution of three-phase boundary (TPB) zones in electrode microstructure models with thermal cycling [32] is studied by implementing an interface damage scheme that includes reduction of both interface strength and fracture energy. It is found that TPB density decreases over a number of normalized heating cycles. Degradation of the mechanical integrity of the cathode model under repeated thermal loading is also observed, in the form of progressively decreasing strain energy content due to cumulative damage dissipation with thermal cycling. These observations indicate that interface damage may be a major mechanism responsible for SOFC performance degradation over time. [1] [20] MATLAB ®R2010a, the MathWorks, Inc., Natick, Massachusetts [21]

Introduction
The thixotropic substances are considered structured fluids and have a rheological behavior (time-dependent) guided by their structural nature [1][2][3][4][5]. In terms of models (constitutive systems), this behavior is represented by a couple of timedependent equations. These equations connect the micro-structural characteristics to the rheological behavior. In many works, this system of equations is based on a qualitative way, without any formal justification on the physical principles. In this chapter, a different approach is presented based on some well-established physical principles (Smoluchowski's coagulation theory [6][7][8][9][10], Brownian motion [6][7][8][9] and de Gennes reptation model [8,11,12]).
The rheological properties evolution, in many thixotropy models [2,5,[13][14][15][16][17], is presented, generically, in terms of a couple of equations, as follows: τ ¼ τ λ; _ γ ð Þ, constitutive relationship based on linear viscoelastic models; (1) where τ is the shear stress (from Cauchy stress measure), _ γ is the shear rate (infinitesimal strain measure), λ is the structural parameter (a positive scalar quantity associated to the structural level of substance) and G _ γ is a functional form. In large number of works [5,[18][19][20][21], Eq. (1) is represented by a linear viscoelastic model/mechanism (Maxwell, Jeffreys, Kelvin-Voigt, etc.) or combinations of that (in parallel or in series). However, it exists a problem with these approaches, they do not consider, the time/micro-structural dependence of the shear modulus (G) and viscosity (η μ and η ν ) in physical principles from which the constitutive equation system is obtained. Differently, the approach presented in [8] (the modified Jeffreys model) does not present this mistake, considering time/micro-structural evolution in the constitutive system development, and in this way, showing more coherence to represent thixotropic substances behavior.
Eqs. (1) and (2) connect micro-structural characteristic and viscoelasticity nature of the thixotropy. In the specialized literature, a considerable number of micro-structure types [22][23][24][25][26][27][28][29] can be found, although only one single microstructure type is considered in the development of the modified Jeffreys model. In this way, it is natural to imagine that it can exist more than one micro-structure type in a same thixotropic substance. In this chapter, a new perspective in terms of constitutive approach for thixotropic substances, with apparent-yield-stress nature and based on the modified Jeffreys model [8], is presented in terms of two different types of micro-structure (i.e., two structural parameters), and some of their properties are investigated in a rheological and thermodynamic sense. In the presented approach, each micro-structure is associated to a specific mechanism considered in the model (Maxwell-like element in parallel with a pure viscous Newtonian-like element), that is, one micro-structure type is associated to the Maxwellian element and another one to the pure viscous element. In this way, it is important to stand out that each of them (mechanisms) can react, under the shear loads, in distinct (but integrated) forms one from the other.
Aiming to extend the modified Jeffreys model for two distinguished microstructures, in the next sections, the constitutive model is formally presented in terms of equations, the analysis of its thermodynamics consistence is done, points related to the characterization of the transition region are discussed and some illustrative numerical simulations of rheological tests are presented as well.

Main ideas on two types of micro-structures approach
The approach proposed in this work is based on a Jeffreys model [8] for the thixotropic system representation. In this case, a sketch of the model (Maxwell element in parallel with viscous element) is presented in Figure 1.
It is important to point out the presence of two micro-structure types (λ μ and λ ν ) that are intimately related to the behavior of the shear modulus G ¼ G λ ν ð Þ, and viscosity coefficients η ν ¼ η ν λ ν ð Þ and η μ ¼ η μ λ μ À Á . In this way, under a certain shear load (τ), for the Maxwell element (: ν ), one has where _ γ ν v is the viscous strain rate and _ γ ν e is the elastic strain rate, respectively, and the total strain rate for the Maxwell element is _ After some algebraic manipulation, it follows: From the viscous element (: μ ), and remembering that _ It is important to stand out that the time variation of the shear modulus and the viscosity coefficient are taken into account, which is in agreement with the expected physical behavior of thixotropic fluids [8,30,31].
It is important to keep in mind the fact that the parameters λ μ (0 ≤ λ μ ≤ 1) and λ ν (0 ≤ λ ν ≤ 1) are used to characterize the structural level of the substance. The structural parameter closer to 1 implies highly building-up state of micro-structure, and when it is closer to 0, the micro-structure is close to a fully broken-down state, in respective mechanism (: ν or : μ ). In this sense, it is presented the two microstructures version of constitutive model based on the work [8]. Figure 1.
Sketch of the two micro-structured thixotropic constitutive model.

On the thermodynamic consistence
In this section, the thermodynamic consistence of constitutive model is analyzed in terms of the entropy-producing processes [33] associated with the concept of natural configuration [8,[34][35][36]. The main objective is to verify the consistence of the constitutive equation to the Clausius-Duhem inequality.  Configuration spaces.
and the extensive thermodynamic Helmholtz potential (Ψ) can be represented as where follows from chain rule Ψ, ξ i Bk ), that is, for an appropriate rotated natural configurationk n t ð Þ , it follows [8]: where D is the symmetric part of velocity gradient (L). Returning to Eq. (20), one has Then, from the Clausius-Duhem inequality follows the specific rate of dissipation Ξ 27 On the Thermodynamic Consistency of a Two Micro-Structured Thixotropic Constitutive Model DOI: http://dx.doi.org/10.5772/intechopen.75987 where τ is the Cauchy stress tensor. This approach makes possible the analysis of the constitutive law and the nonnegativity of the entropy production in the context of irreversible thermodynamic processes. In this sense, it is assumed that total rate dissipation Ξ can be divided into three parts, each one associated to a specific process, as follows: where Ξk t is the dissipation rate due to changes ink t B ð Þ), Ξ λ is the dissipation rate due to changes in micro-structure and Ξk n t ð Þ is the dissipation rate due to changes ink n t ð Þ B ð Þ. Following similar steps to the work [8], one has where Thus, one has for each ξ i -relaxation mechanism.
where ðÞ ⊙ ¼ _ ðÞ � LðÞ � ðÞL T is the "Oldroyd time derivative". Considering the two relaxation mechanisms, with r ν ¼ r Therefore, as τ ¼ τ ν þ τ μ one has Note that the isochoric motion constraint was not required in this analysis. In this way, associating each relaxation mechanism to a micro-structure, it is important to point out that depending on the nature (level of complexity) of the analyzed thixotropic substance, the reasoning presented here can be extrapolated for more than two relaxation mechanisms (i.e., more than two micro-structure types) aiming a consistent/coherent (in thermodynamic and rheological senses) approach of the model.

On the transition region
The objective of this topic is to analyze the transition/yielding region criteria under a theoretical and formal point of view. In this sense, it is important to stand out that the constitutive equation can be rewritten in a more appropriate form as follows: where It can be seen that the form of Eq. (36) is quite similar to the standard Jeffreys model. The quantities Π, _ Γ, η, θ 1 and θ 2 can be interpreted as a dimensionless stress, dimensionless strain, dimensionless apparent viscosity, relaxation time and retardation time, respectively. Note that since Eq. (36) corresponds to the standard Jeffreys model. In this way, it is clear the relation between the nonlinear viscoelastic region (associated to the thixotropic effect) and _ λ. Many works propose that an indicative for the beginning of the transition/yielding region is a region where occurs the modification in the behavior of the thixotropic substance, from the linear viscoelastic region to nonlinear viscoelastic region [37][38][39][40][41] associating an specific strain value for that. However, this specific value is not constant [41][42][43][44][45].
In this sense, taking into account the mapping in terms of variables in the actual configurations (x ¼ x i g i ) and the respective displacements (u ¼ u i g i ), it follows:  , and in this context, one has the Jacobian (J) and with where δ ij kl and δ ijk rst are the generalized Kronecker delta [46,47]. Thus, the relationship between dV (differential of volume in the reference configuration) and dv (differential of volume in the actual configuration) can be described in the following way: Taking into account that the yield region can be characterized as a transition region. This region, in fact, can be treated as a singularity region. As the continuity axiom holds a suitable asymptotic behavior [48][49][50], the concept of transition phenomenon has an asymptotic nature, and can be treated by limiting approaches. Physically, the transition region (transition/yielding) can be noted when a substance loses part of their original intrinsic properties, and from this point a new set of properties are raised in this substance, resulting in a new behavior. For example, in an elasto-plastic transition behavior, the transition state is related to a strain ellipsoid degeneracy, in a geometric sense, which can be transformed to an infinite cylinder, point or a pair of planes [51].
Considering the abovementioned lines, it can be imagined the nonlinear viscoelastic region as a direct mapping image from the linear viscoelastic region. Thus, the Jacobian (J), of this mapping, should translate the singularity on the transition region vicinity as an asymptotic extremely large deformation of the original microscopic element. In other words, taking into account Eq. (47), it can be seen that for a finite initial volume dV, withJ ! 0 (singularity on the neighbor of the transition region) implicates in dv ! ∞, and in this way, it can write for the transition region the following condition: or in an asymptotic sense on the transition region (t.r.) lim ε ! t:r: It is important to stand out that this condition was stated in terms of the strain invariants, which can be rewritten in terms of the stress invariants or in terms of energy, by the constitutive relationships. Note the difficulty involved in expressing the strain tensor invariants in terms of the stress tensor invariants, due to the constitutive model's form. Remark 1. If it is taken the time rate of Eq. (48), one has Supposing the case that the strain rate is obtained as a response to known stress load, returning to Eq. (36), it follows: with Thus, one has where Note that and in this way, from Eq. (53), it follows: Selected Problems of Contemporary Thermomechanics or in other words where L ¼ ∇v ("Eulerian gradient"), and where v is the velocity field. Consequently, since λ μ and λ ν are taken as constants (pre-transition/yield region). In this way, the criteria Eq. (48) and/or Eq. (50) can be appropriately analyzed for a stress load excitation.

Numerical results
In this section, some illustrative numerical results are presented. Two types of rheological tests are analyzed, the constant shear rate test (CR) and the constant shear stress test (CS). In both cases it is necessary to define consistent initial conditions in relation to physical and mathematical principles. The finite difference method was implemented in MATLAB, with the following set of parameters: ℵ ¼ 1, η 0 ¼ 0:08 Pa:s, Pa, t ν ς ν ¼ 10 4 s �1 and t μ ς μ ¼ 10 3 s �1 . It is important to point out that for both tests it was used Δt ¼ 10 �3 s. It is also important to comment that time-steps Δt ¼ 10 �4 s, 10 �5 s and 10 �6 s were analyzed, but modifications in responses were not noted. The tolerance used for the Newton's method was 10 �6 . Numerical results for real thixotropic substances (numerical/experimental comparative responses, regularization methods associated to nonlinear identification parameter problem, etc.) can be found in [52].

Constant shear rate test
In this section, it is considered the constant shear rate test, that is taken into account load conditions as _ γ t ð Þ ¼ H t ð Þ _ γ 0 , where H t ð Þ represents the standard Heaviside function and _ γ 0 it is a positive real constant. It is reasonable to think that in the begin of the test, the micro-structure is in a fully structured state (λ μ 0 ð Þ ¼ λ ν 0 ð Þ ¼ 1). Figure 3 shows an interesting aspect, for the loads _ γ 0 ¼ 10 �2 s �1 and 10 �1 s �1 the behavior is close to a purely viscoelastic response. These behaviors are in agreement with Figures 4 and 5, for the same loads. Note that these two shear rates correspond to a low level of modification in λ μ and λ ν . In these both cases, the stress of steady state is the maximum stress reached.  It is important to point out that for _ γ 0 ¼ 1 s �1 the thixotropic effect can ever be noted in stress response (Figure 3). The loads _ γ 0 ¼ 10 s �1 and 100 s �1 show a typical thixotropic behavior (Figure 3). These responses are related to considerable modifications in the structural nature of the substance, as can be seen in Figures 4  and 5. In these both cases, the maximum stress occurs before the steady state. It can be proved, in a theoretical/mathematical sense, that the point where the maximum breakdown rate occurs is before the stress peak. This fact is in agreement with expected behavior of thixotropic substance.
Other interesting result is related to the energy behavior of the thixotropic substance (Figures 6 and 7). Note the presence of a specific slope change in the region where higher rates of decrease on the micro-structures (λ μ , λ μ ) levels occurred. The same behavior can be observed on experimental results. In fact, this can be explained, in a theoretical sense, via the following relationship  On the Thermodynamic Consistency of a Two Micro-Structured Thixotropic Constitutive Model DOI: http://dx.doi.org/10.5772/intechopen.75987 that comes from the rate Eqs. (11) and (12).

Constant stress test
This chapter presents some points and results for the constant shear stress test. In this sense, it is considered loads as τ t ð Þ ¼ H t ð Þτ 0 , where H t ð Þ represents the standard Heaviside function and τ 0 is a positive real constant. It is assumed that in the begin of the test, the micro-structure is in a fully structured state (λ μ 0 ð Þ ¼ λ ν 0 ð Þ ¼ 1). It can be seen in Figure 8, for the loads τ 0 ¼ 10, 20, 30, 40, 50 and 100 Pa, a typical behavior for thixotropic substances under low level of constant stress loads. In this test, it can be seen the "Avalanche effect," and their relation with the microstructural level (λ μ , λ ν ) of the substance (Figures 9 and 10). In these cases, it can be observed the strain rate decrease, related to an increase of the construction parcel (cp μ ) and decrease of the destruction parcel (dp μ ) of λ μ , as can be seen in Figures 11  and 12 respectively, where For the Maxwell element (: ν ), in the same region (the strain rate decreases), cp ν is close to null value and dp ν is increasing, but this increase level is not sufficient to stop the decreasing process of _ λ. A posterior increase of _ γ is due to an abrupt decrease of linkages number (λ μ , λ ν ) before the steady state.
It is important to note the relationship between the structural nature and the behavior of the substance (Figures 9 and 10). Small changes in the values of λ μ and λ ν were detected for stress loads 10 and 20 Pa, in relation to the others. In these two cases _ γ presents an nonincreasing behavior. Other interesting result is presented in Figure 13, where it can be seen the energy behavior along the test. Note a changing on slope of the energy lines, tending to horizontal slope, in the _ γ decreasing region and another one related to the _ γ increasing region. Figure 11. cp μ vs: t: Figure 12. dp μ vs: t:

Concluding remarks
The main objective of this chapter is to investigate the constitutive model for thixotropic fluids based on [8] related to the existence of two different types of micro-structure, and their consistence in some rheological tests. The model and analysis presented here are based on well-established physics principles. In this sense, it is important to stand out some nonstandard points of the approach presented in this work for thixotropic modeling in respect to the others models. In this sense, it is important to stand out some points: • the shear modulus (G λ ν ð Þ), the viscosity coefficients (η μ λ μ À Á and η ν λ ν ð Þ) and their dependences from the two different micro-structure types are considered in the development of the constitutive model (Eqs. (9)-(14)); • the set of the rate equations (Eqs. (11)- (12)) are related to well-established physical principles as the "reptation" model and the "Smoluchowski" theory of coagulation. It is important to point out that in the "Smoluchowski" theory of coagulation the effect of the Brownian motion is clearly taken into account; • the thermodynamic consistence of the model was analyzed (Section 3); • a theoretical criterium for the transition region based on the strain gradient mapping degeneration was discussed and exploited (Section 4), however it is important to stand out the necessity of more discussions and analysis on the relationships (Eqs. (48) and (56)) taking into account some additional characteristics (as temperature, …) and their consequences. These characteristics shall be considered in future work goals; • the illustrative numerical examples attest the capability of the model to predict the expected behavior of real thixotropic substances under some typical rheological tests (Section 5). It is also important to comment that the developed ideas and presented in this work can be easily extended to approaches including more than two types of microstructure (related to each specific relaxation mechanism) that are taken into account in the constitutive equation system. It is clear that the complexity level of the considered substance determines the necessity of these incorporations. In this context, it can be noted the versatility of the approach exposed here, presenting some interesting perspectives on thixotropic modeling that can be more explored in future works. Oldroyd time derivative ν associated to the Maxwellian mechanism μ associated to the pure viscous mechanism first-order tensor second-order tensor Á inner product :

Symbology
double inner product Author details

Introduction
The increasing demand for a higher thermal efficiency and a higher operational flexibility of modern power generation plants results in severe mechanical and thermal loading. Modern energy markets put a requirement of a high operational flexibility of power plant units [1]. The flexible operation generates elevated loads and stresses in turbine components leading to material damage due to thermomechanical fatigue. In order to control the fatigue damage, thermal stresses should be computed online and limited to permissible values [2].
In power generation industry, the approach based on the Duhamel integral and the Green functions has been widely used for thermal stress calculations [3][4][5][6]. Such an approach has also been adopted for monitoring power boiler operation [7] and can also be used for calculating stress intensity factors at transient thermal loads [8].
The major issue in using Green's function and Duhamel's integral method in modeling transient thermal stresses in steam turbine components is the time dependence of material physical properties and heat transfer coefficients affecting proper evaluation of Green's functions and making the problem nonlinear. There are known approaches assuming the determination of the influence functions at constant values of these quantities [9]. The inclusion of temperature-dependent physical properties proposed by Koo et al. [10] relies on determining the weight functions for steady-state and transient operating conditions. A more important variation of the heat transfer coefficient can be taken into account by calculating the surface temperature using a reduced heat transfer model and employing Green's function to calculate a stress response to the step change of metal surface temperature [11]. However, numerical solution of a multi-dimensional heat transfer model for complex geometries present in steam turbines is complicated and timeconsuming, and due to this, it cannot be used in online calculations. A full inclusion of time variability of the physical properties and heat transfer coefficients proposed by Zhang et al. [12] relies on the solution of nonlinear heat conduction problem by using artificial parameter method and superposition rule and replacing the timedependent heat transfer coefficient with a constant value together with a modified fluid temperature. The effectiveness of the method has been proved by an example of a three-dimensional model of a pressure vessel of nuclear reactor.
This chapter presents and compares two alternative methods for steam turbine components both employing the Green functions and Duhamel's integral and considering the variations of material physical properties and heat transfer coefficient. The first method uses a numerical solution of one-dimensional (1D) heat conduction problem and determination of equivalent steam temperature for use with the Green function determined with a constant heat transfer coefficient. The second method employs the equivalent Green function determined with variable physical properties and heat transfer coefficient. The methods have been validated and compared for a steam turbine valve with the three-dimensional (3D) numerical solutions.

Thermoelasticity problem
The problem under consideration is heating up and cooling down of a steam turbine taking place during start-up and the subsequent shutdown. Heating and cooling of steam turbine components is caused by a hot steam flow through the turbine with varying temperature, pressure and velocity. Non-stationary heat transfer takes place via forced convection, and the thermal load is a primary load of the turbine. The rotor rotates with a constant rotational speed ω in a steady-state operation. Steam pressure and rotational body forces due to rotation are also considered in the model. Non-uniform temperature distribution in the component induces thermal stresses due to thermal expansion.
The material is assumed isotropic, and its physical and mechanical properties are temperature-dependent. The heating and cooling process is described by the linear theory of heat conduction and nonlinear convective boundary conditions. The resulting thermal stresses are determined from the solution of system of equations of uncoupled thermoelasticity. Due to the geometrical and loading symmetry of the turbine components like rotors or valve casings, the computational region is assumed axisymmetric. Mechanical loads are modeled as surface pressure and volumetric body force.
Based on the above assumptions, a boundary problem of heat conduction and thermoelasticity is formulated [13,14] and solved by means of a finite element method (FEM) [15].
Knowing the axisymmetric temperature field in the component, a solution of the problem of stress state induced by it can be obtained. The equilibrium conditions written for stresses are transformed taking into account the relations between stresses, strains and displacements, to obtain differential equations with unknown displacements u and w in the radial r-and axial z-direction, respectively [16], where G ij (r, t) is a Green function for thermal stress tensor component σ ij , ij = r,θ,z.
For simple geometries like a plate, cylinder or sphere, Green's functions for temperature and stress can be determined analytically by solving one-dimensional transient heat conduction problem. In case of more complicated shapes, it is necessary to use numerical methods, for example, finite element method in order to obtain accurate Green's functions.
In order to solve Duhamel's integral given by Eq. (11), it is transformed from a continuous into a discrete form: where t d is a cutoff time, G s , ij (r) is a value of the influence function in steady state, while G t , ij (r) is a transient part of the influence function.
Green's functions are determined individually for each supervised region, and in case of thermal stresses, these functions must be determined for each component of the stress tensor σ ij .
An example of Green's function for temperature and stress components at a step increase of fluid temperature by 1°C heating the external surface of a long cylinder is shown in Figure 1. The calculations were performed with a constant heat transfer coefficient α = 1000 W/m 2 /K and constant material physical properties evaluated at temperature 350°C. The temperature ( Figure 1a) and stress (Figure 1b) responses are presented for the cylinder axis and its outer surface heated by steam.
Green's function for surface temperature X_surf initially rapidly increases and reaches nearly 1 after approx. 10 4 s. The temperature at cylinder axis X_cent shows some inertia, and a delay time of temperature response equal to 500 s is seen. After this time, the temperature increases at approximately constant rate, and after c.a. 4000 s, the rate of variation visibly starts to decrease. The cylinder surface response is typical for surface-type sensor, while Green's function at cylinder axis is a response of middle-type sensor. The stress response at both areas is a result of temperature variations: axial and circumferential stresses at cylinder surface are equal to each other and rapidly grow reaching a minimum (compressive stress) and subsequently tend slowly to zero; the stresses at cylinder axis grow more slowly, their maxima are lower than at surface (axial stress is two times higher than circumferential stress) and are tensile stresses (sign +). All stress components within the cylinder tend to zero and a stress-free state is reached after approx. 1.4 � 10 4 s, which corresponds to the instant when temperature within the whole cylinder is uniform.

Green's function for variable physical properties and heat transfer coefficient
The above-presented Green functions were determined with constant physical properties and constant heat transfer coefficient. The physical properties of steel influencing temperature (specific heat c p , thermal conductivity λ) and stress (Young's modulus E, Poisson's number ν, thermal expansion coefficient β) distribution strongly depend on temperature, while the heat transfer coefficient varies with temperature and flow velocity (via the Reynolds number Re). Among the considered physical properties, the most varying is the specific heat which, for a typical low-alloy heat-resistant steel, increases by more than 50% in the temperature range of 20-500°C. The largest drop is observed for the Young modulus which decreases by more than 30% in the same temperature range. Heat transfer coefficient depends mostly on the flow velocity and during turbine start-up can change even by two orders of magnitude.
A consequence of the described variations of physical properties and heat transfer coefficient is a dependence of the Green function on these parameters. The dependence is illustrated in Figure 2 showing the Green function variations for axial stress at cylinder surface calculated by the finite element method at different temperatures and heat transfer coefficients. Each individual curve corresponds to stress variation determined at constant temperature T ¼ const and constant heat transfer coefficient α ¼ const. As it is seen from the plots, the largest effect on the Green function is observed for the heat transfer coefficient, which significantly affects the stress maxima, time of their occurrence and decay rate. The influence of temperature variation is much smaller and most visible from the instant of maximum stress to the moment of their vanishing.

Practical methods for thermoelastic stress calculation
In order to consider the variation of the Green function in Eq. (11), two methods are proposed: • Use of an equivalent steam temperature determined with a constant heat transfer coefficient α ¼ const ¼ α c • Use of an equivalent Green's function determined with variable heat transfer coefficient and physical properties

Equivalent steam temperature T eq
A flow chart of the proposed stress calculation algorithm using an equivalent steam temperature is presented in Figure 3 [21]. The calculation process is split into two steps: thermal analysis (step 1) and stress analysis (step 2). In step 1, thermal calculations are performed online, and one-dimensional (1D) heat conduction problem with time-dependent boundary conditions and temperature-dependent physical properties is solved using the finite difference method (FDM). Boundary conditions are given by the measured steam temperature and heat transfer coefficient computed on the basis of the steam mass flow rate. By solving the 1D heat conduction equation, one obtains the radial temperature distribution in the component model. The use of FDM enables full consideration of nonlinearities introduced by time-dependent boundary conditions and temperature-dependent physical properties. Knowing the thermal boundary conditions and radial temperature distribution T r; t ð Þ, the equivalent steam temperature T eq t ð Þ for the constant heat transfer coefficient α c is computed from the surface heat flux (see Eq. (22)). Calculations of the equivalent steam temperature are performed online at each time step for the same constant heat transfer coefficient α c .
In step 2, thermal stress calculations are performed online with the help of the Duhamel integral. Input for the integral is given by the equivalent steam temperature from step 1 and the Green functions determined offline by means of two/three dimensional (2D/3D) finite element method (FEM) calculations carried out for the real component geometry and constant heat transfer coefficient α c .
Heat conduction in a homogeneous isotropic solid is described by the Fourier-Kirchhoff differential equation [22] div λgradT where T is the metal temperature, g is the heat source and c p is the specific heat. Neglecting the heat source g and assuming heat conduction in the radial direction r only, Eq. (13) is re-written in a cylindrical coordinate system as follows: A uniform temperature distribution T r ð Þ ¼ T 0 is assumed as initial condition at t ¼ 0. Boundary conditions are defined depending on the component modeled. For rotor sections modeled as an externally heated solid cylinder, the boundary conditions are.

∂T ∂r
where T surf is the metal surface temperature.
For a hollow cylinder model, condition (16) directly applies to the outer surface, and for the inner surface, the boundary condition assumes the form.

∂T ∂r
In case of turbine and valve casings, their thermal behavior is approximated by an internally heated hollow cylinder with the following boundary conditions: ∂T ∂r The one-dimensional heat conduction Eq. (14) is numerically solved using the finite difference method [23]. The material physical properties are temperaturedependent, and their variation is described by polynomial functions. The variation of heat transfer coefficient α t ð Þ is considered by the Nusselt number Nu changing in time as a function of the Reynolds number Re and Prandt number Pr: The Nusselt number is defined as follows: where d denotes the cylinder characteristic diameter. A detailed form of Eq. (20) depends on the component geometry and flow character [24]. Practical methods of determining heat transfer coefficients are presented in [3]. The heat transfer model adopted in this study is based on the well-proven formulae for heat transfer coefficients in steam turbine components and can be employed in online calculations of temperatures and thermal stresses which are not possible when more advanced thermal FSI (fluid-structure interaction) modeling is adopted [25][26][27][28].
The convective surface heat flux defined by Eq. (16) or (18) with timedependent heat transfer coefficient α t ð Þ can be alternatively calculated using an equivalent steam temperature T eq and a constant heat transfer coefficient α c [21]: This relationship can be used to determine the equivalent steam temperature T eq which together with the assumed constant value of heat transfer coefficient α c ensures the same surface heat flux. The equivalent steam temperature evaluated as is transferred to step 2 ( Figure 3) for stress calculations. The algorithm basing on the equivalent steam temperature has a capability of modeling different heat transfer mechanisms. The conventional application of the Duhamel integral is based on predefined Green functions evaluated for a constant heat transfer coefficient assuming convection heat transfer mode, which is not the only heat exchange mechanism occurring in steam turbines operation. When the turbine starts from a cold state, steam condensation can take place on the rotor outer surfaces or casing inner surfaces, and condensation heat transfer should be taken into account. It is assumed that condensation occurs on the wall, if the following conditions are satisfied [21]: • The steam temperature is below the saturation temperature at the local pressure p • The heat flux from the water film to metal is greater than the heat flux from steam to the water film Condensation heat transfer takes place with very high heat transfer coefficients resulting in high heat fluxes on the component surface. This, in turn, brings about high surface temperature rates producing elevated thermal stresses on the surface or in the stress concentration areas. The most favorable conditions for the condensation to take place prevail at the initial phase of cold start-ups when the metal temperature is low, that is, below 100°C, and the steam saturation temperature high enough to produce high heat fluxes. This phenomenon is minimized by prewarming steam turbine components and supplying steam of possibly highest superheating.

Equivalent Green's function
The second method relies on using an equivalent Green's function in Eq. (11) determined, assuming a variable heat transfer coefficient and material physical properties. It can be done if the variation of process parameters is known in advance and used to determine the shape of the equivalent Green's function. This is the case in steam turbines operation as steam parameters are defined in design start-up diagrams for specific types of starts and should be followed within some margins in real operation. The method assumes little deviations of real operating conditions from those defined by the start-up diagrams.
The equivalent Green function is evaluated by performing finite element thermal analysis with the following boundary conditions [29]: 1. Steam temperature step characteristic for a given start-up type.
2. Variable heat transfer coefficient resulting from the variation in time of the steam mass flow rate characteristic for a given start-up.
The flow chart of the thermoelastic stress calculation algorithm based on the equivalent Green function is shown in Figure 4. The equivalent Green function is pre-computed offline by performing finite element calculations in 2D or 3D model of the real component. The finite element model is nonlinear due to the temperature-dependent material properties and time-dependent heat transfer coefficient. The equivalent Green functions are evaluated for each stress component and for various transient events like a cold, warm and hot start-up and shutdown. These events are characterized by specific variations of the heat transfer coefficients and the range of temperature variations affecting the material physical properties. The appropriate Green function is selected at the beginning of the transient event and used to calculate the thermal stress online with the help of the Duhamel integral.

Validation of the methods for turbine valve
First validations of both methods for a simple cylinder heated externally by steam were presented in [21] for the equivalent steam temperature and in [29] for the equivalent Green's function. A more complex example of a steam turbine valve operating at highly variable steam conditions, like temperature, pressure and flow rate, was used in this study to compare and validate the proposed algorithms. Numerical calculations were performed using a FEM model, and the presented algorithms with transient boundary conditions corresponding to a cold start-up are shown in Figure 5. The equivalent steam temperature was obtained for the calculated surface temperature and the assumed heat transfer coefficient α c ¼ 500 W/ m 2 K �1 . The equivalent Green functions were computed assuming the variation of heat transfer coefficient shown in Figure 5. As it is seen, the equivalent steam temperature exceeds the real steam temperature at a time point where the actual heat transfer coefficient rises above the assumed constant value of α c , which is a consequence of the equal surface heat fluxes assumed.  Figure 6 presents the variation of surface heat flux Q surf ¼ α T surf � T st À Á and surface temperature rate dT surf /dt during the cold start-up. Both curves exhibit two local minima: the first one after 2340 s and the second one after 3840 s. They correspond to the instants of the change in slope of the heat transfer coefficient curve shown in Figure 5 and will affect the shape of the Green functions and the variation of the equivalent Huber-Mises stress.
The calculations were performed with variable material physical properties depending on metal temperature. Figure 7 presents their variation during the cold start-up at the point of highest thermoelastic stresses located at the bowl inner surface. The largest drop is observed for the Young modulus which changes by around 30% during the cold start, while the largest rise occurs for the specific heat varying by more than 50% at the same time. Figures 8 and 9 show the temperature and Huber-Mises stress distributions in the valve casing at three instants of the cold start calculated with the FEM model. These stresses were used to validate the proposed algorithms by comparing them  with the predictions of both methods. The temperature field in the valve casing is two-dimensional, which results from both the non-regular geometry and nonuniform heat transfer conditions prevailing on the heated surfaces. The most intensive heat transfer takes place in the bowl and diffuser seat area, and the inner surface temperature is highest in these regions. The highest stress occurs at the inner surface of the valve bowl where both the temperature and stress gradients are predominantly present in the radial direction. Thermal stress distribution in the bowl wall ( Figure 9) is similar to the typical distribution in a cylindrical or a spherical wall with stress maxima occurring on the outer and inner surfaces and the   minimum stress present approximately in the middle of the wall thickness. The FEM model was also used for determining the equivalent Green functions and Green functions for a constant heat transfer coefficient α c = 500 W/m 2 /K in the bowl region and variable physical properties corresponding to the temperature range 30-525°C. The Green functions for each stress component at the bowl inner surface for the constant heat transfer coefficient are shown in Figure 10 and for the varying heat transfer coefficient are presented in Figure 11. All stresses are compressive, and except for the radial component, rapidly increase to their maxima and then slowly decay. The largest stresses are generated in axial and circumferential directions (G 22 and G 33 , respectively), perpendicular to the direction of the highest temperature gradient. These Green functions implicitly include also the effect of 2D temperature distribution resulting from different heat transfer conditions  prevailing on different valve surfaces heated by steam. This effect was accounted for by applying to these surfaces the heat transfer coefficient values different than α c = 500 W/m 2 /K (in the equivalent steam temperature method) or α = var (in the equivalent Green function method) which was achieved by applying the appropriate multiplication factors to the base heat transfer coefficients. There is a clear qualitative difference observed between the Green functions obtained with the constant and variable heat transfer coefficient. When the applied heat transfer coefficient is constant during the heating process, all the Green functions exhibit one maximum at the very beginning of the process, followed by continuous stress decay ( Figure 10). The effect of non-monotonic variation of the surface heat flux is not included here in the Green function but is taken into account in the equivalent steam temperature rate exhibiting two local minima ( Figure 6). When the applied heat transfer coefficient is varying during the heating process, all the Green functions exhibit two maxima: the first one at the beginning of the process and resulting from the initial temperature step and the second one corresponding to the increase of the surface heat flux (Figure 11).
The variation of Huber-Mises stress at the bowl obtained from the FEM analysis is compared in Figure 12 with the predictions of both simplified algorithms. Calculations were carried out using the Green functions shown in Figures 10 and 11, and the steam temperatures presented in Figure 5. As it is seen, both methods predict the occurrence of two stress maxima with high accuracy as compared with the FE results. A slightly better agreement is found for the method using the equivalent Green function which underpredicts the first maximum by 1.1% and the second maximum by 2.5%. The instants of both maxima correspond to those obtained from the FE model. The method using the equivalent steam temperature overpredicts both the first and second stress maximum with an error of 2.8 and 3.3%, respectively. A worse agreement with the FE results is also observed for the instants of stress maximum: the first maximum occurs later while the second one sooner as compared with the FE calculation results.
The occurrence of two subsequent stress maxima predicted in the finite element analysis is a result of variable surface heat flux exhibiting in time local maxima and minima. This effect was reproduced by both simplified methods but in a different way. In the equivalent Green function method, this function had two maxima and in combination with the applied steam temperature profile resulted in the two stress maxima. In the equivalent steam temperature method, this effect was obtained as a result of specific variation of the equivalent steam temperature applied together with the Green function of a typical shape.
Also, the stress decay period persisting after 4000 s of the cold start-up is accurately modeled by the two simplified methods. The relative error in this phase is generally below 10%, and initially, both methods underpredict the FE equivalent stress with a tendency to overpredict it at longer times.

Summary
The presented methods for online calculation of the thermoelastic stresses in steam turbine components are based on the Green function and Duhamel integral and consider the effect of variable heat transfer coefficient and material physical properties on thermal stresses. This effect is taken into account either by using the equivalent steam temperature determined with a constant heat transfer coefficient or by applying the equivalent Green's function determined with variable heat transfer coefficient and physical properties.
The first method is based on a numerical solution of 1D heat conduction problem and the Duhamel integral for thermal stress calculation. The equivalent steam temperature is calculated in step 1 and then used as an input to thermal stress model employing the Green function evaluated with a constant heat transfer coefficient. This approach allows for a full consideration of nonlinearities introduced by variable physical properties and heat transfer coefficient in temperature calculations. The effect of variable heat transfer coefficient significantly affecting the thermal stress results was taken into account by introducing the equivalent steam temperature evaluated on the basis of surface heat flux. The minor influence of variable physical properties on the stress response was considered by determining the Green functions from a step change of steam temperature in the range of its possible variation.
The second method relies on using the equivalent Green's function determined assuming a variable heat transfer coefficient and material physical properties. The equivalent Green function is evaluated by performing finite element thermal analysis for a step change of steam temperature and variable heat transfer coefficient resulting from the variation in time of the steam mass flow rate.
The effectiveness of both methods was shown by comparing their predictions with the results of finite element calculations of a steam turbine valve. The multidimensional geometry and non-uniform thermal boundary conditions result in twodimensional temperature and stress fields calculated using the FEM model. The comparison of equivalent stresses in the casing bowl calculated using the FEM model and the proposed methods revealed their good accuracy. The relative errors of maximum stress prediction are in the range of 1-3% and can be considered as very low taking into account the complex geometry of the valve, nonlinear boundary conditions and a broad range of temperature variation inducing significant variation of material physical properties. The methods can be applied for online calculation of thermoelastic stress utilized in thermal stress control and lifetime monitoring systems of steam turbines.