INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 40, 2679–2692 (1997) A SPACE–TIME FINITE ELEMENT MODEL TO STUDY THE INFLUENCE OF INTERFACIAL KINETICS AND DIFFUSION ON CRYSTALLIZATION KINETICS K. M. KIT∗ AND J. M. SCHULTZ Materials Science Program; University of Delaware; Newark; DE 19716; U.S.A. SUMMARY We present a space–time nite element formulation to study the cooperative growth of adjacent needlelike crystals in a two-dimensional, binary melt. It is assumed that the system is isothermal and that the compositions of the melt and the crystals are dierent. The growth rate of the crystals is taken to be a function of the melt composition in front of the growing crystals, and the composition of the melt as a function of space and time is determined by the diusion equation. The positions of the growth fronts of each crystal are tracked. Good agreement is found between the numerical solution of an approximated onedimensional problem and an analytical solution. Numerical results of the simulation of the growth of isolated and adjacent crystals are presented. ? 1997 by John Wiley & Sons, Ltd. KEY WORDS: space–time nite elements; phase transformation; crystal growth; polymer blend; numerical simulation INTRODUCTION The rate at which a material crystallizes in a supercooled melt is inuenced by interfacial kinetics, diusional kinetics (of heat or mass), or both. As growth occurs, latent heat is produced which must diuse away from the growth front. If the composition of the growing crystal and the melt dier, then rejected mass must also diuse away from the growth front. The rates of both heat production and mass ux are dependent on the velocity of growth. The growth velocity may be dependent on temperature and composition, which, in turn, are aected by the rates of heat production and solute rejection. For many reasons, it is desirable to be able to predict the rate at which crystallization occurs given the conditions of the system. Problems which are diusion limited (e.g. the Stephan problem) have been extensively studied. The one-dimensional diusion limited √ problem usually yields a solution for the position of the growth front, X , of the type X ∝ Dt, where D is a diusion coecient and t is time. Problems which are dominated by interfacial kinetics usually yield a solution of the type X ∝ t. Linear growth rates are very commonly observed during the quiescent crystallization of polymers. Traditionally, the Finite Element (FE) analysis of time-dependent problems has employed a nite element discretization of space and a nite dierence discretization of time.1 However, a space– time FE method in which time is discretized similarly to the spatial variables can be advantageous. This type of method can easily allow for deformation of the spatial domain with time. Bonnerot ∗ Present Address: Department of Materials Science and Engineering, 434 Dougherty Engineering Building, The University of Tennessee, Knoxville, TN 37996-2200, U.S.A. CCC 0029–5981/97/142679–14$17.50 ? 1997 by John Wiley & Sons, Ltd. Received 6 February 1996 Revised 12 November 1996 2680 K. M. KIT AND J. M. SCHULTZ and Jamet 2; 3 have used such a method to simulate the one- and two-dimensional Stephan problems. Weill et al.4 used a similar method to model the temperature eld and freezing front in a domain containing two cryoprobes (heat sinks). These problems are diusion controlled and interfacial kinetics are not incorporated into the models. The space–time formulation can also be useful for other applications. Tezduyar et al.5 have developed a space–time FE model to study uid ow around a moving object. The purpose of this paper is to present a nite element model which can be used to study the eects of both interfacial and diusional kinetics on the crystallization rate and morphological development of polymer spherulites. We present a space–time formulation based on the diusion equation to track the isothermal growth of crystals from a melt of a dierent composition. In this formulation, the position of the interface varies continuously in time. We will assume a functional relationship between growth rate and the composition of the melt at the interface. The crystal is allowed to grow at a constant rate for a certain time interval. This growth causes exclusion of solute from the crystal, changing the composition of the melt at the interface and, hence, the growth rate for the next time interval. In this way, we can track the growth of one isolated crystal or the cooperative growth of many adjacent crystals. This formulation can be used to study diusion-limited and diusion-coupled growth. In the latter, growth conditions can reach a steady state, but these conditions are dependent on diusional and interfacial kinetics. DESCRIPTION OF THE PROBLEM We consider a group of adjacent polymer lamellae growing from an innite polymer melt. The following assumptions are made. (i) The system is isothermal during growth. (ii) The group is dened by the following parameters; N , number of lamellae in the group, L, long spacing (or distance between lamellae centres) Lc , lamellae thickness, and S, a vector of initial lengths of the lamellae. Further, L and Lc are assumed to be constant during growth and to vary with temperature. (iii) The melt consists of a miscible binary blend of polymers A and B. (iv) The lamellae consist only of polymer A. (v) The growth velocity, for a particular A=B blend, is a function of the temperature and melt composition adjacent to the growing crystal front. (vi) The diusion coecient, D, is constant during growth and is a function of temperature. As a lamella grows, B is typically excluded from the crystal and it accumulates in front of the lamella. The extent to which B is excluded is given by the distribution coecient, k, dened as the ratio of CS∗ to CL∗ , where CS∗ and CL∗ are the concentration of B in the solid and liquid, respectively, adjacent to the interface. Total exclusion occurs when k = 0: (This is the case for the present problem.) Total or partial exclusion results in a compositional gradient and diusion of B chains away from the growing interface. Consider a group of lamellae as depicted in Figure 1. Growth occurs in the x and z directions only. (The y dimension of the lamellae is given by Lc .) Generally, the dimensions of the lamellae in the z and x directions quickly grow much larger than Lc . Away from the corners then, diusion of B chains is eectively limited to 2 dimensions; the y direction and the direction normal to the growth face. Therefore, the full three-dimensional diffusion problem can be closely approximated by a two-dimensional problem, as shown in Figure 2. Here the lamellae grow in the +x and −x directions only. It is assumed that all of the lamellae nucleated from the same x position. Assuming that each lamella grows equally in the +x and −x directions, the problem becomes symmetric about the dotted line ab. We will also assume that cd Int. J. Numer. Meth. Engng., 40, 2679–2692 (1997) ? 1997 by John Wiley & Sons, Ltd. A SPACE–TIME FINITE ELEMENT MODEL 2681 Figure 1. 3-D lamellae stack Figure 2. 2-D domain showing lines of symmetry is a line of symmetry. The reduced spatial solution domain is represented by the cross-hatched region in Figure 3 and denoted as . The boundary, , is partitioned into 1 , 2 , and 3 . 1 is the far eld boundary. 2 is a symmetry or insulating boundary which includes the lamellae surfaces normal to the y direction as well as the surfaces labelled in Figure 3. 3 contains the growth surfaces. The regions occupied by lamellae are not part of the solution domain since it is assumed that CB = 0 there. It should be noted that we are interested here in centre-of-mass diusion of the polymer chains. The order of magnitude of the radius of gyration of the polymer coils is generally 10 nm while the order of magnitude of the diusion length is generally not less than 100 nm. Thus, the diusion length is at worst one order of magnitude larger than the molecular diameter, and macroscopic diusion laws should provide an excellent approximation to reality. The mathematical problem to be solved is the time-dependent diusion equation over the spatial domain shown in Figure 3, with appropriate initial and boundary conditions: @CB − D∇2 CB = 0 in G @t 0 ∇CB ·˜n = 0 on 1 ? 1997 by John Wiley & Sons, Ltd. (1) (2) Int. J. Numer. Meth. Engng., 40, 2679–2692 (1997) 2682 K. M. KIT AND J. M. SCHULTZ Figure 3. Reduced x; y domain ∇CB ·˜n = 0 on 2 (3) (1 − k)CB V − D∇CB ·˜n = 0 on 3 (4) CB (x0 ; y0 ; t 0 = 0) = C0 (5) where D is the mass diusion coecient, ∇ is the two-dimensional gradient operator, CB is the volume fraction of B chains, V is the growth velocity, ˜n is an outward normal unit vector, and G is the (x0 ; y0 ; t 0 ) solution space. Equation (4) is a mass balance between the rate at which solute is excluded from a growing crystal and the rate at which the solute diuses away from the growth front. Equation (5) prescribes an initially homogeneous melt composition. Now we wish to non-dimensionalize this problem by dividing each independent variable by a characteristic parameter. For the x direction, D=V0 is chosen as a length scale where V0 is a typical growth velocity of the given system. D=V0 is commonly known as the diusion length and is a rst approximation for the depth of a diusive boundary layer in front of a growing object. For the y direction, Lc is chosen as a length scale. Lc is a measure of the distance that a B chain must diuse in the y direction to get out of the way of the growing crystal. A characteristic time is chosen to be D=V02 . The non-dimensional variables are dened as follows: x= x0 ; D=V0 y= y0 ; Lc t= t0 D=V02 We will dene C = CB − C0 , where C0 is the initial concentration of B chains in the melt and C ∗ to be the value at the growth front. The reduced velocity is dened as Vr = V=V0 where V0 is the velocity associated with C ∗ = 0 and V is a velocity associated with a general value of C ∗ . In most blend systems is possible to write V = V0 f(C ∗ ) or Vr = f(C ∗ ). Equations (1) – (5) can now be written as @(C) 2 − ∇˜ (C) = 0 @t ˜ ∇(C) ·˜n = 0 Int. J. Numer. Meth. Engng., 40, 2679–2692 (1997) in G (6) on (7) 1 ? 1997 by John Wiley & Sons, Ltd. 2683 A SPACE–TIME FINITE ELEMENT MODEL ˜ ∇(C) ·˜n = 0 on 2 (8) @(C) =0 @x on 3 (9) (1 − k)f(C ∗ )(C ∗ + C0 ) + C(x; y; t = 0) = 0 (10) where ∇˜ is the operator, ((@=@x) î + De (@=@y) ĵ); and De (Deborah number) is the dimensionless quantity, D=Lc V0 . The Deborah number can be thought of as the ratio of a characteristic time for crystallization, Lc =V0 , to a characteristic time for diusion in a direction perpendicular to the growth direction, L2c =D. It is noted that the non-dimensional problem as stated above is uniquely determined by the following parameters; k; C0 ; De; L; N; f(C ∗ ); and S. The actual physical problem is dened over an innite spatial domain, and the true far eld boundary condition is C = 0 as x or y → ∞. Here, however the innite domain is approximated by a nite one. The far eld boundary condition must then also be approximated on the articial boundary, 1 . This can be done by ensuring that both C and ∇(C) · ˜n are zero along 1 . However, only one of these conditions can be specied on the boundary. Satisfaction of the other condition must be conrmed after a solution is found. The articial boundary must be suciently far from the growing lamellae for it to have no eect on the diusion eld and for the unspecied condition to be satised. In this work, the specied condition will be the zero ux condition. This is chosen because it leads to a simpler weak form as shown below and because the C = 0 condition is somewhat easier to conrm a posteriori. The initial spatial domain at t = 0 is given by 0 and bounded by . As time progresses, each lamella grows at a rate V which is a priori an unknown function of time. Therefore, the 3 boundaries move with time. 1x will also move with time in order to maintain a minimum distance away from 3 . All other portions of are stationary. The spatial domain (t) will trace a volume G and its boundary, , will trace a surface S in x; y; t space. FINITE ELEMENT FORMULATION A weak form of the dierential equation is formed by multiplying equation (6) by an arbitrary function, w, and setting the integral of the product to zero. ZZZ w G @(C) 2 ˜ − ∇ (C) dG = 0 @t (11) After integration by parts (see the appendix), the weak form becomes ZZ ZZ ZZZ t=t2 @w dG + w(C) d − (C) wCf(C) dS @t t=t1 G S3 ZZ − @Xmax dS + w(C) @t S1x ZZ − S2 ZZ ZZZ G ˜ · ∇(C) ˜ ∇w dG − S1 ˜ w∇(C) ·˜n dS ZZ ˜ w∇(C) ·˜n dS − ? 1997 by John Wiley & Sons, Ltd. S3 ˜ w∇(C) ·˜n dS = 0 (12) Int. J. Numer. Meth. Engng., 40, 2679–2692 (1997) 2684 K. M. KIT AND J. M. SCHULTZ Figure 4. An example 8-node element in x; y; t space Due to the boundary condition on S1 and S2 that ∇CB ·˜n = 0, the sixth and seventh terms above vanish. The third and eighth terms can be combined and the resultant weak form is ZZ ZZZ t=t2 @w dG (C)w d − (C) @t t=t 1 G ZZ ZZZ @Xmax ˜ · ∇(C) ˜ ∇w dG − dS w(C) + @t G S1 x ZZ + w Cf(C ∗ ) − D∇CB ·˜n dS = 0 (13) S3 Now, by applying the ux condition (Equation (4)) on S3 , the nal weak form becomes ZZ ZZZ t=t2 @w dG (C)w d − (C) @t t=t1 G ZZ ZZZ @Xmax ˜ · ∇(C) ˜ ∇w dG − dS w(C) + @t G S1 x ZZ + wk C(C − C0 )f(C ∗ ) dS = 0 (14) S3 In order to solve equations (6)–(10) using the nite element method, the domain G must be discretized. However, the shape of the surface, 3 is unknown a priori because V is an unknown function of time. Therefore, the solution of equations (6) – (10) must be found piecewise in time. The growth velocity, V , will be taken to be constant for some time, t. During this time, each lamella will grow by an amount V‘ t, where V‘ is the growth velocity of lamella ‘. The shape of 3 is then known during t. A slice of G bounded by (t1 ) and (t2 ) (where t2 = t1 + t) is then discretized using a single layer of 8-node isoparametric space–time nite elements. An example element is shown in Figure 4. Nodes 1– 4 lie in the plane t = t2 and nodes 5–8 in the plane t = t1 . Also, nodes 1, 2, 5 and 6 lie in the plane y = y2 and nodes 3, 4, 7 and 8 in the plane y = y1 . The original mesh will deform with time as the lamellae grow. During one timeslice t n , each crystal will grow by a length of V‘ · tn , and the nodes in (t2 ) will be shifted with respect to the nodes in (t1 ). The nodes on 3 will be shifted by exactly V‘ t in the +x direction and Int. J. Numer. Meth. Engng., 40, 2679–2692 (1997) ? 1997 by John Wiley & Sons, Ltd. A SPACE–TIME FINITE ELEMENT MODEL 2685 the nodes on 1x by V‘max t, where V‘max is the maximum of all V‘ . The new x positions of the interior nodes are determined by linearly interpolating between the positions of the nodes on the boundary. The elements to be used are 8-node three-dimensional linear isoparametric elements. There is one degree of freedom per node; this degree of freedom is the nodal concentration, C. Each element in x; y; t space is mapped into ; ; space using the following transformation: x = y= t = 8 P m=1 8 P m=1 8 P xm Nm (15) ym Nm tm Nm m=1 where (xm ; ym ; tm ) are the nodal co-ordinates in G. Nm are the linear interpolation functions (1 − ) (1 − )(1 − ) (1 − ) (1 − ) (1 − )(1 − ) (1 − )(1 − )(1 − ) (1 − )(1 − ) (16) with , and ranging from 0 to 1 within the element. The concentration, C within an element is approximated by 8 ˆ = P Cm Nm C (17) m=1 where Cm are the nodal concentrations. The weak form in equation (14) must be discretized in order to solve for the unknown nodal concentrations. The derivatives of C are approximated by 8 ˆ P @(C) @Nm = Cm @z @z m=1 (18) where z is any one of the three independent variables x, y and t. Following a Galerkin formulation, the arbitrary function w is replaced by the interpolation functions Nn given above, and an approximate weak form is constructed: ZZZ ZZZ ZZ t=t2 ˆ n ) d − ˆ @Nn dG + ˆ dG ˜ n · ∇˜ C ∇N (CN C @t t=t1 G G ZZ ZZ ˆ max dS + ˆ − C0 )f(C ∗ ) dS = 0 Nn CV − kNn (C (19) ‘ ‘ S1x ? 1997 by John Wiley & Sons, Ltd. S3 Int. J. Numer. Meth. Engng., 40, 2679–2692 (1997) 2686 K. M. KIT AND J. M. SCHULTZ This can be put into the matrix form KC = f (20) where K is an 8 × 8 element matrix, C is the vector of nodal concentrations, f is an 8 × 1 vector, ZZ ZZ ZZZ @Nn Knm = dG (Nm Nn ) d − (Nm Nn ) d − Nm @t 1 2 G ZZ ZZZ @Nn @Nm @Nn @Nm + De2 dG − + Nn Nm V‘max dS @x @x @y @y G S1 x ZZ k Nn Nm V‘ dS (21) + S3 and fn = kf(C‘∗ )C0 ZZ S3 Nn dS (22) The derivatives of Nn with respect to x; y; t in terms of ; ; can be found by @N @Ni i @ @x @N @Ni i −1 =J @y @ @N i @Ni @t @ (23) where J is the Jacobian of the isoparametric transformation. The element matrix K is a function of the size and shape of the element. Elements of the same size and shape will have identical element matrices. The terms Knm are numerically integrated using Gauss –Legendre quadrature. The terms in the integrands in the rst three and last two integrals in equation (21) are polynomials of no higher than third order, and these terms can be exactly integrated using two-point Gauss– Legendre quadrature. However, some terms in the integrand of the fourth integral in equation (21) are of the form N=( + constant) where N is an integer between 0 and 4. These terms cannot be integrated exactly using Gauss–Legendre quadrature. However, ve-point quadrature was found to be accurate to within 2 parts in 1 × 107 . Assembly of all of the element matrices into the system matrix is done following standard practices. Before a solution can be obtained, an expression for the lamella growth rate, V , is needed. For a homopolymer, the growth rate is given by the Homan–Lauritzen theory:6 −Kg −Q exp Vn = A exp R(T − T1 ) T (T − Tm ) where Vn is the linear growth rate, A and K g are material constants, Q is is the crystallization temperature, Tm is the melting temperature, and T1 the glass transition temperature. In blends, however, expressions for the empirical. The form used here for A crystals growing from an A=B melt V = Vn (1 − C ∗ ) Int. J. Numer. Meth. Engng., 40, 2679–2692 (1997) or Vr = 1 − C0 − C ∗ 1 − C0 an activation energy, T is a temperature below growth rate are mostly is simply (24) ? 1997 by John Wiley & Sons, Ltd. A SPACE–TIME FINITE ELEMENT MODEL 2687 where CB‘ is the concentration of B in front of crystal ‘ at = 0. This expression assumes that the growth rate of A crystals is reduced only through dilution of A chains by B chains at the growth front. More complex formulations of the concentration dependence of V may be used in the present model; the simple relationship of equation (24) is, however, sucient to test the nite element model. For the rst timeslice, t1 , C for all nodes in 0 is given by the initial condition (equation (10)). The system equations are then solved for the nodal concentrations in 1 . The solution is then checked to see if the condition C = 0 is satised on 1 . If not, the values of Xmax and Ymax must be increased and a new solution found. If the condition is satised, then the computation can proceed to the next timeslice (bounded by 1 and 2 ). New values of C‘∗ and V‘ are calculated, and the positions of the nodes in 2 are determined. The system matrix is then recalculated. The nodal concentrations in 1 are prescribed to be the values obtained from the previous timeslice, and the nodal concentrations in 2 are solved for. Again, the condition C = 0 is checked on 1 and the simulation continues. NUMERICAL RESULTS In order to validate the correctness of this nite element scheme, a problem with a known analytical solution is solved. The problem is one of zone renement of a binary (A/B) alloy in which the component B has lower solubility in the solid than in the liquid. In practice, the alloy is forced to solidify at a constant velocity, V , in one direction only. As solidication occurs, B is partially rejected from the growing solid. The rate at which B is rejected is V (CL∗ − CS∗ ) or VCL∗ (1 − k). The concentration in the liquid, CL satises the one-dimensional time-dependent diusion equation D @CL @2 CL = @x2 @t (25) as x → ∞ (26) subject to the conditions CL = C0 D @CL + VCL (1 − k) = 0 @x CL = C0 at the solid=liquid interface at t = 0 (27) (28) Equations (25)–(28) are the one-dimensional analogue of the problem described in equations (1) – (5). Smith et al.7 have given a solution of the form CL (x0 ; t) where x0 is measured from the solid/liquid interface. The solution for k = 0 is obtained by taking the limit of CL (x0 ; t) as k → 0. 0 1 1 lim CL (x0 ; t) = C0 1 − erfc[R(x0 + Vt)] + e−Ux U (Vt − x0 ) erfc[R(x0 − Vt)] k→0 2 2 r (29) 2 V t −R 2 (x0 −Vt)2 e + erfc[R(x0 − Vt)] +2 D where 1 R= 2 ? 1997 by John Wiley & Sons, Ltd. r 1 Dt and U= V D Int. J. Numer. Meth. Engng., 40, 2679–2692 (1997) 2688 K. M. KIT AND J. M. SCHULTZ Figure 5. 2-D spatial domain which represents an innite growth plane Figure 6. Comparison of numerical and analytical solutions for growth of innite plane at constant velocity; k = 0, De = 0, C0 = 0·5, f(C ∗ ) = 1, t = 0·1 This problem has been solved using the present nite element scheme by setting k = 0 and V = constant and by using the domain in Figure 5. Recall that 1 is the approximated far eld boundary, 2 is the symmetry boundary, and 3 is the growth front. This domain simulates a domain which is innite in the y direction and eectively reduces the problem to one in x and t. A comparison of the numerical solution to the analytical solution is shown in Figure 6 where CL (x0 ; t) is shown at various times. The dotted lines follow equation (29) and the points are the results of the numerical simulation. The agreement is very good, and the maximum relative error is 0·18 per cent. In the case of an innite plane growing at a constant velocity there is no mechanism to keep C ∗ below a value of unity. Therefore, we must consider C to be in units other than volume fraction, such as atoms=cm3 . The numerical results of the problem in which V = f(C ∗ ) are shown in Figure 7 for growth of single lamellae into a two-dimensional innite melt. A timestep of t = 0 · 01, was used for the rst 10 steps and a t = 0· 1 was used for the remaining timesteps. When a t = 0· 1 was used for all timesteps, oscillations in C ∗ were seen at small times. However, at longer times C ∗ , converged to the same curves seen in Figure 7. The simulation where L = ∞ was done on the domain shown in Figure 5. For the one-dimensional problem (Lc = ∞ or De = 0), we expect C ∗ to approach a value of (1 − C0 ) at long times, due to the (1 − C0 − C ∗ ) dependence of Vr . As C ∗ continues to increase, Vr and the ux of rejected B chains will approach zero. This will result Int. J. Numer. Meth. Engng., 40, 2679–2692 (1997) ? 1997 by John Wiley & Sons, Ltd. A SPACE–TIME FINITE ELEMENT MODEL 2689 Figure 7. Results of 2-D simulations of the growth of isolated lamellae; k = 0, C0 = 0·5, S = h0i, t = 0·1 Figure 8. C ∗ versus reduced time for three adjacent lamella. Curves 1, 2 and 3 denote each of the three lamellae, 1 being on the line of symmetry and 3 being the outside member. k = 0, C0 = 0·5, De = 2·5, L = 1·75Lc , S = h6; 5; 4i, t = 0·1 in ever smaller increases in C ∗ . In this way an interface concentration of C ∗ = (1 − C0 ) will be approached but never reached. At long times this becomes a diusion controlled problem. Indeed, the value of d ln (d(C ∗ )=dt)=d ln t approaches a value of −1·5 at long times. However, we expect a single lamella growing into an innite melt to approach a steady state in which C ∗ ¡ (1 − C0 ). A steady state should develop for the same reasons that plate-like precipitates and lamellar eutectoid structures8 in metallic alloys grow under steady state conditions. In all of these situations, molecular or atomic diusion can occur in a direction perpendicular to the growth direction. These processes are diusion-coupled in the sense that diusional kinetics inuence which steady state is reached, but they are not diusion controlled. For a single lamella, the distance over which excluded species must diuse will scale with Lc . Recall that De can be thought of as the ratio of a characteristic time for crystallization, Lc =V , to a characteristic time for diusion in a direction perpendicular to the growth direction, L2c =D. From the data in Figure 7 it appears that the steady-state value of C ∗ which is reached decreases as De increases. An increase in De corresponds to faster diusional kinetics and /or slower crystallization kinetics. ? 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 2679–2692 (1997) 2690 K. M. KIT AND J. M. SCHULTZ Figure 9. CB (x; y) at t = 3; Other parameters as given in Figure 14 Figure 10. Contour plot of CB (x; y) at t = 3; Other parameters as given in Figure 8 The results of a simulation of the growth of ve adjacent lamellae are shown in Figures 14 – 10. Due to the symmetry condition, the solution space contains only 3 lamellae. Lamella 1 is centred on the x-axis and lamella 3 is farthest from the x-axis. The lamellae had initial reduced lengths of 6, 5 and 4. The mean interface concentration for each lamella is plotted in Figure 14. The interface concentration for lamella 1 is higher than for an isolated lamella under the same conditions (Figure 7, t = 2, De = 2·5). This increase is due to the presence of the adjacent lamellae. Figure 9 is a three-dimensional plot of CB (x; y) at t = 2, and Figure 10 is a contour plot of the same data. In Figure 9, the concentration in the crystals is plotted as C0 for clarity purposes only. These regions are not part of the solution domain, due to the assumption that the concentration within the crystals is zero. The concentration is highest in front of the growing crystals, and it decreases smoothly from those values. The whole solution domain (06x620; 06y625) is not plotted here. CONCLUSION A space–time nite element model has been presented which is capable of simulating the growth of many neighbouring needle-like crystals from a blend melt. It was assumed that one blend component cannot be incorporated into the crystals, and that the presence of this component at the growth front has an eect on the crystalline growth rate. The model determines the concentration eld in the melt surrounding the growing crystals by solving the time-dependent diusion equation. The use of space–time nite elements allows the position of each growth front to be tracked continuously in time. Simulation results were shown for the growth of an innite plane, a single needle-like crystal, and a group of needle-like crystals. It was seen that the growth of needle-like crystals approached a steady state and that the steady-state conditions were aected by the magnitude of the dimensionless Deborah number and the presence of other growing crystals. We intend to use this nite element Int. J. Numer. Meth. Engng., 40, 2679–2692 (1997) ? 1997 by John Wiley & Sons, Ltd. A SPACE–TIME FINITE ELEMENT MODEL 2691 formulation to study the eects of interface kinetics and solute diusion on crystallization kinetics and pattern formation during spherulitic crystallization of polymer blends. APPENDIX Derivation of weak form The integration of equation (11) by parts is not straightforward, due to the motion of the boundaries 1x and 3 . Equation (11) can be re-expressed as Z t2Z y2Z x2 (y;t) Z t2Z y2Z x2 (y;t) @(C) 2 dx dy dt − w w∇˜ (C) dx dy dt = 0 (30) @t t1 y1 t1 y1 x1 (y;t) x1 (y;t) where t1 , t2 , y1 = 0, y2 = Ymax , x1 (y; t), and x2 (y; t) dene the boundary of one timeslice of G. The rst term can be rewritten as Z t2Z y2Z x2 (y;t) @(C) dt (31) w dx dy @t t1 y1 x1 (y;t) and integrated by parts on t using Z y2Z U= y1 x2 (y; t) w dx dy and x1 (y;t) dV = @(C) dt @t Note that since the x limits of integration are functions of t, the Liebnitz rule must be used in the determination of dU . The result is ZZZ ZZ t=t2 @w dG w(C) d − (C) @t t=t 1 G Z t2Z y2 @x1 (y; t) (x1 (y; t); y; t) dy dt + w(C) @t t1 y1 Z t2Z y2 @x2 (y; t) (x2 (y; t); y; t) dy dt − w(C) (32) @t t1 y1 The x limits of integration, x1 (y; t) and x2 (y; t), describe the motion of the left and right boundaries of the solution domain, respectively. @x1 (y; t)=@t is non-zero only on S3 where it equals reduced velocity or f(C ∗ ). @x2 (y; t)=@t is the rate of change of the position of 1x which is equal to @Xmax =@t. The second term in equation (30) is now rewritten using the divergence theorem as ZZZ ZZ ˜ · ∇(C) ˜ ˜ ∇w dG (33) w∇(C) ·˜n dS + − S G The weak form can now be written as ZZ ZZZ ZZ t=t2 @w dG + w(C) d − (C) wCf(C) dS @t t=t1 G S3 ZZ ZZZ ZZ @Xmax ˜ · ∇(C) ˜ ˜ ∇w dG − dS + w(C) w∇(C) ·˜n dS − @t S1x G S1 ZZ ZZ ˜ ˜ w∇(C) ·˜n dS − w∇(C) ·˜n dS = 0 − S2 ? 1997 by John Wiley & Sons, Ltd. (34) S3 Int. J. Numer. Meth. Engng., 40, 2679–2692 (1997) 2692 K. M. KIT AND J. M. SCHULTZ ACKNOWLDEGMENTS This research was supported by the National Science Foundation, under Grant DMR-9115308. Additionally, we are grateful to Dr. Andre Benard for his continued interest and advice and particularly for his suggestions regarding nondimensionalization. REFERENCES 1. O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, Vol 2, 4th Edn, McGraw-Hill Book Company, London, 1991. 2. R. Bonnerot and P. Jamet, ‘A second order nite element method for the one-dimensional Stefan problem’, Int. j. numer. methods eng., 8, 811– 820 (1974). 3. R. Bonnerot and P. Jamet, ‘Numerical computation of the free boundary for the two-dimensional Stefan problem by space–time nite elements’, J. Comput. Phys., 25, 163–181 (1977). 4. A. Weill, A. Shitzer and P. Bar-Yoseph, ‘Finite element analysis of the temperature eld around two adjacent cryoprobes’, J. Biomech. Eng. Trans. ASME, 115, 374–379 (1993). 5. T. E. Tezduyar, M. Behr and J. Liou, ‘A new strategy for nite element computations involving moving boundaries and interfaces—the deforming-spatial-domain=space–time procedure: I. The concept and the preliminary numerical tests’, Comput. Methods Appl. Mech. Eng., 94, 339–351 (1992). 6. J. D. Homan and J. I. Lauritzen, ‘Crystallization of bulk polymers with chain folding: theory of growth of lamellar shperultites’, J. Res. Nat. Bur. Stand., 65A, 297–336 (1961). 7. V. G. Smith, W. A. Tiller and J. W. Rutter, ‘A mathematical analysis of solute redistribution during solidication’, Canad. J. Phys., 33, 723–745 (1955). 8. J. W. Christian, The Theory of Transformations in Metals and Alloys, 1st edn, Pergamon Press, Oxford, 1965. Int. J. Numer. Meth. Engng., 40, 2679–2692 (1997) ? 1997 by John Wiley & Sons, Ltd.