Characterization of Mechanical Heterogeneity in Dissimilar Metal Welded Joints

Dissimilar metal welded joints (DMWJs) possess significant localized mechanical heterogeneity. Using finite element software ABAQUS with the User-defined Material (UMAT) subroutine, this study proposed a constitutive equation that may be used to express the heterogeneous mechanical properties of the heat-affected and fusion zones at the interfaces in DMWJs. By eliminating sudden stress changes at the material interfaces, the proposed approach provides a more realistic and accurate characterization of the mechanical heterogeneity in the local regions of DMWJs than existing methods. As such, the proposed approach enables the structural integrity of DMWJs to be analyzed in greater detail.


Introduction
The leading pressure-bearing equipment in the primary water systems of pressurized water reactors (PWRs) is connected to the main pipeline using a safe-end dissimilar metal welded joint (DMWJ). Because of their high corrosion resistance and suitable mechanical properties, nickel-based alloys are widely used as welding metals for connecting stainless steel piping and pressure vessel nozzles [1,2]. In general, DMWJs consist of at least three base materials and weld materials and are subjected to a variety of complex loads while in service; therefore, structural failures often occur at DMWJs [1][2][3]. Given their status as a critical region for structural integrity analysis, the mechanical heterogeneity in the locality of DMWJs needs to be expressed appropriately. However, current approaches for assessing the structural integrity of DMWJs have certain limitations and rely heavily on the use of general material strength analysis and the operational experience of the inspector [4]. As such, a more accurate expression of the mechanical heterogeneity of DMWJs is required to assist in the analysis of their local mechanical properties [5].
The highly heterogeneous distribution of fracture properties, mechanical properties, and microstructure along the DMWJ has restricted the development of mechanical property evaluation systems [6,7]. Existing structural integrity assessments of welded joints, such as those proposed by Fan et al. [8,9] and Xue et al. [10], tend to ignore the influence of the heat-affected zone (HAZ) and the fusion zone (FZ). The simple bi-or tri-material models typically used do not evaluate the mechanical properties of strength-mismatched welded joints directly, potentially affecting the accuracy of calculations. Engineering applications that consider the mechanical heterogeneity of DMWJs usually use sandwich composite structures [11][12][13]. For simplicity, the welded joints are partitioned, with each partition corresponding to the mechanical parameters of a specific material. Despite their outstanding structural properties, such sandwich composite structures are mismatched in terms of their interregional material and geometric properties [14,15]. Significant interfacial stress variations, induced by mechanical loadings, occur at each regional interface [16,17]. Consequently, sandwich composite structures retain several limitations in terms of expressing the mechanical properties of the materials in the locality of the DMWJ [18,19]. To address the abrupt changes in the interfacial stresses encountered in simplified sandwich composite structures, a temperature field can be introduced to analyze the mechanical behavior in the locality of the DMWJ. To approximate the mechanical heterogeneity in terms of the temperature field, different temperatures can be defined at different locations [20,21]. By reproducing this method, we found that it yielded an approximation of the mechanical heterogeneity in the local regions of the DMWJ. To represent the mechanical heterogeneity at arbitrary positions in a simple and robust manner, the mechanical parameters were defined for the different element integration points [22].
This objectives of this study were: (i) to develop a more appropriate method for characterizing the mechanical heterogeneity of the DMWJ based on the User-defined Materials (UMAT) subroutine and (ii) to elucidate the variation of mechanical properties in the local regions of the DMWJ. The heat-affected and fusion zones of mechanical heterogeneity were considered in the analysis of the changes in the mechanical behavior of local regions of the DMWJ. Numerical simulations of the stress-strain curves at different positions in the vicinity of the DMWJ interface indicated that the interface location was likely to be a weak point in DMWJ failure. The proposed method can be applied to analyze the integrity of important engineering structures and to improve the stability of vulnerable structures.

UMAT Theoretical Basis
The UMAT subroutine in the finite element analysis (FEA) software ABAQUS (Abaqus 6.14, Dassault Systemes, Paris, France) can be used to define the mechanical behavior of materials and perform calculations at different element integration points [23]. In this study, the mechanical behavior of the element integration points was matched, pointby-point, to the coordinate positions in order to determine the mechanical properties at different positions in the DMWJ more realistically. Compared with previous methods, this approach was simpler and offered better accuracy and reproducibility. The UMAT subroutine calculated the trial and equivalent stresses, the equivalent plastic strain, and updated the state variables.
When an external force is applied, the subroutine calculates the elastic trial stress σ trial n+1 at the moment (n + 1) based on the total strain increment ∆ε applied at the moment t n ; thus, the total strain increment ∆ε is decomposed as where ∆ε e and ∆ε p are the elastic and plastic components of the total strain, respectively. Based on the assumption that the initial state is purely elastic, the trial stress at the beginning of the subroutine can be expressed as where trial represents the trial state and D e is the Jacobian matrix in the elastic state: where δ ij denotes the abbreviation for the unit matrix (when i = j, δ ij = 1; when i = j, δ ij = 0), E is Young's modulus, and v is Poisson's ratio. Next, the equivalent stress σ can be expressed as [24]  where S represents deviatoric stress. According to the J 2 flow rule [22], which means that the plastic strain ε p n+1 at the moment t n+1 is According to Ref. [24], the equivalent plastic strain increment at the nth step is Updating the stress values in the UMAT subroutine necessitates using different algorithms for the elastic and plastic stages. When the stress lies within the elastic range, the relationship between the stress and strain increments can be expressed as whereas when the stress is beyond the elastic range where D ep = D e − D p in which D p = 9G 2 * S * S T σ 2 (H +3G) [24,25].
The integration algorithm used by the UMAT subroutine comprises two main parts: the elastic trial and the return mapping. When loaded by an external force, the elastic trial stress σ trial n+1 at the moment t n+1 and the difference between it and the yield function f trial n+1 are calculated by the subroutine according to the total strain increment ∆ε applied at the moment t n , with the outcome determining whether the stress enters the plastic phase. When f trial n+1 ≤ 0, the integration algorithm is in the elastic phase a . In addition, because the plastic strain increment ∆λ = 0, the total strain increment ∆ε represents the elastic strain increment and the elastic stiffness matrix D e is used to update the stress at time t n to that at time t n+1 . Conversely, when f trial n+1 > 0, the integration algorithm is in the plastic phase b ; therefore, the plastic strain increment ∆λ > 0, the total strain increment ∆ε is decomposed into elastic strain increment ∆ε e and plastic strain increment ∆ε p , and the stress entering the plastic stage uses the elastic-plastic stiffness matrix D ep to update the stress at time t n to that at time t n+1 . Figure 1 shows the corresponding algorithm flow diagram.

Material Model
The DMWJs were modelled using Alloy52M as the filler metal and low-alloy steel (LAS) SA508 and austenitic stainless steel 316L as the base metal [26]. The chemical compositions of various materials are listed in Table 1 [26]. The DMWJ sample is shown in Figure 2. The weld metal consisted of buttering Alloy52Mb (52Mb) and weld Alloy52Mw (52Mw). The inner wall of the tube possessed an overlayer of austenitic stainless steel 304. This study focused mainly on the base metal (SA508 and 316L) and weld (Alloy52M) zones, where the mechanical heterogeneity was clearest. Table 2 lists the room-temperature mechanical properties of various DMWJ materials according to Refs. [20,27].

Material Model
The DMWJs were modelled using Alloy52M as the filler metal and low-alloy steel (LAS) SA508 and austenitic stainless steel 316L as the base metal [26]. The chemical compositions of various materials are listed in Table 1 [26]. The DMWJ sample is shown in Figure 2. The weld metal consisted of buttering Alloy52Mb (52Mb) and weld Alloy52Mw (52Mw). The inner wall of the tube possessed an overlayer of austenitic stainless steel 304. This study focused mainly on the base metal (SA508 and 316L) and weld (Alloy52M) zones, where the mechanical heterogeneity was clearest.   Table 2 lists the room-temperature mechanical properties of various DMWJ materials according to Refs. [20,27]. Table 2. Mechanical properties of various DMWJ materials at room temperature [20,27].

Material Model
The DMWJs were modelled using Alloy52M as the filler metal and low-alloy steel (LAS) SA508 and austenitic stainless steel 316L as the base metal [26]. The chemical compositions of various materials are listed in Table 1 [26]. The DMWJ sample is shown in Figure 2. The weld metal consisted of buttering Alloy52Mb (52Mb) and weld Alloy52Mw (52Mw). The inner wall of the tube possessed an overlayer of austenitic stainless steel 304. This study focused mainly on the base metal (SA508 and 316L) and weld (Alloy52M) zones, where the mechanical heterogeneity was clearest.   Table 2 lists the room-temperature mechanical properties of various DMWJ materials according to Refs. [20,27]. Table 2. Mechanical properties of various DMWJ materials at room temperature [20,27].  When σ ≤ σ y , the material deforms elastically, whereas when σ > σ y , the material deforms plastically; thus where σ is the stress, ε e is the elastic strain, σ y is the yield strength of the material, and H is the material hardening coefficient in the plastic deformation phase. The continuous transition model developed in this study considers the continuous variation of the mechanical properties of the materials along the DMWJ. Taking the axial direction as the y-direction, the yield strength σ y of the welded joint varies along the axial direction according to the function σ y (y); upon entering the plastic deformation stage, the hardening coefficient H changes according to the function H (y). The yield strength σ y at different positions and the plastic-phase hardening coefficient H are expressed as Substituting Equation (11) into Equation (10) yields the principal equation for the continuous variation of the mechanical properties along the DMWJ: In the UMAT subroutine, the mechanical heterogeneity of the DMWJ is characterized by defining E, v, σ y (y), and H (y). at different positions in the DMWJ to establish a finite element model and then analyzing it.
Following the approach described in [28], the relationships between hardness and the yield and ultimate strengths in the base metal, HAZ, and weld metal are described by Equation (13). In addition, Equation (14) describes the plastic-phase hardening coefficient: According to References. [3,11,27], the yield strengths of the DMWJ constituent materials (SA508, Alloy52M, and 316L) at room temperature were determined to be 426 MPa, 400 MPa, and 345 MPa, respectively. However, the local mechanical properties of the HAZ and the FZ during welding differed substantially from those of the base material (BM) and the weld metal (WM). Therefore, for the interface region, considering only the mechanical properties of the BM or WM may produce nonconservative (unsafe) or overly conservative results. By applying the empirical conversion equation to microhardness experimental data [28], the yield strength and hardening coefficient distributions at different positions along the DMWJ can be obtained, as shown in Figure 3.   Combined with the yield strength and hardening coefficient distributions along the DMWJ, this study used a polynomial function with a curve fitting correlation coefficient of R 2 ≥ 0.95 to fit the individual regions of the yield strength and hardening coefficient distributions obtained using the continuous transition material model. The continuous changes in the mechanical properties of the materials in the DMWJ interface region were considered in the constitutive equation and integrated into the UMAT subroutine: y ∈ [0, 6.0] 429.7 − 3.5y + 0.5y 2 y ∈ (6.0, 7.5] −11, 077 + 2913y − 184y 2 y ∈ (7.5, 8 For analyses using the sandwich structure, the model was partitioned, with different regions assigned different yield strengths as per Equation (17):

Geometric Model
The DMWJ consisted of a base metal (A508 and 316L) away from the interface and a weld metal (buttering Alloy52Mb and weld Alloy52Mw). The transition from the base metal to the weld zone inevitably caused lattice distortion. Therefore, a narrow (0.5-1 mm) fusion zone (FZ) occurred near the interface-the composition, microstructure, and hardness of which all changed significantly. In addition, the structure and properties of the HAZ (1-5 mm) were unevenly distributed, owing to the heat-related effects of welding [3,11,17]. Figure 4 shows a simplified schematic diagram of the structure used for the FEA. It is important to note that material interfaces will not be perfectly straight in reality; however, for simplicity, straight interfaces were used in all simulations. fusion zone (FZ) occurred near the interface-the composition, microstructure, and hardness of which all changed significantly. In addition, the structure and properties of the HAZ (1-5 mm) were unevenly distributed, owing to the heat-related effects of welding [3,11,17]. Figure 4 shows a simplified schematic diagram of the structure used for the FEA.
It is important to note that material interfaces will not be perfectly straight in reality; however, for simplicity, straight interfaces were used in all simulations.

Mesh and Load Model
Although welding is a three-dimensional process, it is widely acknowledged that axisymmetric models are appropriate for simulating the welding of cylindrical structures to

Mesh and Load Model
Although welding is a three-dimensional process, it is widely acknowledged that axisymmetric models are appropriate for simulating the welding of cylindrical structures to determine the local mechanical properties of DMWJs [29]. Therefore, axisymmetric finite element models, which are fast and easy to use, are widely used to simulate heterogeneity. As shown in Figure 5, the mesh was composed of 4-node bilinear axisymmetric elements (CAX4) with a total of 1080 elements (0.5 mm × 0.5 mm). Displacement and rotation constraints were applied to the left end of the model in the y-direction and around the x-direction, respectively. The mechanical properties of the model were compared and analyzed by applying U y = 0.5 mm, pressure = 0 MPa, and U y = 0.5 mm, pressure = 100 MPa, respectively.

Mesh and Load Model
Although welding is a three-dimensional process, it is widely acknowledged that axisymmetric models are appropriate for simulating the welding of cylindrical structures to determine the local mechanical properties of DMWJs [29]. Therefore, axisymmetric finite element models, which are fast and easy to use, are widely used to simulate heterogeneity.
As shown in Figure 5, the mesh was composed of 4-node bilinear axisymmetric elements (CAX4) with a total of 1080 elements (0.5 mm × 0.5 mm). Displacement and rotation constraints were applied to the left end of the model in the y-direction and around the xdirection, respectively. The mechanical properties of the model were compared and analyzed by applying = 0.5 mm, pressure = 0 MPa, and = 0.5 mm, pressure = 100 MPa, respectively.

Continuous Transition Model Simulation Results
The numerical simulation results for the model subjected to axial loading are shown in Figure 6. Specifically, Figure 6a,b show the stress results for the developed continuous transition and sandwich models, respectively. The results obtained for applying the continuous transition model revealed that the stress in the interface region gradually increased at first before gradually decreasing. In contrast, for the sandwich composite structure model, the stress at the interface changed suddenly. As shown in Figure 7, the models were subjected to axial loading and internal pressure simultaneously. As a consequence, the

Continuous Transition Model Simulation Results
The numerical simulation results for the model subjected to axial loading are shown in Figure 6. Specifically, Figure 6a,b show the stress results for the developed continuous transition and sandwich models, respectively. The results obtained for applying the continuous transition model revealed that the stress in the interface region gradually increased at first before gradually decreasing. In contrast, for the sandwich composite structure model, the stress at the interface changed suddenly. As shown in Figure 7, the models were subjected to axial loading and internal pressure simultaneously. As a consequence, the internal stresses of the model changed, although they remained most significant at the interface region. This highlights the importance of the interface region when assessing the structural integrity of DMWJs. The next section discusses the stressstrain relationship in the interface region for a representative example in which an axial loading of 0.5 mm was applied. internal stresses of the model changed, although they remained most significant at the interface region. This highlights the importance of the interface region when assessing the structural integrity of DMWJs. The next section discusses the stress-strain relationship in the interface region for a representative example in which an axial loading of 0.5 mm was applied.    Figure 8a shows the stress-strain curves in the vicinity of the SA508/52Mb interface in response to an axial loading of 0.5 mm. The curves indicate the presence of high stress and low strain in the HAZ of the SA508 material. In contrast, the 52Mb region exhibited low stress and high strain. Figure 8b shows the stress-strain curves in the vicinity of the 52Mw/316L interface: compared with the 52Mw region, the 316L region, which was situated farther from the interface, exhibited lower stress and higher plastic strain. Nevertheless, both interface regions showed high stress and plastic strain, which promotes crack initiation and growth. Therefore, from the numerical simulation results, we inferred that the interface represents a high-probability zone for DMWJ failure, which is consistent with experimental results reported in the literature [27,30]. Figure 8 shows that the maximum stress at position e in the SA508/52Mb interface zone was 485 MPa, while the stress at point a in the base material was 435 MPa, representing a difference of 50 MPa. The maximum stress at position d in the 52Mb/316L interface region was 450 MPa, while the maximum stress at point a in the 52Mw region was 415 MPa, representing a stress change of 35 MPa. Owing to the larger magnitude of the stress change in the SA508/52Mb interface region, it was deemed more susceptible to destruction. As such, it should be focused on.

Local Stress-Strain Simulation
Results at the SA508/52Mb and 52Mw/316L Interfaces Figure 8a shows the stress-strain curves in the vicinity of the SA508/52Mb interface in response to an axial loading of 0.5 mm. The curves indicate the presence of high stress and low strain in the HAZ of the SA508 material. In contrast, the 52Mb region exhibited low stress and high strain. Figure 8b shows the stress-strain curves in the vicinity of the 52Mw/316L interface: compared with the 52Mw region, the 316L region, which was situated farther from the interface, exhibited lower stress and higher plastic strain. Nevertheless, both interface regions showed high stress and plastic strain, which promotes crack initiation and growth. Therefore, from the numerical simulation results, we inferred that the interface represents a high-probability zone for DMWJ failure, which is consistent with experimental results reported in the literature [27,30]. Figure 8 shows that the maximum stress at position e in the SA508/52Mb interface zone was 485 MPa, while the stress at point a in the base material was 435 MPa, representing a difference of 50 MPa. The maximum stress at position d in the 52Mb/316L interface region was 450 MPa, while the maximum stress at point a in the 52Mw region was 415 MPa, representing a stress change of 35 MPa. Owing to the larger magnitude of the stress change in the SA508/52Mb interface region, it was deemed more susceptible to destruction. As such, it should be focused on.

Effect of Mechanical Heterogeneity on the DMWJ Local Region Mechanical Properties
From Figure 9, it can be seen that the stress component in the x-direction was compressed in both the sandwich structure model and in the continuous transition model-a phenomenon most likely due to the axial loading-induced displacement. Moreover, this phenomenon was amplified by the effects of internal pressure. It was evident that the stress in the interface region changed significantly.

Effect of Mechanical Heterogeneity on the DMWJ Local Region Mechanical Properties
From Figure 9, it can be seen that the stress component in the x-direction was compressed in both the sandwich structure model and in the continuous transition model-a phenomenon most likely due to the axial loading-induced displacement. Moreover, this Materials 2021, 14, 4145 9 of 13 phenomenon was amplified by the effects of internal pressure. It was evident that the stress in the interface region changed significantly.

Effect of Mechanical Heterogeneity on the DMWJ Local Region Mechanical Properties
From Figure 9, it can be seen that the stress component in the x-direction was compressed in both the sandwich structure model and in the continuous transition model-a phenomenon most likely due to the axial loading-induced displacement. Moreover, this phenomenon was amplified by the effects of internal pressure. It was evident that the stress in the interface region changed significantly. As shown in Figure 10, when axial loads alone were considered, the y-directional stress component was significantly larger than when both the axial and the internal pipe pressure loads were considered. This is due to the reduction of the y-directional stress component by the internal pressure. The y-directional component exhibited the same behavior under both loading conditions, with the stress in the interface region changing abruptly. This further supports the need to prioritize the mechanical properties of the interface region. As shown in Figure 10, when axial loads alone were considered, the y-directional stress component was significantly larger than when both the axial and the internal pipe pressure loads were considered. This is due to the reduction of the y-directional stress component by the internal pressure. The y-directional component exhibited the same behavior under both loading conditions, with the stress in the interface region changing abruptly. This further supports the need to prioritize the mechanical properties of the interface region. As shown in Figure 11, under an axial load, the sandwich structure model showed a significant increase in stress at the SA508/52M and 52M/316L interfaces. For an internal pressure loading, the opposite trend was observed and a sudden change in stress occurred at the interface region. The continuous transition model eliminated abrupt changes in stress, allowing for a continuous transition in the mechanical properties of the materials in the interface regions.  As shown in Figure 11, under an axial load, the sandwich structure model showed a significant increase in stress at the SA508/52M and 52M/316L interfaces. For an internal pressure loading, the opposite trend was observed and a sudden change in stress occurred at the interface region. The continuous transition model eliminated abrupt changes in stress, allowing for a continuous transition in the mechanical properties of the materials in the interface regions.
As shown in Figures 12-15, by combining the analysis results for the stress components in the x-, y-, and z-directions, it was concluded that the sandwich structure model exhibited an abrupt change in stress and in the equivalent plastic strain (PEEQ) in the interface region. This can be attributed to the fact that the model assigns different mechanical properties for different partitions, while ignoring the existence of the HAZ and FZ. As shown in Figure 12, the results upon applying the developed continuous transition model revealed that, in the SA508/52M interface region, the stress gradually increased to 485 MPa before gradually decreasing to 405 MPa, whereas in the 52Mw/316L interface region, the stress gradually increased to 450 MPa before gradually decreasing to 376 MPa. In the sandwich structure model, the stress values did not transition continuously in the interface region but instead changed abruptly. This was evident in the sudden stress changes (from 430 MPa to 404 MPa in the SA508/52Mb interface region and from 404 MPa to 377 MPa in the 52Mw/316L interface region). For an internal pressure of 100 MPa, the analysis method was the same. In addition, the method proposed in this study can characterize the variation of mechanical properties between overlay layers composed of 52Mb and 52Mw. Thus, the method proposed herein is more appropriate for expressing the mechanical properties in the interface region and better for assessing the structural integrity of DMWJs. As shown in Figure 11, under an axial load, the sandwich structure model showed a significant increase in stress at the SA508/52M and 52M/316L interfaces. For an internal pressure loading, the opposite trend was observed and a sudden change in stress occurred at the interface region. The continuous transition model eliminated abrupt changes in stress, allowing for a continuous transition in the mechanical properties of the materials in the interface regions. As shown in Figures 12-15, by combining the analysis results for the stress components in the x-, y-, and z-directions, it was concluded that the sandwich structure model exhibited an abrupt change in stress and in the equivalent plastic strain (PEEQ) in the interface region. This can be attributed to the fact that the model assigns different mechanical properties for different partitions, while ignoring the existence of the HAZ and FZ. As shown in Figure 12, the results upon applying the developed continuous transition model revealed that, in the SA508/52M interface region, the stress gradually increased to 485 MPa before gradually decreasing to 405 MPa, whereas in the 52Mw/316L interface region, the stress gradually increased to 450 MPa before gradually decreasing to 376 MPa. In the sandwich structure model, the stress values did not transition continuously in the interface region but instead changed abruptly. This was evident in the sudden stress changes (from 430 MPa to 404 MPa in the SA508/52Mb interface region and from 404 MPa to 377 MPa in the 52Mw/316L interface region). For an internal pressure of 100 MPa, the analysis method was the same. In addition, the method proposed in this study can characterize the variation of mechanical properties between overlay layers composed of 52Mb and 52Mw. Thus, the method proposed herein is more appropriate for expressing the mechanical properties in the interface region and better for assessing the structural integrity of DMWJs. The mechanical property variations shown in Figures 12-15 are in general agreement with the room temperature stress-strain results obtained for different regions of welded components in previous studies [3,7,11]. Sudden changes in stress in the interface regions caused subsequent changes in the PEEQ. Therefore, DMWJs possess significant mechanical heterogeneity in multi-material interface regions, and the use of conventional sandwich structure models to study DMWJs may produce nonconservative (unsafe) or overly conservative results. Consequently, the mechanical property distributions for DMWJ interface regions must be treated properly when assessing the structural integrity of DMWJs.        [3,7,11]. Sudden changes in stress in the interface regions caused subsequent changes in the PEEQ. Therefore, DMWJs possess significant mechanical heterogeneity in multi-material interface regions, and the use of conventional sandwich structure models to study DMWJs may produce nonconservative (unsafe) or overly conservative results. Consequently, the mechanical property distributions for DMWJ interface regions must be treated properly when assessing the structural integrity of DMWJs.

Conclusions
This paper presented an approach for characterizing the mechanical heterogeneity in the local regions of DMWJs. Variations of mechanical properties in the local regions of the DMWJ composed of Alloy52M were analyzed, leading to the following conclusions: (1) The use of a user-defined material utility, UMAT, enables the continuous variation of mechanical properties in all regions of the DMWJ to be characterized to obtain more accurate and realistic results.
(2) By considering the local mechanical heterogeneity, the proposed method avoids mismatching interregional material and geometric properties and eliminates abrupt changes in interfacial stresses. Accurate finite element simulations were performed using ABAQUS commercial software. (3) Owing to their serious mechanical heterogeneities, particularly in the interface regions, the local mechanical properties of DMWJs change considerably. Compared with the 52Mw/316L interface region, the SA508/52Mb interface region is more susceptible to mechanical heterogeneity and should be regarded as a high-risk region in terms of structural integrity.
Finally, the proposed method can also be used to analyze the micromechanical fields and growth paths of cracks at different positions in mechanically heterogeneous DMWJ interface zones. Such analyses will be the subject of future research.