close

Вход

Забыли?

вход по аккаунту

?

656

код для вставкиСкачать
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.
Документ
Категория
Без категории
Просмотров
2
Размер файла
480 Кб
Теги
656
1/--страниц
Пожаловаться на содержимое документа