Global analysis and optimal harvesting for a hybrid stochastic phytoplankton-zooplankton-fish model with distributed delays

1.   Introduction

Bainitic press hardening is a key production technology for advanced high strength steels, specifically in the automotive industry ^[1]. Structural parts like A- and B-pillars, roof rails or bumpers are produced based on this method. Currently, it is of specific interest to the industry to take this important technology to the next level in terms of process optimisation. This technique, using a combination of interrupted cooling and quasi-isothermal holding in the bainitic parameter regime, offers highly desirable strength-failure-strain combinations ^[2]. The possibility to control the process tool temperature offers a vast amount of possible advances in process technology, and work pieces with locally varying thermoelastic properties for a wide range of applications are one of the central goals in this context. A challenge of central importance to the bainitic press hardening process is the understanding of the bainitic transformation.

The mechanism of the bainite transformation is the subject of ongoing debate for decades, see e.g., ^[3,4,5,6]. The main controversy focuses on the following two arguments: (ⅰ) the nature of the bainitic ferrite/retained austenite interface and its capacity to support the displacive mechanism of bainitic ferrite growth and (ⅱ) the role of carbon transport during the bainite transformation, i.e., whether or not bainitic ferrite is initially supersaturated with carbon and then carbon partitioning takes place from the bainitic ferrite. Within this publication, we will focus on the transport-governed regime, i.e., the transformation is limited by carbon redistribution. The main branches of the scientific discussion about the details of the transport-governed regime can be summarised as follows.

Considerations about the carbon partitioning published in ^[7,8] suggest that the ferritic component of bainitic microstructures can form with maximal carbon supersaturation, while the escape of carbon from ferrite to austenite after transformation is the essential aspect of the carbon transport. Alternatively, the transport of carbon during diffusional transformation of austenite to ferrite under metastable equilibrium conditions at the interface is also considered (e.g., ^[5]). An intermediate scenario, where ferrite grows with a partial supersaturation of carbon, with the remaining carbon partitioning into austenite or forming carbides, is described e.g., in ^[9,10,11]. All of these hypothesises effectively consider the competition between carbon escape and interface migration at high supersaturations, thus the transformation might exhibit a higher velocity than is expected by equilibrium partitioning of carbon. For our investigations, we focus on a regime where the diffusional transport of carbon from the supersaturated austenitic phase at local equilibrium is the rate-limiting mechanism.

The following parts of this publication are organised in four sections; a basic theoretical treatment of solid-solid transitions in the displacive and the diffuse-displacive case is given in Section 2. We also introduce the basics of the boundary integral method, which is used for the investigations on the terminal subunit growth regime.

In Section 3.1, we discuss the influence of coherency stresses and capillarity stresses on the growth of banitic subunits of small radii in the initial growth regime of the diffusional-displacive transition.

In Section 3.2, we present calculations for the specific scenario of bainitic subunit growth frustration, which results form the decreasing temperature during the continuous cooling of the bainitic transformation.

In Section 4, we discuss the context of other related growth regimes also relevant to elastically influenced solid-solid transformations.

2.   Materials and method

The work in this article is based on theoretical work and aims at elementary aspects of growth behaviour in bainitic subunits. As representative geometry, we present Figure 1.

The bainitic transformation is an example for diffusional or mixed diplacive-diffusional solid-solid transitions strongly influenced by elastic effects. An overview of elastical-diffusional transformations can be found e.g., in ^[12,13,14], in our studies reported in ^{[15,16,17,18,19,20,21,22]}, we focused on transformations exhibiting coherent interfaces. Though the coupling of displacive and diffusional transitions has been investigated recently also emphasising engineering aspects, see ^[23], our focus in the course of this publication is on the possible dominance of elastic effects in the initial and terminal regime of subunit growth. When we refer to the terminal growth regime of a bainitic subunit, we consider the freezing-out carbon diffusion upon the continuous decrease of temperature during the progressing bainitic transformation.

Within this manuscript, we focus on two distinct scenarios. First, the initial growth regime of the subunit is considered. The corresponding model is a sharp interface model, which aims at the prediction of the geometric properties of the growing pattern. Second, we investigate the terminal regime of subunit growth, where we apply the boundary integral method. As the main part of the work is associated to the terminal growth regime, we present here the theoretical basics of the boundary integral method in the context of the physical model we assume for the terminal growth regime. As the boundary integral method is especially useful in the study of steady-state moving fronts, we focus on these steady-state scenarios, here, as systematically integrated in the non-isothermal context of the bainitic transformation. The growth predictions obtained from such steady-state scenarios can be interpreted as valid during the persistence time scale which is basically limited by the variation of the governing parameters of the process.

To introduce the boundary integral method, we focus on the diffusional aspect of the transformation. We represent the diffusion equation for a field $u$ in a frame of reference which is co-moving along the y-direction symbolically via a linear differential operator $\mathcal{L}$:

and use the Green's function $\mathcal{G}$ of this operator ^[15] to express the symbolic solution for the diffusion field via Lagrange's identity:

where $\mathcal{L}^\dagger$ is the adjoint operator of $\mathcal{L}$ and $\bar{\mathcal{J}}(u, \mathcal{G})$ is a function of the diffusion field and the Green's function with derivations of order not higher than one in the case of the diffusion equation. Apparently, when we integrate the Lagrange identity and have no volume inhomogeneities, the integral contributions reduces to boundary terms. In our case, the inhomogeneity is the mass conservation of carbon at the moving interface, as we discuss in the subsequent description of the physical model. Finally, the whole description reduces to a fully interfacial problem as the obtained symbolic solution expresses the value of the diffusion field on each point at the interface when the local equilibrium condition is described, see also the following description of the physical model.

For the terminal growth regime, the basis of our considerations is a coupled model of diffusional carbon transport and static elastic misfit stresses within the two adjoining phases, i.e., the parent austenite and the growing ferrite phase, together with the physically necessary boundary conditions at the progressing phase boundary. The diffusional transport is governed by a diffusion equation in the bulk, complemented by the mass conservation condition and the Gibbs-Thompson condition of local phase equilibrium. In total, we obtain the following three equations for the dimensionless carbon concentration field $u = (c-c_0)/\Delta c$, with $c$ being the non-equilibrium carbon concentration field, $c_0$ is a reference concentration and $\Delta c = c_{\gamma}^{eq}-c_{\alpha}^{eq}$ denoting the difference in the equilibrium carbon concentrations of the $\gamma$-austenite and the $\alpha$-ferrite phase at the respectively imposed constant temperature

By $D$ we denote the carbon diffusion constant, $T_{eq}$ is the transition temperature, $L$ is the effective latent heat ^[24,25]. The capillary length $d$ is defined here as $d = \gamma T_{eq}/ L m \Delta c$, where $\gamma$ is the surface energy, $\kappa$ is the interface curvature and the dimensionless undercooling is given by $\Delta = (c_{eq}-c_0)/(c_{\gamma}-c_{\alpha})$. The slope of the coexistence line is denoted by $m = \| dT/dc \|$. In accordance with our previous studies on displacive transformations, we choose the class of invariant plain strain transitions here, as it allows for extended coherent interfaces. When we indicate the mother phase via $\alpha$ and the onsetting phase via $\beta$, the elastic contribution to the local equilibrium equation is defined as

Here the indices ${ik}$ refer to the components of the strain and eigenstrain tensors, with ${ss}$, ${s\tau}$ and $\tau\tau$ denoting the two tangential components in natural coordinates. This is obtained by introduction of the appropriate free energy and assuming phase and mechanical equilibrium at the coherent interface, see ^[15] for details. We require that the eigenstrain tensor $\epsilon_{ik}^{(0)}$ allows to satisfy coherent interfaces, the eigenstrain stress tensor $\sigma_{ik}^{(0)}$ is the stress tensor corresponding to Hooke's law for the eigenstrains $\epsilon_{ik}^{(0)}$. Therefore, the first contribution on the r.h.s of Eq 4 involves the eigenstresses and the strains in the austenitic phase. Following this brief introduction of our theoretical basis, we continue with the discussion of the initial growth regime.

The physical picture emerging from these equations is clear. The propagation of the interface releases carbon proportional to the normal component of the growth velocity $\upsilon_n$, while the local transition point itself depends on the interfacial stresses arising from the coherent lattice transformation and the curvature. The transport of the carbon from the interface limits the propagation velocity. We note that this model represents the coupling of the elastostatic and the diffusional problem via the local equilibrium condition for the carbon concentration at the transformation interface alone.

3.   Results and discussion

3.1.   Pattern selection in the initial growth regime subject to elastic softening

The transformation of austenite to bainite involves displacive deformations of the matrix, diffusion controlled growth of bainitic sheaves and the precipitation of carbides inside the ferrite (lower bainite) or inside the austenite (upper bainite). A successful prediction of the transformation kinetics requires consideration of the complex interaction of these processes: elastic/plastic deformation, carbon partitioning and carbide precipitation. The local driving force for the transformation depends sensitively on the carbon concentration, the local misfit-stresses and the temperature. During the transformation, the carbon in the supersaturated bainitic-ferrite can either be released in the surrounding austenite phase or be captured by the embedded carbides, which are competitive for the solute atoms.

When we focus on the very initial regime of the solid solid transformation, involving a nucleus of nanometer size, the size dependence of the elastic parameters in this stage can lead to substantial elastic softening of the growing phase. Furthermore, on this scale, the bulk crystal symmetry structure can be suppressed due to the dominance of interfacial effects on the overall energy of the nucleus. Therefore, we approximate the growing inclusion by concerning only the hydrostatic, isotropic elastic response to stresses from the surrounding bulk phase. The assumption of isotropic elasticity exhibits the lowest expected error on the real elastic behaviour of the nucleus prior any information on the early $nm$ growth regime. Thus, we consider this initial stage in analogy to the argumentation given in ^[12], which provides an approximative and analytically accessible description of the growth problem.

We assume that the new phase, which contains N particles and occupies a volume W, is homogeneous with chemical potential $\mu$ and a pressure P (we assume that the external pressure is zero). Because of the conservation of mass, we have the following expression for the volume of the melt:

where $v_{\alpha}$ and $v_{\gamma}$ are the atomic volumes of the ferritic and austenitic phases, and the last term describes the change in the volume of the nucleus upon deformation ($u_n$ is the normal component of the displacement vector at the interface). The subunit is assumed to be a very oblate lentil of radius R and height h $\ll$ R. In a first approximation, we can ignore the height of the lentil in solving the elastic problem. In this case the pressure P, exerted on the nucleus by the retained austenite, is given on a circular plane cut within the radius R. This problem is equivalent to a crack problem and its solution provides the normal component of the displacement vector at the interface, see ^[26]. In this scenario, the equilibrium Gibbs-Thomson condition at the interface reads as ^[27,28]

Here $f^0_{\gamma}$ is the free energy density of the undeformed (initial) crystal and we have introduced cylindrical coordinates, $\bar{\sigma}_{ik}$ is a normalised stress tensor, $\nu$ and $E$ are the Poisson and Young coefficients, $\gamma$ is the interfacial energy and $\kappa$ is the curvature of the interface (considered positive for a convex solid). To focus on the role of the elastic effects of the transformation, Eq 6 approximates the iron chemical potential by reducing contributions depending on the carbon content to an effective lattice expansion which is taken into account via the atomic volumes $v_\alpha$, $v_\beta$.

Following the argumentation in ^[29] and ^[12], the radius $R_c$ of the critical nucleus is related to the pressure $P$ by the Griffith formula, leading to a relation $\Delta(R) = \Delta (1-\sqrt{1-R_c/R}) = u_{Int}$ for the interfacial equilibrium condition. This relation together with mass conservation provides the characteristic length scale $\rho$ and the height $h$ of the inclusion

Apparently, these expressions remain valid only when the elastic softening is pronounced and the assumption that neglecting the nonydrostatic aspect of the elastic behaviour in the inclusions is acceptable. However, within the approximative validity of these assumptions, we can consider the onset to be elastically dominated. The crucial point is the absence of any thermal quantities in $\rho$. This indicates that the pattern is dominated by elastic effects in this early growth regime. After the shape undergoes a Mullins-Sekerka instability, which takes place as derived in ^[30], it grows via the coupled diffusional-displacive mechanism, qualitatively comparable to the dendritic growth regime as we described it in ^[19]. In the next section, we introduce the regime of frustrated growth of bainitic subunits, which is represented by a limiting case of the elasto-chemical diffusion-limited growth, i.e., diffusion-limited growth mediated by an elastic selection.

3.2.   Elastic pattern selection in frustrated bainitic subunit growth: zero velocity limit

Within this section, we focus on the elasto-chemical diffusion-limited growth of a subunit in the terminal regime of carbon transport frustration, i.e., where an effective reduction of the thermodynamic driving force due to the size of the grown subunit can be assumed. The frustration of the carbon transport within the retained austenite is related to the continuous decrease of the overall temperature during the progression of the bainitic transformation. Upon a linear temperature decrease the carbon diffusion coefficient undergoes an exponential decay, as it obeys an Arrhenius-type temperature dependence. The combination of a reduced driving force and an exponentially decreasing diffusion coefficient suggest to study the zero velocity or zero diffusion limit.

To take full advantage of our assumptions, the boundary integral method, as introduced in the methods section, see Eq 3 and in ^[19], is most suitable. The resulting formulation of the problem combines the governing equations of the elastostatic and transport problem into a closed representation for the unknown boundary, $y(x)$, and the tip radius of curvature $R$:

To simplify the representation, we introduce a splitting of the elastic contribution to the local equilibrium into a constant, hysteresis contribution $\Delta_{el}$ and the remaining, spatially varying contribution, for details we refer to ^[19]. Important here is to note that the Peclet number of the problem depends on the net driving force which takes the elastic hysteresis into account, i.e., $\tilde{\Delta} = \Delta - \Delta_{el} = \sqrt{p\pi} exp(p) erfc(p)$. This includes that a strong elastic hysteresis leads to a substantially reduced driving force for the transformation, and we use approximations for the dependence of the effective driving force on the Peclet number and the Green's function of the problem as described in ^[19] to obtain the following representation:

Here we define $\Phi = 2(1-\nu^2)/(E(\epsilon^2_{yy}+\epsilon_{zz}^2+2\nu\epsilon_{yy}\epsilon_{zz}))$, function $B(\vec{x})$ includes the local contributions to the elastic free energy difference. The stability parameter of classical dendritic growth is $\sigma = d_0/pR$. We note that the velocity and the radius of the asymptotically matching Ivantsov parabola are defined as

As mentioned, the specific scenario we apply the boundary integral formulation for is a scenario of very slow diffusion, so that the interface propagates slowly. Then, the diffusion is sufficiently low to have the elastic effects dominate the transition. In the case of very small $\Delta_{el}/p$, one shall exploit the possibilities of analytical methods, e.g., apply the complex matching method which yielded the solution of the classical dendritic problem, but for the limit of p/$\Delta_{el} \ll$ 1, which represents the scenario of transport and frustrated transformation in the subunits, analytical approaches in the aforementioned spirit are not at hand, suggesting specific numerical consideration here.

Therefore, we approximate this frustrated transformation regime via the limit of zero velocity and diffusion. Formally, we can reach the regime of a possibly negligible release of latent heat at the interface when $\upsilon \to 0$. The Peclet number is given for a fixed $\tilde{\Delta} \ll 1$ as

which shows that for a given radius of the subunit, i.e., if $R(\tilde{\Delta}) \sim 1$, the limiting case of a decreasing diffusion coefficient $D \to 0$ is equivalent to $\upsilon \to 0$. As the exponential dependence of the carbon diffusion coefficient on the temperature provides a pragmatically reasonable approximation to the strong decay of carbon diffusion, we will consider if a scenario with $R(\tilde{\Delta}) \sim 1$ can be predicted.

We introduce a representation of the problem which suits this approach, namely

and the corresponding solvability parameter reads as

We note that we have introduced a solvability parameter which describes the growth regime without involving any diffusional transport parameters. This reflects the possible elastic selection on the formal level. We implemented a symmetric diffusion model instead of the one-sided model which is typically used for highly asymmetric carbon transport. Therefore, we note that the dependence of the stability parameter on the control parameter in the symmetric and one-sided model is qualitatively identical. The symmetric model implementation is numerically easier to handle in the regime of $p/\Delta_{el} \to 0$. The resulting calculations for pure shear can cover the regime up to $p/\Delta_{el} = 0.05$, but even at this already small value, an extrapolation to $p/\Delta_{el} \to 0$ remains inconclusive. As the plot in Figure 2 shows, we can at least assume that the hidden behaviour of the selection does not obey a linear law. Exemplarily, we discuss the result of a potentially present finite $\sigma_{el}$ and the result for the case $\sigma_{el} \sim (p/\Delta_{el})^{1/2}$.

A finite value of $\sigma_{el}$ at zero $p/\Delta_{el}$ would lead to an elastic selection as

where we note that the net driving force relative to the elastic hysteresis in our model is assumed to be reduced by the frustration of the transformation. A possible square-root dependence of the solvability parameter on $p/\Delta_{el}$ would lead to scaling laws for the transformation velocity and the radius as

Despite of the inconclusive numerical results for the zero velocity limit concerning a finite $\sigma_{el}$, the results show that an elastically dominated terminal growth regime of the subunits is possible.

4.   Conclusion

We presented our recent work on the possibilities of elastically dominated regimes in diffusion-limited solid state transformations. Referencing the bainitic transformation, specifically the scenario of sheaf growth via subunit growth into the retained austenite, the very initial and the terminal growth regime are considered. For the early growth regime, we find that under strong assumptions for the elastic softening and effective isotropy of the onsetting phase, the problem can be considered analogous to a crack formation problem. This is in line with previously published investigations on melt inclusion growth published in ^[29]. When the Mullins-Sekerka instability sets the transition from this early, elastically governed growth regime to the coupled diffusive-displacive growth regime, described in ^{[15,16,18,19]}, the transport of carbon in the retained austenite is assumed to be the limiting factor of subunit growth. Taking into account the decrease of the carbon diffusion coefficient with decreasing temperature suggests to consider the effect of a reduced diffusion coefficient. In combination with elastically frustrated growth, i.e., reduction of the effective driving force, the reduced diffusion coefficient suggests to study a zero velocity limit. The presented numerical results show the possibility of an elastically dominated terminal growth regime of subunits.

Acknowledgments

N. T. would like to acknowledge financial support from Innovation Foundation For Postgraduate and Fundamental Research Funds of Central South University (Grant No. 1053320171074), Changsha, China. K. W. would like to acknowledge financial support from Chinese Scholarship Council (Grant No. 201706370221), China. X. Y., M. F. and C. H. acknowledge the support of the Deutsche Forschungsgemeinschaft under Project No. SPP 1713.

Conflict of interest

All authors declare no conflict of interest in this paper.