Partitioned coupling framework to combine a kinematic hardening plasticity model and a creep model for structures in a high‐temperature environment

The present study proposes a partitioned coupling framework for the implementation of a combined kinematic hardening plasticity and creep model. In the present framework, plastic and creep implicit stress integrations are performed alternately and iteratively until convergence in a partitioned coupling manner. This approach enables us to tractably implement a constitutive model that combines a plasticity model and a creep model in the program. The present framework was applied to the combination of the Ohno–Wang plasticity model and the Norton–Bailey creep model for numerical examples. The numerical examples demonstrated the capability of the present framework in analyzing realistic structures at high temperature with the consideration of both plasticity and creep.

temperature using the Ohno-Wang kinematic hardening model. 3,4 Abdel-Karim and Ohno 5 proposed a kinematic hardening model to model the steady state in ratcheting of the modified 9Cr-1Mo steel at 550 • C and IN738LC alloy at 850 • C. This kinematic hardening model is in the framework of strain hardening and dynamic recovery, which is similar to the Ohno-Wang model. Ohguchi et al. 6 and Sasaki and Ohguchi 7 analyzed soldering problems by a combined model of the Prager kinematic hardening model and the Norton creep model. These studies also provided experimental results for solder alloys that exhibited both plastic and creep behaviors. Kobayashi et al. 8 and Akamatsu et al. 9 also analyzed soldering problems by a nonunified constitutive model, in which both transient and steady state strain rates were considered. Using this model, Nishimoto et al. 10,11 simulated the combustion chamber of a liquid rocket engine in a high-temperature environment. Chen et al. 12 added the consideration of viscoplasticity and creep in the model of Kobayashi et al. 8 to model a nickel-based super alloy and 40Sn-60Pb solder at high temperature. Panteghini and Genna 13 used the model of Kobayashi et al. 8 in conjunction with the Garofalo hyperbolic-sine creep law and the viscoplasticity model of Perzyna 14 to calculate the residual stress of wire drawing processes. Hosseini et al. 15 analyzed a compact tension specimen of the 1%CrMoV steel at 550 • C using five different creep models, including the Norton, two-regime Norton, Norton-Bailey, Bartsch, and hyperbolic-sine models. The studies mentioned in this paragraph indicate that a combined plasticity-creep model is necessary in some situations rather than a single plasticity or creep model.
The most classical and simplest kinematic hardening model would be the Prager model, in which the back stress rate is proportional to the plastic strain rate. Frederick and Armstrong 16 added a dynamic recovery term to the evolution equation of the Prager model. Then, Chaboche 17 introduced the approach of multiple back stresses by superposing the Armstrong-Frederick model. After that, Ohno and Wang 3,4 proposed an elastic-plastic constitutive model with a kinematic hardening model in which strain hardening and dynamic recovery with a critical state are considered. The Ohno-Wang model has been used for high-temperature mechanical problems, such as ratcheting at high temperature, 1,2,5 soldering, 8,9 a liquid rocket engine, 10,11 and wire drawing. 13 Creep is frequently modeled by an evolution equation of the equivalent creep strain. The Norton creep model, which is also known as the power law creep model, is commonly used to represent secondary creep. The Norton creep model is a stress hardening model in which the equivalent creep strain rate is proportional to the equivalent stress to the power of a material parameter. Furthermore, for primary creep, the Bailey creep model, which is a time hardening model, is popular. Therefore, the Norton-Bailey model is appropriate for modeling both the primary and secondary creep. In the Norton-Bailey model, the equivalent creep strain rate is proportional to the equivalent stress to the power of a parameter multiplied by the time to the power of another parameter. The Norton or Norton-Bailey model has been applied to high-temperature problems such as soldering, 6,7 a nickel-based super alloy and 40Sn-60Pb solder at high temperature, 12 and 1%CrMoV steel at 550 • C. 15 Furthermore, for detailed modeling of transitions of mechanism depending on temperature and stress, a creep deformation mechanism map, such as the Ashby map, 18 is considered. For creep crack growth, a creep damage model, such as the Nikbin-Smith-Webster (NSW) model, 19 is used.
For the combination of plasticity and creep, several constitutive models have been proposed. These models can be categorized into unified and nonunified models. In unified models, 20 plastic and creep strains are considered to arise from a single physical mechanism. However, unified models have difficulty in determining the material parameters. In contrast, plastic and creep material parameters can be measured independently in nonunified models, [6][7][8][9]12,13 where the plastic and creep strains are assumed to arise from two independent constitutive laws. Nonunified models have been applied to problems at high temperature, such as solder joints, 6-9 a liquid rocket engine, 10,11 and wire drawing. 13 In the implementation of a nonunified model, a nonlinear system of equations consisting of plastic and creep stress integrations is solved monolithically for each element integration point. 8,9 This procedure requires an algorithm and a program that are specific to the combined plasticity-creep model used in the analysis.
To overcome this issue, the present study proposes a framework for the implementation of nonunified models. The present framework is based on a partitioned coupling technique, which can be seen in many fluid-structure interaction analyses. [21][22][23][24][25] Two existing plastic and creep models are combined in the framework of coupling iteration. This approach enables us to tractably reuse existing programs of a single plasticity model and a single creep model. The present framework is the generalization of the authors' previous study presented at a conference proceeding. 26 In the present study, numerical examples using a combined Ohno-Wang and Norton-Bailey model are presented. These numerical examples demonstrate the capability of the proposed partitioned coupling technique.
In the present article, a general description of the constitutive equations of the plasticity model and the creep model is first introduced. Then, we present partitioned stress integration algorithms and their consistent tangent moduli of the single plasticity model, the single creep model, and the combined plasticity-creep model. After that, the algorithm and the implementation of the partitioned coupling framework with the combined plasticity-creep model are explained. In particular, the advantage of the partitioned coupling framework in separate implementation of the plasticity and creep models is emphasized. Using the partitioned coupling framework, we introduce specific equations of a combined plasticity-creep model, namely, the combined Ohno-Wang and Norton-Bailey model. Finally, the capability of the proposed partitioned coupling technique with the combined Ohno-Wang and Norton-Bailey model is demonstrated by numerical examples.

CONSTITUTIVE MODELS
In this section, constitutive models, such as a general kinematic hardening plasticity model, a general creep model, and a combined plasticity-creep model, are introduced. First, vector representation of stress and strain tensors is explained, followed by operators to describe constitutive models. Then, a single kinematic hardening plasticity model and a single creep model are presented. Finally, a combined kinematic hardening plasticity and creep model is explained. Although a number of assumptions, such as the von Mises criterion and the associated flow rule, are introduced in the present study, the evolution equations of the back stresses and the equivalent creep strain are described in general forms.

Vector representation and operators
In the present study, vector representation of stress and strain tensors is used to describe constitutive models. A stress, such as the stress, , the deviatoric stress, s, the back stress, , or the relative stress, kin , is represented by a vector having six components, that is, In this representation, symmetry is assumed, for example, xy = yx . A strain such as the strain, , the elastic strain, e , the plastic strain, p , or the creep strain, cr , is represented by Shear components of a strain vector are doubled for convenience. In order to double the shear components of a stress vector, a linear operator, is introduced. Using Q, a norm operator, is defined.
A volumetric operator, T vol , and a deviatoric operator, T dev , are introduced: Note that T vol + T dev = I.

General kinematic hardening plasticity model
First, the strain, , is introduced. The strain is decomposed additively into elastic and plastic parts, that is, Here, e and p denote elastic and plastic strains, respectively. Then, the constitutive equation of a plasticity model is where denotes the stress. Here, D e is the elasticity matrix, which can be represented by where K and G are the bulk modulus and the shear modulus, respectively. In the present study, the plastic strain, p , is assumed to evolve by the associated flow rule, aṡ wherėis the plastic multiplier, and kin denotes the relative stress, which is the stress, , minus the back stress, , that is, In addition, g is the yield function, for which the von Mises yield criterion, is adopted. Here, kin is the von Mises equivalent stress with respect to the deviatoric part of the relative stress, kin , that is, where s is the deviatoric part of the stress, that is, s = T dev , and is the back stress. In Equation (11), Y is the yield stress, which is a function of the equivalent plastic strain, p , that is, Here, the equivalent plastic strain, p , is defined as Using Equations (9), (11), and (15), we havėp Based on the generalized Armstrong-Frederick model by Chaboche, 17 the back stress, , can be decomposed into M parts, that is, For the evolution of each part of the back stress, i , the contribution of the plastic strain rate,̇p, and that of the back stress component, i , are considered:̇i where h i is a material parameter, anḋi is a rate-type scalar variable, the expression of which is determined by the selection of a kinematic hardening model. The terṁi may depend on both material parameters and state variables. A nonlinear hardening can be represented by replacing the constant, h i , with a function of p .

General creep model
Similar to the previous subsection, the strain, , is decomposed additively into elastic and creep parts, that is, where e and cr denote elastic and creep strains, respectively. Then, the constitutive equation of a creep model is where and D e denote the stress and the elasticity matrix, respectively. The creep strain, cr , is assumed to evolve in a similar manner to the plastic strain (Equation 16), aṡ where the von Mises equivalent stress, , and the equivalent creep strain rate,̇c r , are Note that the creep strain evolves in the direction of the deviatoric stress, s, whereas the plastic strain evolves in the direction of the relative deviatoric stress, s − . Then, the equivalent creep strain rate,̇c r , in Equation (21) is described in a general form, aṡc

Combined kinematic hardening plasticity and creep model
Similar to Sections 2.2 and 2.3, the strain, , is decomposed additively into elastic, plastic, and creep parts, that is, where e , p , and cr denote elastic, plastic, and creep strains, respectively. Then, the constitutive equation is where and D e denote the stress and the elasticity matrix, respectively. The evolution equations of the plastic strain, p , and the creep strain, cr , are the same as the strains presented in Sections 2.2 and 2.3, respectively.

IMPLICIT STRESS INTEGRATION
In this section, stress integration procedures in the time interval from t to t + Δt are derived. First, stress integration of a single kinematic hardening plasticity model is presented, followed by that of a single creep model. Finally, partitioned stress integration of a combined kinematic hardening plasticity and creep model is explained.

Implicit stress integration of a single kinematic hardening plasticity model
The stress at t + Δt, t+Δt , can be decomposed into volumetric and deviatoric parts: The inelastic deformation is assumed to be incompressible. In this case, the volumetric part (hydrostatic stress), t+Δt p, is purely elastic, that is, where Δ is the strain increment. On the other hand, the deviatoric part, t+Δt s, is calculated simultaneously with the unknown plastic strain increment, Δ p , which is also related to t+Δt s and t+Δt . The backward Euler integration method combined with the elastic predictor and radial return method 27 is used for the stress integration to calculate stress and strain increments. The deviatoric stress, t+Δt s, can be calculated as where t+Δt s * is an elastic predictor that is calculated by If t+Δt s * is placed inside the yield surface (g < 0), then t+Δt s * is in an elastic state and t+Δt s = t+Δt s * with Δ p = 0. If t+Δt s * is outside the yield surface (g ≥ 0), then This equation was derived from Equation (16) using the backward Euler integration method. Here, Δ p is the equivalent plastic increment. In Equation (31), t+Δt and Δ p should be determined implicitly.
Integrating Equation (18) yields the ith back stress increment: From this equation, t+Δt in Equation (31) can be determined by For determining Δ p in Equation (31), a scalar equation will be derived. From Equations (29), (31), and (33)-(35), we Then, moving t+Δt s − t+Δt to the left-hand side yields Finally, Equation (38) is applied to Equation (12), resulting in The solution of this nonlinear scalar equation is Δ p .

Implicit stress integration of a single creep model
The stress can be decomposed into the hydrostatic and deviatoric stresses in the manner described in the previous subsection (Equation 27). Moreover, the hydrostatic stress follows Equation (28). The deviatoric stress, t+Δt s, can be calculated by where t+Δt s * is an elastic predictor that is given by Equation (30), and Δ cr is the creep strain increment, which is Here, Δ cr is the equivalent creep strain increment, which should be determined implicitly. Equation (41) Then, moving t+Δt s to the left-hand side yields Finally, Equation (43) is applied to Equation (22), resulting in The solution of this nonlinear scalar equation is . The equivalent creep strain increment, Δ cr , is the integration of Equation (24), that is,

Partitioned implicit stress integration of a combined kinematic hardening plasticity and creep model
The stress can be decomposed into the hydrostatic and deviatoric stresses by Equation (27). Moreover, the hydrostatic stress follows Equation (28). The deviatoric stress, t+Δt s, can be calculated as where t+Δt s * is an elastic predictor that is given by Equation (30). The plastic strain increment, Δ p , and the creep strain increment, Δ cr , can be computed using Equations (31) and (41), respectively. In Equations (31) and (41), t+Δt , Δ p , and Δ cr should be determined implicitly. Here, t+Δt follows Equations (33) and (34).
For determining Δ p in Equation (31) Then, moving t+Δt s − t+Δt in the former equation and t+Δt s in the latter equation to the left-hand sides yields where Finally, Equations (49) and (50) are applied to Equations (12) and (22), respectively, resulting in Although the nonlinear system of Equations (53) and (54) can be solved monolithically by the high-dimensional Newton-Raphson method, one cannot tractably reuse the stress integration algorithm of a single constitutive model. Moreover, when not yielding, the combined constitutive model becomes a single creep model. In this case, solving Equation (54) alone would be the simplest way. Therefore, we adopted a partitioned coupling approach in which Equations (53) and (54) are solved in a partitioned manner under an iterative method, such as the successive substitution method. First, Equation (54) is solved for t+Δt using the Newton-Raphson method. The convergence at the jth iteration step is checked by where is the tolerance. Then, after calculating several variables from the solution of Equation (54), the yield function, g in Equation (11), is evaluated. If t+Δt g ≥ 0, then Equation (53) is solved for Δ p using the Newton-Raphson method. The convergence is checked by These procedures are repeated until both Equations (55) and (56) are satisfied.

CONSISTENT TANGENT MODULI
In this section, tangent moduli that are consistent with the stress integration procedures in the previous section are derived. First, a consistent tangent modulus of a single kinematic hardening plasticity model is presented, followed by that of a single creep model. Finally, a consistent tangent modulus of a combined kinematic hardening plasticity and creep model is explained.

Consistent tangent modulus of a single kinematic hardening plasticity model
The consistent tangent modulus, , of a single kinematic hardening plasticity model is derived in this section. Here, d t+Δt can be decomposed into volumetric and deviatoric parts: The volumetric part is assumed to be purely elastic, that is, For the deviatoric part, differentiating Equation (29) yields Then, differentiation of Equation (31) gives where With t+Δt n kin ⋅ Q d t+Δt n kin = 0 that is derived from t+Δt n kin ⋅ Q t+Δt n kin = 1, multiplying both sides of Equation (60) by t+Δt n kin yields In addition, differentiation of Equation (61) leads to where Substituting Equations (59) and (62) into Equation (63), we obtain Then, replacing d t+Δt n kin in Equation (60) with Equation (65), we have where I is a 6 × 6 unit matrix. Here, t+Δt H i is the constitutive matrix obtained by differentiating t+Δt i with respect to Δ p , that is, t+Δt The detailed derivation of t+Δt H i for the Ohno-Wang model, as an example, is given by Kobayashi and Ohno. 28 Equation (66) establishes the relation between dΔ p and dΔ , which can be written as where Solving Equation (67) yields By adding the volumetric and deviatoric parts, the consistent tangent modulus of a single kinematic hardening plasticity model becomes Algorithm 1. Nonlinear finite element analysis 1: Initialize the nodal displacement: 0 u ← 0 2: for t ← 0 to t max do 3: Initialize the nodal displacement increment: Δu ← 0 4: Generate the external force, t+Δt f ext , from boundary conditions 5: repeat 6: Calculate the strain increment: Δ ← BΔu, where B is the strain-nodal-displacement matrix 7: Calculate the stress, t+Δt , by the implicit stress integration (Section 3) 8: Calculate the internal force: Calculate the consistent tangent modulus: 10: Calculate the coefficient matrix: t+Δt K ← ∫ V B Tt+Δt DBdV 11: Update the displacement increment by solving a linear system of equations: until convergence 13: Update the nodal displacement: t+Δt u ← t u + Δu 14: Update the time: t ← t + Δt 15: end for

Consistent tangent modulus of a single creep model
Differentiating Equation (40) yields Then, differentiation of Equation (41) gives where With t+Δt n ⋅ Q d t+Δt n = 0 that is derived from t+Δt n ⋅ Q t+Δt n = 1, multiplying both sides of Equation (72) by t+Δt n gives In addition, differentiation of Equation (73) leads to Solving Equation (77) yields Finally, the consistent tangent modulus of a single creep model becomes

Consistent tangent modulus of a combined kinematic hardening plasticity and creep model
Differentiating Equation (46) yields Substituting Equations (81) and (62) where I is a 6 × 6 unit matrix. Following a similar derivation, the following creep-related equations can be obtained: Equations (83) and (84) establish the relationship among dΔ p , dΔ cr , and dΔ , which can be written as where dΔ p + dΔ cr can explicitly be expressed in terms of dΔ by solving Equation (85), that is, Finally, the consistent tangent modulus of a combined kinematic hardening plasticity and creep model becomes

ALGORITHM AND IMPLEMENTATION OF PARTITIONED IMPLICIT STRESS INTEGRATION AND CONSISTENT TANGENT MODULUS
First, the algorithm of a nonlinear finite element analysis that is based on the Newton-Raphson method is briefly explained in Algorithm 1. Then, we focus on the algorithm of the implicit stress integration in the seventh line of Algorithm 1. In particular, that of the partitioned implicit stress integration (Section 3.3) is listed in Algorithm 2.
In the implementation using the partitioned coupling technique, existing programs for the stress integration of plastic constitutive models and creep constitutive models can be reused with a few modifications because these programs are used separately in the partitioned coupling technique. In addition, an integrated program of material models can be developed easily. For instance, by using macros of the C preprocessor, functions of stress integration for a single plasticity model, a single creep model, and a combined plasticity-creep model can be generated from a single program.
In the stress integration procedures, the plastic part (lines 8-12 of Algorithm 2) can be implemented separately from the creep part (lines 14-20 of Algorithm 2). If the stress is placed inside the yield surface, then only the creep part is executed. Furthermore, the plasticity model or the creep model can be changed tractably by calling another function of the program. Each function is dedicated to a specific plasticity or creep model, such as the Ohno-Wang model or the Norton-Bailey model. A function for a specific plasticity-creep model, such as the combined Ohno-Wang and Norton-Bailey model, is not needed.
Although it is also possible to calculate consistent tangent moduli in a partitioned manner, a monolithic block matrix is adopted in the present study (see Equation 85). However, the advantage of separate implementation remains, even in the monolithic tangent modulus system. The plastic submatrix, t+Δt L p , and the creep submatrix, t+Δt L cr , can be generated separately for each plasticity or creep model. If the element integration point does not yield, ignoring the plastic submatrix, t+Δt L p , reduces the monolithic 6 × 6 matrix to a 3 × 3 matrix, namely, t+Δt L cr . Furthermore, when another plasticity or creep model is to be used, the plastic submatrix, t+Δt L p , or the creep submatrix, t+Δt L cr , can be replaced by that of the model to be used.

COMBINED OHNO-WANG AND NORTON-BAILEY MODEL USING THE PARTITIONED COUPLING FRAMEWORK
For the evolution law of the ith back stress in Equation (18), the most classical choice would be the Prager model, in whicḣ where h i is a material parameter. Note that M = 1 for the classical Prager model, whereas M ≥ 2 for the generalized Prager model with the approach of multiple back stresses. This equation indicates thaṫi in Equation (18) is zero. The next choice would be the Armstrong-Frederick model. 16 In this model, a dynamic recovery term is considered in addition to the strain-hardening term (Equation 88), aṡi where h i and i are material parameters for the strain hardening and the dynamic recovery, respectively. Note that M = 1 for the classical Armstrong-Frederick model, whereas M ≥ 2 for the generalized Armstrong-Frederick model. In this equation,̇i = − i̇p . Here, we adopted the Ohno-Wang model. 3,4 This model replaceṡp in Equation (89) with a more general form that corresponds to the dynamic recovery of each i separately, that is, Here, H denotes the Heaviside step function: where ⟨ ⟩ are the Macaulay brackets, that is, ⟨x⟩ = xH (x), and r i is a material parameter given as r i = h i ∕ i . Moreover, f i is defined as In the Ohno-Wang model,̇i in Equation (18) Another expression of t+Δt i is also available in the Ohno-Wang model: The latter expression is used in computation.
For the equivalent creep strain rate,̇c r , in Equation (24), we have various choices, such as the Norton law, the Bailey law,̇c the Norton-Bailey law,̇c and other creep models that use an exponential or hyperbolic function. In these equations, A, n, and m are material parameters. Among them, we used the Norton-Bailey law. The Norton-Bailey law considers both stress-hardening and time-hardening effects, enabling us to model both the primary and secondary creep. In particular, the time-hardening effect has not been considered in the previous time-dependent constitutive models with the Ohno-Wang model. 8,9,12,13 The equivalent creep strain increment, Δ cr , in Equation (45) can be calculated by integration of Equation (98) using the backward Euler integration method:

NUMERICAL EXAMPLES
The proposed partitioned coupling approach with the Ohno-Wang model and the Norton-Bailey model was implemented in an open-source parallel finite element solver, ADVENTURE_Solid. 29,30 ADVENTURE_Solid is based on the hierarchical domain decomposition method 31 and has an advantage in solving large-scale structural problems. Due to the symmetry of the consistent tangent modulus of the combined Ohno-Wang and Norton-Bailey model, 28 the conjugate gradient (CG) method was used to solve the linear system of equations in the hierarchical domain decomposition method. The Ohno-Wang model in the developed program was verified by numerical comparisons with a commercial finite element solver, ADVENTURECluster, which has the function of the Ohno-Wang model. The Norton-Bailey model was verified using a commercial finite element solver, ABAQUS. Although ADVENTURECluster has the function of a rate-dependent nonlinear kinematic hardening model, 9 it is not the same model as the combined Ohno-Wang and Norton-Bailey model in the present program. Therefore, the combined model could not be verified by numerical comparison with other solvers. The partitioned coupling framework was verified by checking the convergence criteria of Equations (55) and (56). It was also verified that the scalar equations (53) and (54) are satisfied when the convergence is attained. Using the developed program, first, a cube problem was solved to show the basic behavior of the present constitutive model. Then, a thruster problem was analyzed to demonstrate the capability of the present constitutive model and algorithm for a realistic structure in a high-temperature environment. In addition, the convergence performances of the Newton-Raphson method, the CG method, and the partitioned stress integration were investigated.

Cube problem
A cube model consisting of one hexahedral element was analyzed in order to investigate the basic behavior of the present constitutive model, such as the differences among the Ohno-Wang model, the Norton-Bailey model, and the combined model. The dimensions and the boundary conditions of the cube are illustrated in Figure 1A. The bottom surface of the cube was totally constrained, whereas the top surface was subjected to a cyclic load, the transition of which is plotted in Figure 1B. Time increment, Δt, was set to 1 s. The adopted material parameters of the Ohno-Wang and Norton-Bailey models are listed in Table 1. To demonstrate the differences among the three models, these parameters were calibrated so Norton-Bailey parameter, m 0 that both plasticity and creep appear simultaneously. Note that, in the present numerical example, we adopted the linear isotropic hardening rule: where Y 0 and H are the initial yield stress and the isotropic hardening coefficient, respectively. The computed relationships between stress zz and strain zz of the Ohno-Wang model, the Norton-Bailey model, and the combined model are plotted in Figure 2A Figure 2B. The horizontal and vertical axes represent the time and strain zz, respectively. The Ohno-Wang and combined models exhibited much larger strain than the Norton-Bailey model. Moreover, the results for the combined model differed gradually from those of the Ohno-Wang model.

Thruster problem
We analyzed a thruster model considering heat conduction and thermal expansion. The objectives of this numerical example are to show the importance of considering both a kinematic hardening plasticity model and a creep model in comparison with a single plasticity model and a single creep model and to demonstrate the capability of the present partitioned coupling framework in simulating a realistic structure at high temperature. A finite element heat conduction analysis was performed by ADVENTURE_Thermal, 29,30 followed by a finite element structural analysis by ADVEN-TURE_Solid in a one-way coupling manner. ADVENTURE_Thermal is a heat conduction solver that is also based on the hierarchical domain decomposition method. 31 A single mesh with its hierarchical domain decomposition was used in both the heat conduction and structural analyses. The temperature at each time step computed by the heat conduction analysis was used for the input of the structural analysis via files. Each time step was divided by linear interpolation. In the structural analysis, thermal strain, th , was computed from the temperature, T, by where and T ref are the coefficient of thermal expansion and the reference temperature, respectively. Moreover, we assumed instead of Equation (25).
The dimensions of the thruster model are described in Figure 3A. Based on these dimensions, we created a three-dimensional solid mesh, the sectional view of which is depicted in Figure 3B. This mesh consisted of quadratic tetrahedral finite elements, each of which had 10 nodes. The numbers of elements and nodes were 1,608,768 and 2,445,696, respectively. For parallel computing, this mesh was decomposed hierarchically into 16 parts with 1250 subdomains for each part, as visualized in Figure 3C. In this figure, the parts are placed separately, whereas the subdomains are colored. The numbers of message passing interface (MPI) processes and OpenMP threads for each MPI process were set to 16 and 24, respectively. Eight computing nodes of the supercomputer Fugaku were used.
The boundary conditions of the heat conduction analysis are illustrated in Figure 3D. They consisted of a heating process of 50 s and a cooling process of 550 s. In the heating process, a fixed temperature was prescribed on the inner surface of the combustion chamber. The fixed temperature on this surface was 800 • C, whereas the initial temperature at other locations was 20 • C. This fixed-temperature boundary condition was a simplified model to simulate the contribution of combustion. The temperature in the vicinity of the inner surface was assumed to increase in a very short time. After that, the fixed temperature remained 800 • C during the heating process (50 s). At the end of the heating process, this boundary became adiabatic. In the heating process, temperature is expected to rise, leading to high thermal stress and large plastic and creep deformation. In the cooling process, zero or very small plastic deformation and small creep deformation are expected. The time increment was set to 5 s in both the heating and cooling processes. Hence, the number of time steps was 120. In addition, the heat-transfer boundary condition was applied to the outer surface of the thruster throughout the analysis. The heat transfer coefficient and the ambient temperature were assumed to be 1.00 × 10 −2 W/(mm 2 K) and 20 • C, respectively. Furthermore, the left-hand side of the combustion chamber was constrained in the structural analysis.
As the material, Ti-6Al-4V was assumed. The adopted material parameters are listed in Table 2. For the heat conduction analysis, a single set of material parameters was used. In contrast, multiple sets of material parameters were tested for the structural analysis to demonstrate possible behaviors. Moreover, we tested the Ohno-Wang model, the Norton-Bailey model, and the combined model to clarify the differences among the three models. at the end of the cooling process (600 s) set of material parameters was used. These parameters are from Kourousis et al. 33 Then, for the Norton-Bailey model, we adopted three sets of material parameters. Hence, the combined model involved three sets of material parameters. The three sets are referred to as materials #1, #2, and #3. The parameters for material #1 were calibrated to approximately reproduce the history of the equivalent creep strain of Ti-6Al-4V 32 under a constant equivalent stress of 593.5 MPa for 0 h to 200 h. Although the difference between the curves obtained by material #1 and the experimental result was small for 593.5 MPa, the difference became large when the magnitude of the equivalent stress was changed. This is a limitation of the Norton-Bailey model, and the combined model is not fully validated for Ti-6Al-4V. A creep model that is more suitable for Ti-6Al-4V can be used for Equation (24) and is easily implemented using the proposed framework. Materials #2 and #3 are modified versions of material #1. Figure 4A shows histories of the equivalent creep strain for materials #1, #2, and #3, and those obtained by the empirical equation, 32 which was fitted to the experimental results. Note that the histories are for 0 s to 600 s. The lines for materials #1, #2, and #3 follow the integral of Equation (98). The differences among materials #1, #2, and #3 are the effects of stress hardening. If the equivalent stress is low, the equivalent creep strain rate of material #3 is the smallest among them. If the equivalent stress is high, then it is the largest. Similarly, Figure 4B shows histories of the stress under constant strain. In this example, material #3 exhibits the largest amount of stress relaxation. Figure 5 visualizes the distributions of the computed temperature at the end of the heating process (50 s) and the cooling process (600 s). Moreover, the transition of the temperature is plotted in Figure 6. The horizontal and vertical axes represent the time and the temperature, respectively. The solid and dashed lines denote the temperature at points A and B, respectively, in Figure 3D. The temperature of the combustion chamber part was very high in the heating process and decreased in the cooling process. The results of the heat conduction analysis were used in the structural analysis. To set the initial temperature to 20 • C, one time step was added to the beginning of the 120 time steps. Hence, the number of time steps of the structural analysis was 121. The time increment of the first step was assumed to be 1 × 10 −3 s. Figures 7-9 show the distributions of the equivalent stress of the Ohno-Wang model, the Norton-Bailey model, and the combined model, respectively, at the final time step. In all of the models, the equivalent stress at the fixed end, which is plotted in Figure 10, was high. The Ohno-Wang model, the Norton-Bailey model with material #1, and the combined model with material #1 exhibited similar distributions of the equivalent stress, although the Norton-Bailey model showed higher equivalent stress at the fixed end. However, in the Norton-Bailey model and the combined model, the equivalent stress of material #2 was lower than that of material #1, due to stress relaxation. Moreover, that of material #3 was much lower. Figures 11 and 12 depict the distributions of the equivalent plastic strain of the Ohno-Wang model and the combined model, respectively. The equivalent plastic strain at the fixed end is plotted in Figure 13. Similar to the equivalent stress, the distributions of the equivalent plastic strain of the Ohno-Wang model and the combined model with material #1 were similar. However, those with materials #2 and #3 were smaller. Figures 14 and 15 visualize the distributions of the equivalent creep strain of the Norton-Bailey model and the combined model, respectively. The equivalent creep strain at the fixed end is plotted in Figure 16. In both the Norton-Bailey model and the combined model, material #1 exhibited very small equivalent creep strain. Those of materials #2 and #3 were larger. In summary, material #1 exhibited a plasticity-dominant behavior, whereas material #3 exhibited a creep-dominant behavior. Material #2 showed both plasticity and creep behaviors, in which the equivalent creep strain was slightly The results for material #2 were investigated in detail. The transitions of the equivalent stress, the equivalent plastic strain, and the equivalent creep strain at point A in Figure 3D are plotted in Figure 17A

ACKNOWLEDGMENTS
The authors would like to thank Dr.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.