Implementation of a reaction-diffusion process in the Abaqus Finite Element software

,


Introduction
Simulating the effect of impurities on the integrity of structures leads to account for several interactions between, e.g., the mechanical fields, the impurities transport and trapping, the thermal fields, etc.The simulation of all these phenomena simultaneously is a complex task, especially when strong couplings are involved or investigated: in the hydrogen embrittlement of metals [1], or in the hydrolysis of polymers [2,3], for instance, mobile species are adsorbed, transported through the material, and trapped on specific sites whose density is time and space dependent [4,5] (e.g,, through the development of plasticity for hydrogen in metals [6]).Furthermore, mechanical fields can be affected by these species because of induced deformations or through modifications of mechanical properties [7,8].Numerous studies account for such interactions in finite element codes, in various application fields (metal/hydrogen, water/polymer, metal/lithium ions, see [9][10][11][12][13][14][15][16][17] among others), but very few developments include several phenomena in the computations [18,19], especially in the commercial finite element codes, due to their inherent limitations in terms of available degrees of freedoms at each node.Such an inclusion may, however, be of importance, e.g., to model the behavior of structures in the presence of both impurities and evolving thermal boundary conditions [20].The aim of this work is thus to introduce some developments performed in Abaqus to solve coupled mechanicalmultidiffusion finite element problems.This paper is limited to a reaction-diffusion process between two species, which is solved by using a coupled mechanical-diffusion scheme ('coupled tempdisplacement' in Abaqus) that allows further developments to account for the mechanical fields as well.First, the multidiffusion implementation strategy is presented, and then an application to the Gray-Scott reaction-diffusion model is presented to illustrate the new capabilities [21,22].

Introduction of a multidiffusion process in Abaqus
In order to solve a complex problem with mechanics and multidiffusive fields in a finite element (FE) software, it is mandatory (i) to have a finite element formulation that includes as many degrees of freedom (DOFs) per node as the number of unknown fields, and (ii) to introduce the correct weak formulation for all of these DOFs for solving the problem.Introducing extra DOFs is complex; one may exploit the unused mechanical DOFs (rotations, numbered from 3 to 6, or the third displacement component in 2D problems), or add extra features to the elements (see [23,24] for phase field implementation in Abaqus) through a user element (UEL) routine [25].One approach of particular interest has been proposed by Chester [26] to solve coupled thermo-chemo-mechanical problems in polymers (this work has been applied in [27] for a simple adsorption process).In this work, a UEL has been developed that activated an extra DOF (not numbered between 1 and 6, for displacements and rotations, nor NT11), in addition to the introduction of a relevant weak formulation as specified in [25].Such DOFs are included by default in the Abaqus element library for 'coupled tempdisplacement' procedures, but they are hidden and cannot be accessed through the CAE interface or input files1.These DOFs, numbered from 12 to 30, correspond to NT (for 'Nodal Temperature') variables.Once activated by the UEL routine, their boundary conditions can be imposed in the input file and their values (NT, HFL, etc.) can be required in the output database file.
It is worth noting that all the studies mentioned above, where an UEL was used to redefine the problem, have also superimposed additional layers of elements taken from the Abaqus library in order to visualize the results.As demonstrated in [29], it is possible to go further and extend the Abaqus finite element formulation by superimposing a user element to an Abaqus element: the terms that are not included by default in the formulation are introduced through the UEL routine and the Abaqus material library.The approach chosen in the present study combines the advantages of keeping the features of the Abaqus libraries (materials, elements, etc.) and of adding extra terms and DOFs in the finite element formulation by using a superimposed UEL.Thus, the implementation work is optimized because the mechanical behavior needs not being redefined.Even if a multidiffusion process only is considered here, the ultimate goal of a fully coupled mechanical-multidiffusion problem has been kept in mind during the developments.

Implementation process
Our strategy is presented in Figure 1: several element layers sharing the same nodes are defined, and a 'coupled temp-displacement' procedure is used.In this example, the three UEL layers have the same numbers of DOFs and, assuming that DOFs 11, 12, and 13 represent diffusion DOFs between which a reaction may occur, all user elements layers share the same UEL routine with different parameters.Each layer, in this example, has a specific role: (i) Layer 1: the Abaqus element (with mechanical DOFs 1 to 3, and 11) involves the mechanical behavior, one diffusion phenomenon (related to DOF 11), and its effects on the mechanical behavior.The problem is strongly coupled (i.e., the diffusion and the mechanical problems are solved simultaneously), but no effect of mechanics on diffusion is possible here (except with developments beyond the scope of this work).(ii) Layer 2: this UEL layer activates DOF 12 through its related weak formulation (here, diffusion, but it could be any other physical or chemical process), and the coupling between mechanics and DOF 11 (for no complete strong thermo-mechanical coupling is included by default in Abaqus).It is worth noting that other approaches can be considered in the superimposition process (for instance, a single UEL can be used to activate DOFs 12 and 13, and to introduce all the ingredients needed in Abaqus).Each element layer leads to the computation of a specific stiffness matrix, performed either by Abaqus or by the UEL, as shown below: (1) In the case presented in Figure 1, the stiffness matrices are 16×16 and 20×20 for the Abaqus element and for the UEL, respectively.At the end of the superimposition process, the stiffness matrix of the global problem is 24×24: due to the activation of the extra nodes, the initial Abaqus element stiffness matrix has increased significantly, without any other user manipulation than the activation of hidden DOFs.This strategy is applied below, where only DOFs 12 and 13 are considered, for illustration.The transient 'coupled temp-displacement' procedure is used, even if there is no coupling between DOFs (12,13) and (1,2,3,11) in the present work.

Application
The Gray-Scott model is considered here as a test reaction-diffusion process to be implemented.

The Gray-Scott model
The Gray-Scott (GS) reaction-diffusion model represents a particular case of Turing systems [30], where the reactions of three chemical species are focused on.These species, , , and , define an autocatalytic system so that [21,22] The space-time evolution of species  and  can be obtained by solving the following system of differential equations: where  and  denote the concentrations of species  and , respectively,  ! and  !represent their diffusion coefficients,  is the feed rate for  and  the kill rate for .This reaction has been widely studied as a simple model to reproduce the patterns observed in several chemical reactions (not to speak of natural patterns [31]), as illustrated in Figure 2.
ve Modeling of Turing Patterns J. Phys.Chem.A, Vol.103, No. 12, 1999 1797 observed experimentally, but at such high values of [ClO 2 ] 0 , the pattern amplitude becomes so small that any pattern would be difficult to detect experimentally.This problem is actually encountered when [ClO 2 ] 0 > 1.4 × 10 -3 M; for these values, no patterned states were detected, even for large PVA concentrations.
There seems to be some discrepancy between the model  patterns which we strictly associate with non oscillatory states of the CSTR, and the nonstationary states including traveling waves and various oscillatory states which may or may not be associated with oscillations in the CSTR.
Considering the extreme simplicity of the LE model compared to the actual chemical kinetics, there is a striking agreement between computational and experimental results, at least when observed experimentally, but at such high values of [ClO 2 ] 0 , the pattern amplitude becomes so small that any pattern would be difficult to detect experimentally.This problem is actually encountered when [ClO 2 ] 0 > 1.4 × 10 -3 M; for these values, no patterned states were detected, even for large PVA concentrations.
There seems to be some discrepancy between the model calculations and the experiments at high [PVA] 0 and low [ClO 2 ] 0 below the oscillatory domain.In the experiments, when [MA] 0 increases, the Turing state is followed by a new uniform steady state, characterized by a lower value of complex [SI 3 -] (clear color).No similar computed bound appears in Figure 9.As a matter of fact, when [MA] 0 increases, the computed stationary pattern becomes unstable and is replaced by an oscillatory state.After this ͓H 2 SO 4 ͔ 0 jump, many F-state bubbles are trapped and the dynamic of the pattern is then ruled by that of the F-state bubbles.These may divide, pulsate, and eventually die by overgrowth with the spontaneous development of the M-state in their middle or a local asymmetric commutation of a ͑−͒front to a ͑+͒front, which deforms the circular symmetry of the bubble.The system no longer settles into a stationary pattern mode.Note that the spot pattern is inserted in a large discoid M-state domain disconnected from the mask ͑except at one point, at 12 o'clock in the snapshots-Fig.10, probably due to some local "boundary imperfection"͒.The diameter of this discoid M-state domain undergoes small-amplitude oscillations with a period of about 24 min, similar to that of the oscillatory F-state spots.The sequence of snapshots in Fig. 11 details the breathing dynamics of a small F-state spot of Fig. 10.Note the color changes inside the spot when it grows or shrinks, indicating that the pH decreases when the spot shrinks and vice versa.
We have also explored the effect of ferrocyanide feed concentration on pattern formation in the disk OSFR.For simplicity reasons, let us only consider the case of stationary patterns at ͓PA͔ 0 = 4 mM and ͓H 2 SO 4 ͔ 0 = 2.86 mM.In these conditions, the gel is in the stationary F-state when ͓K 4 Fe͑CN͒ 6 ͔ 0 Ͻ 5 mM.At ͓K 4 Fe͑CN͒ 6 ͔ 0 = 5 mM, the gel switches to the stationary uniform M-state; no nontrivial patterns are observed.Above a critical ͓K 4 Fe͑CN͒ 6 ͔ 0 around 10 mM, the M-state become unstable and a stationary M-state labyrinthine pattern forms, as illustrated in Fig. 12͑a͒ for ͓K 4 Fe͑CN͒ 6 ͔ 0 = 15 mM.The area covered by the M-state and the interconnections in the network increase with increasing ͓K 4 Fe͑CN͒ 6 ͔ 0 ͓Fig.12͑b͔͒.At ͓K 4 Fe͑CN͒ 6 ͔ 0 = 30 mM, similarly to high ͓H 2 SO 4 ͔ 0 ͑Fig.9͒, the stationary stability and the system is dominated dynamics of F-state spots ͓Fig.12͑c͔͒.Test experiments were also made perature.At ͓PA͔ 0 = 2 mM and ͓K 4 Fe tionary lamellae patterns could be obs by adjusting the value of ͓H 2 SO 4 ͔ 0 .T domain shifts to higher acid feed wit ture.Noteworthy, at 22 °C, no CSTR but only bistability. 33,34

DISCUSSION
In this paper, we revisit the condi tions of reaction-diffusion patterns ferrocyanide-iodate-sulfite reaction in reactor.Despite the fantastic patternin tem ͑due to morphological front insta Bloch front bifurcation͒, it was left years, contrary to the case of the CI some unnoticed and thus uncontrolled gels, the original experimental observa reproduce.Lets us point out the simil between the original and our present Both sets of experiments were made chemical feed parameters and tempera dence time in the original experimen 240 s, which is shorter by a factor 6 Conversely, the volume of our stirred cantly greater to minimize feedback eff the gel on the actual composition of the mize global coupling effects.In this t coupling cannot be totally avoided but ments it does not drive the overall pa main differences come from the polym sets of experiences and the way they a the CSTR contents.The original ex weight polyacrylamide gels films ͑0.2tact with the CSTR contents through a a nitrocellulose membrane to provide, rigidity and a white backing for visu such composite of porous membrane w adsorption properties.In our experim material is used-agarose ͑4% weigh ized gel that should make experim easier.Our gels are thicker ͑w = 0.75 m etry OSFR and 1.0 mm in the annular ment can be very sensitive to the thickn original set of publications 38 and in rec bistable systems. 47-49Li et al. 38  CDIMA (chlorine dioxide, iodine, malonic acid) reaction after 12, 20, 35 and 52 minutes [34].

Numerical implementation
The patterns induced by the GS model have been the subject of numerous studies from the seminal work by Pearson [36] (see, e.g., [37][38][39][40][41][42]), including many for entertainment purposes 2 , and a classification of the GS patterns has been proposed (see Figure 3), depending on the (, ) values.Consequently, many implementations of the GS reaction can be found, based on finite differences and forward Euler integration scheme for efficiency reasons ( [43][44][45], among others, and the very complete webpage of R. Munafo [46]), mainly in 2D.Very few [47][48][49] apply the finite element method, especially Abaqus.One study [50] includes mechanical coupling, but no indication on the implementation process is given, nor if extra DOFs have been introduced, unfortunately.
We have implemented the GS reaction in Abaqus by introducing DOFs 12 and 13; the details of the RHS vector and of the AMATRX matrix have been adapted from [48] by considering constant diffusion coefficients, in particular.Computations have been performed with the 'coupled tempdisplacement' procedure, even if no mechanical nor temperature field is computed.In order to evaluate the ability of our implementation to simulate a GS process accurately, all the results are compared with those given by the Python script written by D. Bennewies [44].
) as defined in [36].For (, ) points where no pattern is specified, a constant homogeneous field for  as well as for  is expected.
2 For instance, 'Gray-Scott reaction diffusion' keywords in YouTube gives 779 results.the CSTR corresponding to the vicinity of this critical point. (36)Furthermore, patterns develop behind a front that propagates into the previously existing uniform state that has now become unstable.What is remarkable is that this front exhibits morphological instabilities giving rise to growth modes involving spot division or finger splitting. (36)Figure 6 illustrates the finger tip splitting growth mechanism for Turing pattern which ultimately leads to stripes.

Bistable CSTR
When the CSTR evolves in its bistable region, as it can be the case for the CDI (37) (or also FIS (38) ) reaction, a first important aspect is the determination of the possible corresponding states in the gel.Let one consider the situation along the direction orthogonal to the CSTR-gel reactor boundary, i.e., along the depth of the gel.If the CSTR is in the F branch, at each point along the gel, fresh reactants are brought by diffusion from the feeding edge, where the concentration is kept fixed.Close to this edge, the extent of the reaction is small and the chemical composition remains close to that of the flow branch.As we move away from this edge, the extent of the reaction becomes larger because the amount of fresh reactants that reaches the corresponding space point is limited by its transport through molecular diffusion.So, if the gel film is thick enough, the regions of the gel far from the feeding edge may eventually belong to a state laying on the T branch.In such a case, the composition changes from branch F to branch T somewhere inside the gel.Thus, for the same F state in the CSTR, one may observe two quite different composition profiles as a function of w, the thickness of the gel.If w is very small, the chemical composition in the gel Fig. 7.4 Spiral waves in the BZ reaction.The pattern develops after a break of a circular wave.The same intervals between snapshots (a-h) and the same initial concentrations of reactants as in Fig. 6.6.Reproduced with permission from [17] ter of the grid was then U = ln,V = 114).These e then perturbed with f 1% n order to break the square system was then integrated e steps and an image was ases, the initial disturbance tward from the central patterns in its wake, until was affected by the initial tion.The propagation was h the leading edge of the oving with an approximately ty.Depending on the paramok on the order of 10,000 to ps for the initial perturbation he entire grid.The propagaf the initial perturbation is r of 1 x space units per er the initial period during urbation spread, the system ymptotic state that was either ent or time-dependent, deparameter values.(8), which occurs in the vicinity of a Hopf bifurcation to a stable periodic orbit.The medium is unable to synchronize so the phase of the oscillators varies as a function of position.In the present case, the small-amplitude periodic orbit that bifurcates is unstable.Pattern y is time-dependent.It consists primarily of stripes but there are small localized regions that oscillate with a relatively high frequency (- The active regions disappear, but new ones always appear elsewhere.In ).These conditions were then perturbed with f 1% random noise in order to break the square symmetry.The system was then integrated for 200,000 time steps and an image was saved.In all cases, the initial disturbance propagated outward from the central square, leaving patterns in its wake, until the entire grid was affected by the initial square perturbation.The propagation was wave-like, with the leading edge of the perturbation moving with an approximately constant velocity.Depending on the parameter values, it took on the order of 10,000 to 20,000 time steps for the initial perturbation to spread over the entire grid.The propagation velocity of the initial perturbation is thus on the order of 1 x space units per time unit.After the initial period during which the perturbation spread, the system went into an asymptotic state that was either time-independent or time-dependent, depending on the parameter values.
Figures 2 and 3 are phase diagrams; one can view Fig. 3 as a map and Fig. 2 as the key to the map.The 12 patterns illustrated in Fig. 2  Pattern a is time-dependent and consists of fledgling spirals that are constantly colliding and annihilating each other: full spirals never form.Pattern is time-dependent and consists of what is generally called phase turbulence (8), which occurs in the vicinity of a Hopf bifurcation to a stable periodic orbit.The medium is unable to synchronize so the phase of the oscillators varies as a function of position.In the present case, the small-amplitude periodic orbit that bifurcates is unstable.Pattern y is time-dependent.It consists primarily of stripes but there are small localized regions that oscillate with a relatively high frequency (- The active regions disappear, but new ones always appear elsewhere.In Outside the region bounded by the solid line, there is a single spatially uniform state (called the trivial state) (U = 1, V = 0) that is stable for all (F, k).Inside the region bounded by the solid line, there are three spatially uniform steady states.Above the dotted line and below the solid line, the system is bistable: There are two linearly stable steady states in this region.As F is decreased through the dotted line, the nontrivial stable steady state loses stability through Hopf bifurcation.The bifurcating periodic orbit is stable for k < 0.035 and unstable for k > 0.035.No periodic orbits exist for parameter values outside the region bounded by the solid line.

Configuration studied
The configuration studied is a square domain 2.5×2.5 mm 2 , which is meshed with 250×250 fully integrated linear square elements (i.e., with an element size of 0.01×0.01mm 2 ), over which as many user elements are superimposed for the activation of DOFs 12 and 13 (representing the concentrations of  and , respectively) and for the integration of the reaction-diffusion process.A transient 'coupled temp-displacement' procedure is applied.Periodic boundary conditions are prescribed to DOFs 12 and 13 along the border of the domain, as set in [44].The following initial conditions for  and  are defined using a DISP user subroutine: where Ω is a rectangular domain 0.125(1+ )×0.125(1+) with  ∈ [0,1]a random perturbation.Finally,  ! and  !have been set to 10 -5 mm 2 /s and 2.10 -5 mm 2 /s, respectively.Several ,  parameters have been considered, as listed in Table 1.

Results
The Abaqus results for  (NT12) and  (NT13) are presented in Figure 4(a) to 8(a), with the corresponding Python reference results for  shown on (b).All Abaqus computations have been performed with a constant time increment of 10 s, while the python's one is equal to 1 s.It can be observed that our implementation in the Abaqus code is able to reproduce quite well the results obtained with another software, for various configurations.It may be noted that the U-skate geometries exhibited by Munafo [41,46]) could not be generated, as in [44].

Discussion
An important feature observed in our simulations is a non constant velocity of the pattern front, with a strong influence of the ,  parameters.This behavior is consistent with results obtained by other methods, especially in [44].From the Figure 4 to 8, it might be observed that the front velocity computed by Abaqus has the same order of magnitude than the one obtained using Python.
We have also investigated the effects of the element size and of the time increment (see [46] for a more complete investigation of the time increment influence).The influence of the element size is illustrated in Figure 9 for ,  = 0.006,0.037 .When the element size increases, the generated pattern is strongly influenced by the mesh structure and tends to a square rather than a circle.Moreover, the velocity of the pattern front is increased because of a rapidly vanishing  field that annihilates the reaction process.In contrast, decreasing the time increment has no influence on the Abaqus results and on their consistency with [44], except for ,  = 0.022,0.049where the intensities of the pattern oscillations decrease and a steady state is finally reached for  about 3400 s.For this configuration, the influence of the time increment is shown in Figure 10: when it is decreased from 10 s to 1 s, no steady state is reached with Abaqus up to 5000 s and chaotic oscillations are observed, as in [44].

Conclusion
An appropriate application of user elements allows the extension of Abaqus capabilities, including the modification of library elements, the activation of hidden DOFs, and the addition of various physical processes with or without couplings.This study has been focused on the activation of DOFs and on the addition of chemical reactions in Abaqus.An application to the Gray-Scott model has been made successfully.However, this model, though spectacular, has very complex features in term of spatiotemporal evolution, intimately linked with the used parameters.This complexity leads to some difficulties in the definition of the finite element setup in terms of time increment and mesh.Further work will extend the proposed approach to 3D simulations, reactions involving 3 species or more, and mechanical coupling.
To include mechanical interactions, especially, it will be only necessary to introduce in the UELs the related contribution to the weak formulation.Furthermore, a 4 th layer might be added to include the coupling between DOF 11 and mechanical fields.Equation (1) thus becomes (iii) Layer 3 has the same role as layer 2, but for DOF 13. (iv) Layer 4 defines only the relation between DOFs 12 and 13.

Figure 1 .
Figure 1.Principle of the implementation of a multidiffusion process.
FIG. 11.Breathing small F-state domain.Period of oscillations 24 min.Conditions as in Fig. 10.Sampling time of snapshots from left to right and down 4 min.Size of the frames 3.5ϫ 3.5 mm.

Figure 3 .
Figure 3. (a) Types of patterns obtained with the GS reaction, and (b) their position in the (, ) plane (using != 2 != 2×10 !!) as defined in[36].For (, ) points where no pattern is specified, a constant homogeneous field for  as well as for  is expected.
d 3 are phase diagrams; one as a map and Fig. 2 as the key he 12 patterns illustrated in nated by Greek letters.The the concentration of U with g U = 1 and blue representellow is intermediate to red ig.3, the Greek characters ttern found at that point in parameter space.There are two additional symbols in Fig. 3, R and B, indicating spatially uniform red and blue states, respectively.The red state corresponds to (U = l,V = 0) and the blue state depends on the exact parameter values but corresponds roughly to (U = 0.3,V = 0.25).Pattern a is time-dependent and consists of fledgling spirals that are constantly colliding and annihilating each other: full spirals never form.Pattern is time-dependent and consists of what is generally called phase turbulence

Flg. 3 .
-. .... , ,I ___...--._..-- k agram of the reaction kinetics.ion bounded by the solid line, spatially uniform state (called (U = 1, V = 0) that is stable for he region bounded by the solid three spatially uniform steady he dotted line and below the stem is bistable: There are two eady states in this region.As F rough the dotted line, the nondy state loses stability through .The bifurcating periodic orbit 0.035 and unstable for k > dic orbits exist for parameter e region bounded by the solid .Ig. 2. The key to the map.The patterns shown in the figure are designated by Greek letters, which are used in Fig. 3 to indicate the pattern found at a given point in parameter space.The map.The Greek letters indicate the location in parameter space where the patterns in Fig. 2 were found; B and R indicate that the system evolved to uniform blue 0.06 and red states, respectively.SCIENCE VOL.261 9 JULY 1993 about the center of the grid was then perturbed to (U = ln,V = 114 are designated by Greek letters.The color indicates the concentration of U with red representing U = 1 and blue representing U = 0.2; yellow is intermediate to red and blue.In Fig. 3, the Greek characters indicate the pattern found at that point in parameter space.There are two additional symbols in Fig. 3, R and B, indicating spatially uniform red and blue states, respectively.The red state corresponds to (U = l,V = 0) and the blue state depends on the exact parameter values but corresponds roughly to (U = 0.3,V = 0.25).

Fig. 1 .
Fig. 1.Phase diagram of the reaction kinetics.Outside the region bounded by the solid line, there is a single spatially uniform state (called the trivial state) (U = 1, V = 0) that is stable for all (F, k).Inside the region bounded by the solid line, there are three spatially uniform steady states.Above the dotted line and below the solid line, the system is bistable: There are two linearly stable steady states in this region.As F is decreased through the dotted line, the nontrivial stable steady state loses stability through Hopf bifurcation.The bifurcating periodic orbit

2 .
The key to the map.The patterns shown in the figure are designated by Greek letters, which are used in Fig.3to indicate the pattern found at a given point in parameter space.Flg.3. The map.The Greek letters indicate the location in parameter space where the patterns in Fig.2were found; B and R indicate that the system evolved to uniform blue 0.06 and red states, respectively.

Figure 9 .
Figure 9. Influence of the element size on the Abaqus results at  = 800  for u (left) and v (right), with F, k = 0.006,0.037 .

Figure 10 .
Figure 10.Influence of the time increment on the Abaqus results at  = 5000  for  (left) and  (right), with F, k = 0.022,0.049 .

( 5 )
carried out within the framework of the French Federation for Magnetic Fusion Studies (FR-FCM) and of the Eurofusion consortium, and has received funding from the Euratom research and training program 2014-2018 and 2019-2020 under grant agreement No 633053.The views and opinions expressed herein do not necessarily reflect those of the European Commission.