Effect of matrix solidification on the structure formation in electromagnetic suspensions

For suspensions with electromagnetic particles exposed to an external field, we examine the effect of the solidification of the suspending medium on the formation of particle structures, representative of the curing of a photo-reactive resin during stereolithography. To that end, Brownian Dynamics (BD) simulations are examined in which the solidification of the suspending medium can be incorporated by increasing its viscosity in the course of time. For illustrative purposes, it is assumed that the viscosity function is known apriori in explicit and parametrized form, however, experimental data can be used as well. It is demonstrated that one can study the effects of the viscosity increase by a transformation of time, akin to the time–temperature superposition principle, but here also in the presence of thermal noise on the suspended particles. Therefore, instead of performing BD simulations with a continuously increasing viscosity (computationally inefficient), we advocate performing simulations at a constant (low) viscosity and subsequently transform time (nonlinearly) for re-interpretation of the simulation results. So doing, one can predict the formation of particle structures during on-going solidification of the suspending medium. In practice, the viscosity increase is so drastic that further evolution of the particle structure can be considered as arrested after the characteristic transition-time of the viscosity is reached. Semi-quantitative rules of thumb are formulated for the 3D-printing practitioner.


Introduction
Additive manufacturing is a topic of scientific interest for a long time [1,2]. Various materials are used including metals [3], polymers [4], ceramics [5] and combinations of them [2]. In this paper, we are interested in polymer-matrix composite materials used in stereolithography (SLA) [6,7] of UV-curable resins. Polymer-based materials are printed for various purposes such as rapid prototyping [8], or topology optimization [9]. Applications include aerospace (cabin interiors) [10], medicine (mimics of living tissue) [11], anthropology, where medieval skulls are reconstructed [12], and design (spatially-dependent elasticity) [13]. The use of composites in SLA is possible as multi-materials including suspensions are already being printed [14], however, challenges are still present [15]. The induction of dielectric, conductive and/or magnetic properties to the printed object [8,[16][17][18] has been investigated by the usage of electromagnetically active particles. The internal structure of the composite can be tailored by manipulating, while printing, the particles with an electric or magnetic field. These systems can be used in applications like personalised hearing aids [19], lenses being flat with a micro-structure consisting of a gradient concentration of particles that bestow the functionality of their curved counterparts [20], piezoelectric or Hall effect sensors [21], and direction-specific thermally or electrically conductive composites [22]. The formation of structures of particles in external fields is well documented in the literature, e.g., the formation of strings in field direction [23,24], as well as the merging/thickening and percolation of particle chains [25,26].
The behaviour of such systems is highly influenced by the suspending medium. When an external field is used, the particles create structures [27]. The behavior during structure formation is mostly investigated in a liquid suspending medium, as the structure formation is accessible under experimental conditions [28,29]. However, in SLA the medium is irradiated and as a result it solidifies. This process could be used to fixate the structure, so that the formed structure persists even in the absence of an imposed field. Photo-reactive resins are used in SLA [30], and also as suspending medium [31,32]. The polymerization of such resins is initiated by energy input from a UV-source, and solidification can be envisioned as an increasing viscosity, which results in a slowing-down of the structure formation and eventually might even completely arrest the evolution, whereby the particle arrangement is fixated. The formation of (super-)structures becomes feasible within the usual time-scales ( O(1s) ) of 3D-printing due to the resins having a low viscosity before curing. In this paper, the effect of the solidification of the suspending medium on the structure formation is investigated.
To this end, Brownian Dynamics (BD) simulations are used with dipole-dipole interactions, to simulate the motion of particles (dipoles) in a fluid under an external field. This method provides the evolution of the particle arrangement, and the structure created is characterised. One can achieve various structures of particles depending on the particles, suspending medium and field conditions, e.g. strings [33], planes [34], and networks [22]. The main features of the structure can be identified and quantified. Structures like strings can be described by their length, orientation, and thickness, while planes by their size and orientation, and networks by the number density of branch-points (BP), the degree of BP, the thickness of the branches, the number of clusters present, and the existence or lack of percolation. Various techniques have been used for the characterisation of a structure of particles or points, like the characteristics of the gyration tensor [35], or the Voronoi tessellation technique [36,37]. Voronoi polyhedra have proven useful in the characterisation of many systems, e.g. disordered systems [38], polymers [39][40][41], colloidal gels [42], and granular materials [43]. The (radial or cylindrical) pair-correlation function, bond-angle distribution [44], and quermass integrals [45] are alternative methods for the characterisation of the structure of particles. In previous work, we presented a combination of graph theory tools with the image analysis technique called skeletonization [46][47][48] for the quantitative characterisation of the structures, particularly networks. These, together with the previously presented cluster measures [49,50] for orientation, anisotropy and size, are used to illustrate the effect of solidification on structure formation.
The method developed in this paper is closely related to the well-known time-temperature superposition principle [51]. In our case, the viscosity is employed to transform between different representations of time. In this way, processes on time scales that are accessible numerically/experimentally can be related to processes on much longer time scales when matrix solidification slows down the dynamics dramatically.
The goal of this paper is to investigate the effect of the solidification of the suspending medium on the structure formation, so that the most relevant parameters of a prototypical viscosity profile are identified and semi-quantitative rules for the 3D-printing practitioner can be formulated. To this end, in Sect. 2 we prove analytically with two different procedures that the effect of a increasing viscosity over time on the structure formation of particles can be studied by means of time transformation, even in the presence of thermal fluctuations. In Sect. 3, a viscosity profile representative of the solidification of a photo-reactive resin [52] is introduced. In addition, the simulation details for suspensions of particles are presented, together with our morphology measures for evaluating the structure formation over time. Previously obtained results under constant (low) viscosity [49] are used as an illustrative example of transforming time such that the solidification is apparent. The paper is summarised and the results are discussed in Sect. 4.

System and dimensionless quantities
In this section, the methodology is discussed, including the simulation details, and morphology measures. The system consists of particles with radius R p at volume fraction of the suspension, and the monomeric fluid that serves as suspending medium with viscosity (t) , where t is time; the particles are contained in a finite box of edge-length L, where periodic boundary conditions are applied in all directions. BD simulations are used to simulate the time evolution of the system. For the rest of this paper, we present quantities in dimensionless terms: the positions of the particles are scaled with the characteristic length-scale, r c [49], and time is scaled with the characteristic time-scale, t c [49]. The concrete expressions for r c and t c will be specified in Subsect. 3.2. In Subsect. 2.2, a viscosity profile of the form will be used, where c is a characteristic value for the (dimensional) viscosity . It will be shown in the following, that the solidification of a resin, i.e., of the suspending medium, described by an appropriate time-dependent viscosity can be studied without the need for additional simulations. The viscosity function used is assumed to not depend on the particle positions. Any positive and monotonically increasing function (t) can be used to this end, if it is turned to dimensionless form, see Eq. (1).

Nonlinear transformation of time
When one refers to BD, one considers a stochastic differential equation (SDE) describing the dynamics of the microscopic state of the system specified by the quantities . The system can be described also by the corresponding Fokker-Planck equation (FPE), which describes the evolution of the distribution function p( ) of the microscopic state of the system. These two approaches are equivalent [53].

Procedure 1: Stochastic Differential Equation
If one considers a SDE of the form where t is the time with a varying viscosity, t is the state at time t , ̂ t and ̂ t are functions depending on the state, d t are increments of Wiener processes representative of white (i.e., uncorrelated) noise, and * (t) is the timedependent viscosity. In the following, it is assumed that the viscosity function is known, does not depend on t , is positive, and is increasing monotonically with time. Equation (2) can be simplified by introducing where s is the rescaled time. If Eq. (2) was an ordinary differential equation, the transformation of time would be straightforward. However, the stochastic term, i.e. the Wiener processes, complicates the procedure. When one rescales time, one should consider the condition of Wiener processes being uncorrelated in time, where ⟨…⟩ denotes the statistical average over all realizations of noise, (…) is Dirac's -function, is the unity tensor, and t and t ′ are two moments in time. If one performs the time rescaling of Eq. (5) with the help of Eq. (3), and uses If one rescales Eq. (2), with the help of Eqs. (3) and (6), one finds . Equation (7) resulted from the transformation of time from a time-dependent viscosity to a constant viscosity. Therefore, Eq. (7) is equivalent to BD simulations in dimensionless form like the ones presented in Ref. [49] for constant viscosity. One can simulate a system of varying viscosity (e.g. solidifying resin, see Eq. (2)) by means of nonlinear time-mapping of constant-viscosity simulations (described by Eq. (7)).

Procedure 2: Fokker-Planck Equation
In the previous procedure, we elaborated the procedure on the SDE; here, the same conclusion is derived on the basis of the FPE. According to [53], the SDE Eq. (2) can also be represented in the form of the following FPE for the distribution function p( ), where (…) T denotes the transpose. Since * is independent of , multiplication of Eq. (8) by * is straightforward and leads to where p∕ s = (dt∕ds)( p∕ t) = * p∕ t has been used. Eq. (9) is compatible with the SDE of Eq. (7).

Application of the general procedure
In Sect. 2, we have presented two procedures with which the time-dependence of the viscosity function can be rationalized conveniently by transforming the time. In simpler terms, one can study the dynamics of particles during solidification of the medium (increasing viscosity over time) by transforming the time accordingly. In this section, a concrete example of a prototypical viscosity function is used to illustrate this principle. To this end, results obtained through constant viscosity simulations are re-interpreted as if simulations of a solidifying medium were performed. The prototypical viscosity, the simulation details for the results, the measures for the structure quantification and the transformed results are presented and semi-quantitative rules are formulated. In general, it can be anticipated that the formation of particle structures induced by the imposed field is slowed down by the increase of the viscosity with time; at drastically increased viscosity, structure evolution may even be completely arrested, i.e., the particle arrangement is frozen in. The goal of this section is to elaborate on these general thoughts in more detail.

Time-dependent viscosity during solidification
For transforming time, one needs a viscosity function representative of the curing process. Experimental observations [28,29] indicate that this function has an initial and final plateau, and a transition of certain width in between. The initial plateau value corresponds to the viscosity of the monomeric resin, and the final plateau to the cured sample [52]. In order to capture these main characteristic features in a form suitable for a parameter study (rather than giving a quantitative representation of a certain experimental data set), we choose a hyperbolic tangent of the form where * = * 0 + * ∞ ∕2 and Δ * = * ∞ − * 0 ≥ 0 , with * 0 and * ∞ the values of the initial and final plateaus, respectively, t 0 indicates the mid-point of the transition, and represents the width of the transition, see Fig. 1. The entire viscosity function has already been scaled with c , see Eq.
(1). For the rest of this paper, * 0 = 1 will be used, meaning that the dimensional viscosity at t ∼ −∞ is equal to c .
If one uses Eq. (3) for the transformation of time in combination with Eq. (10) for the explicit viscosity, one obtains The ratio t − t 0 ∕ is increasing with increasing t , so one can anticipate that the second contribution in the numerator of the fraction in the logarithm will cause numerical problems at late times. However, precisely in the limit where this contribution gives numerical problems, the first term in that numerator can be neglected, which results in the following approximation of the s(t)-relation, However, Eq. (12) is not valid for the whole range of t . We choose to use a combination of Eqs. (11) and (12): Eq. (12) is used for large t if the relative error s − s app ∕s is smaller than 10 −4 . In Fig. 2, one can see the relative error plotted against the ratio t − t 0 ∕ .
Our intention is to transform the dynamics from the s -representation to the t-representation, i.e. from an existing constant-viscosity simulation to a time-dependent viscosity situation. Therefore, one needs to determine t from s , which means solving Eqs. (11) and (12), using s-values corresponding to constant-viscosity simulations. In the case of Eq. (12), one can solve towards t analytically, however, for Eq. (11) Fig. 1 Viscosity function based on a shifted hyperbolic tangent with parameters * 0 , * ∞ , t 0 , and this is non-trivial, and the Newton-Raphson method [54] is thus used for solving Eq. (11) numerically. In the following, the dependence of the t − s relation on the parameters of the viscosity function, see Eq. (10), is explored. The parameters to be varied are * ∞ , t 0 , and . The results can be found in Fig. 3. First, we choose the values for the study of the * ∞ -variation, Fig. 3a, such that they are representative of the real process, namely * ∞ = 10 2 , 10 5 , 10 8 , which correspond to the values obtained for different light intensities [52]. The values for the other parameters are = 0.1 and t 0 = 3.0 . Afterwards, we choose the values for the variation in transition time t 0 with respect to the crossing point of G ′ and G ′′ observed experimentally ( t 0 ∼ 2.5 ) [52]; the values are t 0 = 1.0, 3.0, and 10.0 , and the values for the two additional parameters are = 0.1 and * ∞ = 10 5 . The transition time has limitations in our case, as we attempt to use a compact mathematical function to represent the real viscosity function. The limitation has to do with the existence of the initial plateau depicted in Fig. 1; if one uses a low t 0 -value and high -value, then at time t = 0 the transition will be taking place already. This effect is unphysical and for this reason we choose t 0 such that at time t = 0 the viscosity has a value as close as possible to 1.0. The values for studying the variation of are = 0.1 , 0.5, and 1.0, and the values for the other parameters of the function are * ∞ = 10 5 and t 0 = 5.0.
In Fig. 3, one can observe the relation between the time t and the time s for a viscosity function given by Eq. (10). In  large. In Fig. 3b, one can see the dependence on t 0 , which defines the time at which the behaviour deviates from the no-transition case t = s . In the inset of Fig. 3b, one sees the steepness of the transition, which is mainly affected by the transition time , which is low, and the viscosity difference * ∞ − * 0 , which is large. In Fig. 3c, one can see the dependence on . This parameter controls the duration of the transition in the * -function, however, here it seems to affect also the deviation from the no-transition case t = s , as the transition starts earlier if the is larger. The fact that = 1.0 starts from a higher value of t is based on not being significantly lower than t 0 , and the viscosity has a value quite higher than * 0 already at s = 0 . One can observe that for = 1.0 , the transition exhibits a smoother behaviour (smoother increase) compared to the the other values, as it lasts longer.

Simulation details
Considering the many-particle system, one uses = ( 1 , 2 , ..., n ) , where i is the center position of particle i. These particles have induced dipoles due to the mismatch in the dielectric permittivity of the matrix, m , and the particles, p , under an external field. The interactions governing the structure formation are the dipole-dipole interactions [27,55]. Here, BD simulations [53,56] are used in dimensionless form; including the nonlinear re-mapping of time in relation to the time-dependent viscosity-introduced above-one obtains (see Ref. [49] for the t-representation) with uncorrelated and normal-distributed random vectors i , and where em * i describes the dipole-dipole interaction, exv * i accounts for the excluded volume of the particles, * i and p * i are the scaled dipole moment of particle i and its magnitude, r ij and ̂ ij are the length and unit vector, respectively, of ij = i − j , with ̂ ij = ij ∕r ij , and R p denotes the radius of the particles. The parameter = 30 is a constant that determines the range of the excluded volume interaction, while the pre-exponential factor defines the overall strength of the interaction [57]. The Boltzmann constant is denoted by k B , the temperature by T, and the friction coefficient by = 6 c R p , with c the value of viscosity. The positions of the particles are scaled with the characteristic length-scale, r c = 1∕ 3 √ n d for given number density of particles n d , and the dipole moments with the characteristic dipole moment, p c = 4 m KR 3 p E c with E c the constant value of the field strength and K = ( p − m )∕( p + 2 m ) the Clausius-Mossotti constant [58]. The time in the t-representation is scaled with the characteristic time-scale introduced in ref. [49], t c = 8 2 m c R p r 5 c ∕p 2 c . It is mentioned that the characteristic time scale t c (for two particles to come into contact) changes as the viscosity changes in the course of time; therefore, increasing viscosity will result in an overall delay of the structure formation, with potentially severe consequences (e.g., complete arrest of the structure formation) if the viscosity increases by several orders of magnitude.

Morphology characterisation
In this subsection, we briefly summarize the measures for morphology characterisation that have been introduced in our earlier work [49] and that we are going to use for the rest of this paper. The first measure for clusters/structures concerns the relative orientation of the direction of the interparticle vectors with another direction, in that case the direction ̂ of the external field, where the the cluster index is denoted by I, the average being taken over all pairs of particles that belong to the same cluster I. Only the primary image is considered, as periodic boundary conditions are applied.
The second measure quantifies the anisotropy of each cluster. For that purpose, we calculate the gyration tensor [35], where cm,I is the center of mass of the particles of cluster I. An ellipsoid is defined by the gyration tensor [35] and the relative eigenvalues give the relative discrepancy in size of the three axes of the ellipsoid. We introduce a measure that quantifies the anisotropy, where max is the maximum eigenvalue of the gyration tensor, and det I is the determinant of the gyration tensor.
The third measure for clusters has to do with the size of the cluster, where N I is the number of particles of cluster I, and N L is the minimum amount of particles that are needed to be stacked together for spanning the box-length in one direction.
The averages over all clusters of the measures S 2,I , * I , and N * I are reported in the following, and these averages will be denoted by same symbol, but without the cluster index. The number of particles per cluster serves as a weight for these averages, where A is any quantity of the ones referring to clusters, N cl is the number of clusters, N I is the number of particles in cluster I, and A I is the value of the quantity A for the cluster I; normalization is achieved by N p = N cl N I .

Results for unidirectional constant field
As an illustrative example, the nonlinear time-mapping will be applied here to identify the most important parameters for structure formation in a solidifying resin. The results concern suspensions of particles under a unidirectional and spatially constant field [49]. We use the measures introduced in Subsect. 3.3 for tracking the structure formation. The study of the parameters in the viscosity function (10) allows for the construction of rules of thumb concerning their effect on the structure formation.
The results of BD simulations with (low) constant viscosity obtained in previous work [49] are interpreted here as s -representation results, and transformed in order to obtain the structure formation with a prototypical time-dependent viscosity, i.e. the t-representation. The system consists of 10 3 particles at volume fraction = 10% , and three simulations with different random numbers sequences, and different random initial configurations are used to generate the error-bars [49]. In all simulations, B * and t c are kept constant. For typical snapshots of the particle structures obtained, the reader is referred to our earlier work [50], Fig. 3 and Fig. 15 therein.
We use the results obtained in Ref. [49], which are constant-viscosity simulations representative of s-dynamics, and re-interpret them by transforming to the t-representation for In Fig. 4, one can observe the three measures versus t for different values of * ∞ , namely * ∞ = 1.0 , 10 2 , 10 5 , 10 8 , with t 0 = 3.0 and = 0.1 . It is clear that when the transition happens the time is extended by some orders of magnitude with respect to the constant viscosity case. This deceleration depends strongly on the * ∞ -value, since it affects the effective freezing of the structure formation.
In Fig. 5, the influence of the transition time t 0 is explored on the three measures for the morphology for t 0 = 1.0 , 3.0, 10.0, with * ∞ = 10 5 and = 0.1 . The transition time seems to be a rather important parameter, since after this timevalue the structure formation is slowed down significantly (orders of magnitudes). The difference with the * ∞ -variation is that the amount of deceleration of the dynamics is constant (constant * ∞ ), while t 0 controls the moment when this deceleration will occur.
The variation of the parameter can be seen in Fig. 6, for values = 0.1 , 0.5, and 1.0, with t 0 = 5.0 and * ∞ = 10 5 . This parameter determines how swiftly the deceleration occurs. One has to note that this parameter is not affecting how much the dynamics are decelerated or when it happens, however, it affects the way the structure formation is decelerated before and after the transition time.
The data shown in Figs. 4, 5, and 6 is available via a repository [59].

Summary and discussion
In this paper, it has been shown that a suspension of particles in a solidifying matrix can be studied by means of a nonlinear transformation of time of constant-viscosity dynamics. The only assumption is that the solidification is described by a known (monotonically increasing) viscosity-function. The described procedure is equivalent to the time-temperature superposition principle [51], however, in our case in the presence of thermal noise. Our method is suitable for predicting the behavior under experimentally unfeasible or computationally inefficient conditions (i.e., on prohibitively long timescales), by transforming time nonlinearly from results obtained under accessible conditions. While the 3D-printing of composites served in this paper as the initial motivation for developing the rescaling of time, it is to be seen as a prototypical example for a wider class of systems where characteristic time-scales change drastically in the course of time and where stochastic dynamics are involved.
Instead of studying the constant-viscosity scenario followed by a subsequent transformation of time, one might choose to examine the case of t-dependent viscosity directly, e.g., in numerical simulations. If one chose to do so with a constant time-step for the entire simulation, that timestep would have to be chosen small enough to capture the rapid dynamics at early time (low viscosity) well, whereas towards the end a lot of (time-) iteration steps would be done when basically nothing happens because the viscosity has increased by many orders of magnitude and the dynamics is slowed down drastically. As an alternative to a constant time-step, one could adapt the size of the time step in proportion to the increasing viscosity, which is analogous in spirit to the rescaling of time we use in this paper, see Eq.
(3). In Sect. 3, a viscosity function representative of the experimentally observed solidification of a photo-reactive resin [52] has been used for re-interpreting constantviscosity simulations of suspensions as simulations with an increasing viscosity. The viscosity function has been assumed to have the form of a shifted hyperbolic tangent. The relation between the two scalings of time, t and s , has been investigated with respect to variations of the parameters of the viscosity function, * ∞ , t 0 and . The functional form used for the viscosity allowed us to examine, by way of a parameter study, the consequences of the characteristic features on the structure formation. Obviously, if experimental data for the viscosity are available that are not represented well by the time-dependence described in Eq. (10), one may use the experimentally obtained function (t) for the transformation of time as described by the relations Eq. (3) and Eq. (4), respectively.
In this entire paper, it is assumed that the temperature T is constant, because our interest is solely in the effect of solidification, due to a chemical/polymerization reaction. For clarification purposes, and in view of stimulating future work, the following comment about the influence of temperature is added here. In the rescaled evolution equations (13), the temperature T enters explicitly only in the parameter B ⋆ (see Eq. (16)). In addition, the temperature T enters, in general, implicitly via the viscosity . Typically for polymerbased systems, varying the temperature by a few (tens of) degrees has a significant effect on the viscosity. Therefore, transformation of time is significantly more relevant for the T-related changes in the viscosity rather than the changes in B ⋆ . Therefore, if changes in temperature are of interest, it is recommended to tackle only the viscosity-effects in the transformation of time, while the changes in B ⋆ can be left as is in Eq. (13) with Eq. (16).
The competition between the structure formation and the curing kinetics realised by an increasing viscosity has given rise to some rules of thumb for the 3D-printing practitioner. The effect of * ∞ is found to alter how much the dynamics are decelerated; the effective freezing of the structure can be controlled. The parameter t 0 has the most important effect, to the best of our knowledge, as it defines the moment when the deceleration occurs; it can be tuned such that the desired