Patent Application
20190236225
The present invention relates to a method of deriving the buckling condition of a fiber moving in a fluid, a method of calculating the breakage condition, and a method of forecasting the fiber length distribution.
Fiber-reinforced plastics are being used in a variety of fields in recent years. Because the strength of a fiber-reinforced plastic is dependent on its microstructural fiber configuration, including fiber dispersion, orientation, and length, it is important to assess and control the microscopic states of fibers in plastics produced through molding processes. A wide variety of efforts have been made to forecast the dispersion and orientation of fibers in high-viscosity fluid through simulations in order to assess such microscopic states of fibers (see Nonpatent Documents 1 and 2). Meanwhile, although it is necessary to analyze fiber breakage processes in flows in order to forecast fiber length, quantitative evaluation of fiber breakage is difficult, and prediction of fiber length through simulation has not yet been attained.
Predictive analysis of orientation and breakage of fibrous objects moving in a flow has been attempted by Yamamoto et al. (see Nonpatent Documents 3 to 7). In the analysis of Yamamoto et al., a fiber is represented as a combination of spherical particles, and fluid resistance acting on individual spherical particles in a flow is summed up to calculate the force acting on a fibrous object.
Phelps et al. have derived the condition under which a cylindrical fiber buckles in a contraction flow by using an approximate solution with regard to fluid stress acting on a slender cylinder in a contraction flow (see Nonpatent Documents 8 and 9). The derived buckling condition model has been incorporated into molding software (see Nonpatent Document 10).
In the analysis of Yamamoto et al., fibers are expressed as a combination of spherical particles, and the fluid resistance acting on individual spherical particles is summed up to calculate the force of flow acting on the fibrous objects. However, the force acting on an object in a flow generally cannot be expressed by the sum of fluid resistance acting on individual spherical particles that are combined to mimic the object. For example, in the case of an object represented by a combination of spherical particles in a uniform flow, fluid resistance acting on each spherical particle is in proportion to the diameter of the particle in a Stokes region with a low Reynolds number; but meanwhile, fluid resistance acting on the combined particles overall is in inverse proportion to the particle diameter squared, because the number of spherical particles composing the object is in inverse proportion to the particle diameter cubed. In the limit case of particles of zero diameter, fluid resistance acting on the combined particles overall would be infinite. Although a method using such particle combinations can forecast the orientation of fibers in a flow, the problem is that it cannot forecast fiber breakage because it cannot determine the fluid forces acting on fibers.
The derivation process of Phelps et al. involves some theoretical questions, as described later.
An object of the present invention, which has been made in view of these problems, is to calculate the flow field around a cylinder in a contraction flow as addressed by Phelps et al. by a stricter method, and then solve the eigenvalue problem with regard to buckling that is caused by the resulting fluid stress distribution around the cylinder, in order to derive the buckling condition of a cylinder placed in a contraction flow. Another object is to provide a method of using the buckling condition to calculate the breakage condition, and a method of forecasting fiber length distribution.
A first aspect of the present invention to solve the problems described above provides a method of calculating a buckling condition, which is a condition under which buckling occurs, of a cylindrical fiber moving in a fluid when placed in a contraction flow, which is a flow toward a shrinkage in the longitudinal direction of the fiber. The buckling condition derivation method includes the step of multiplying a dimensionless fluid stress distribution, uniquely determined with regard to the aspect ratio r0′ of the cylindrical fiber, by μk (where μ represents fluid viscosity and k represents the contraction rate (velocity gradient in the length direction at the location where the fiber is present)) at a location where the cylinder is present, to obtain the fluid stress distribution acting on the cylinder surface in a contraction flow, and the step of using a minimum eigenvalue λ0, which is the threshold for buckling derived from the fluid stress distribution on the cylinder surface, to derive the buckling condition. The buckling condition is represented by the following equation (where E is Young's modulus and r0′ is the aspect ratio of the cylinder).
A second aspect of the present invention to solve the problems described above provides a method that uses the buckling condition described in the first aspect of the present invention to calculate a breakage condition, which is a condition under which breakage occurs, of a cylindrical fiber. The breakage condition derivation method includes the step of writing the value of μk when the buckling condition holds equality as μkbu, the step of writing the value of μk when the cylindrical fiber breaks as μkbr and writing the ratio of μkbr to μkbu, as a threshold rbr=μkbr/μkbu, the step of taking the threshold rbr as a fitting parameter and setting rbr to match experimental results or setting rbr based on structural analysis of fiber breakage incorporating a breakage model, and the step of determining the value of gk at breakage by replacing μk with μkbr/rbr and replacing the inequality sign with the equality sign in the buckling condition. The breakage condition is represented by the following equation.
A third aspect of the present invention to solve the problems described above provides a method that uses the breakage condition described in the second aspect of the present invention to forecast fiber length distribution, which is the distribution of cylindrical fiber lengths. The fiber length distribution forecasting method includes the step of determining a maximum value μkmax among historical data of μk up to a predetermined measuring point, and the step of calculating the cylindrical fiber's aspect ratio r0′max by assigning the maximum value μkmax to the breakage condition μkbr. Cylindrical fiber diameter df is then used to determine fiber length before breakage df/r0′max and fiber length after breakage dt/(2r0′max).
In the present invention, the flow field around a cylinder in a contraction flow as addressed by Phelps et al. is calculated by a stricter approach, and the eigenvalue problem with regard to buckling that is caused by the resulting fluid stress distribution around the cylinder is solved, providing a method of calculating the buckling condition for a cylinder placed in a contraction flow, and thereby providing a method of calculating the buckling condition, a method of calculating the breakage condition, and a method of forecasting fiber length distribution for a fiber moving in a fluid. Other advantageous effects of the present invention will be described in the embodiment below.
A preferred embodiment of the present invention will now be described in detail, with reference to the accompanying drawings. Throughout the specification and drawings, components having substantially the same functions are identified by the same symbols, without redundant description.
The following description refers to Nonpatent Documents 1 to 10 and, additionally, to Nonpatent Documents 11 to 20.
In general, the Stokes approximation holds for a flow around a fiber moving in a high-viscosity fluid, and viscosity is regarded as approximately constant near the fiber. Under such conditions, the governing equation of flow around the fiber is linear, as follows (see Nonpatent Documents 11 and 12).
∇p=μ∇2v, ∇·v=0 (I)
Here, v is the flow velocity vector, p is fluid pressure, and μ is viscosity. Flow is examined in a coordinate system that moves together with the fiber, and the time-derivative term is neglected in Equation (1). Here, a Taylor series is used to expand the flow field around the slender fibrous object, centered on the fiber center x0, as shown in Equation (2) below. Here, the subscripts i (i=1, 2, 3), j (j=1, 2, 3) represent the (x, y, z) components of the vector. Since the flow governing equation (1) is linear, flow fields can be superposed, and the flow fields can be examined separately for each term of the Taylor expansion.
The first term on the right-hand side in Equation (2) represents a uniform flow moving at the moving velocity v0 of the fiber, and can be neglected in view of the coordinate system that moves together with the fiber. Within the second term on the right-hand side, the term representing shear deformation among nine velocity gradient tensor components represents rotational flow around the center x0, and can be neglected in view of the rotational system that rotates together with the fiber. The remaining velocity gradient tensor components correspond to extension flow or contraction flow in the longitudinal direction of the slender fiber. Contraction flow is thought to be a factor causing fiber breakage, as others have already pointed out (see Nonpatent Documents 13 to 16). This embodiment examines fluid stress acting on a slender cylinder mimicking a fiber, placed in a contraction flow where the uniform flow and rotational flow components are excluded from the flow field around the fiber.
The contraction (extension) flow is generally represented by the following flow velocity distribution (see Nonpatent Document 17).
Here, (vx, vy, vz) are the (x, y, z) components of the flow velocity, and the flow is considered to contract from z=±∞ toward z=0. The velocity gradient in the z direction is k, which corresponds to the contraction rate (if k<0, the extension rate). The velocity components (vx, vy) perpendicular to the contraction direction (z) are dependent on a parameter m (0≤m≤1), where there is axisymmetric flow at m=0 and plane flow at m=1. Since buckling is caused by fluid stress acting on the cylinder surface in the z direction, and the effect of (vx, vy) is considered to be secondary, differences in (vx, vy) due to differing m values are neglected. Based on these considerations, flow around an slender cylinder having the z axis as its central axis placed in an axisymmetric contraction flow (Equation 4) at m=0 is discussed in Section (3), “Distribution of fluid stress acting on a cylinder in a contraction flow.” The formulation of Phelps et al. also assumes a similar axisymmetric contraction flow.
Flow around a cylinder having radius r0 and length 21 placed in an axisymmetric contraction flow represented by Equation (4) is considered, where the central axis of the cylinder is the z axis having a middle point at z=0. The governing equation of flow is Equation (1), and the boundary conditions at the cylinder's surface and infinity are expressed as follows.
Although the following development of the formula can be applied to either contraction flow (k>0) or extension flow (k<0), the following description is based on contraction flow (k>0), which is thought to be a direct factor in fiber breakage. Equation (1) and Boundary conditions (5) are made dimensionless, as follows.
Therefore, Equation (6) is as follows.
Since the only parameter included in Equation (6) is the aspect ratio r0′=r0/l of the cylinder, including the boundary conditions, the dimensionless flow field (flow velocity and pressure) obtained by solving Equation (6) is only dependent on r0′.
Equation (6) can be solved by computational fluid dynamics (CFD).
Since Equation (6) is only dependent on r0′, the resulting dimensionless flow velocity distribution (vr″, vz″) or (vr′, vz′) remains unchanged until r0′ varies. After flow velocity and pressure distribution around the cylinder are determined as shown in
For the viscous stress μδvz/δr acting on the side of a cylinder in an actual flow, the dimensionless velocity gradient of
Similarly, the fluid stress acting on the two ends of the cylinder can be determined by the product of the dimensionless pressure, the velocity gradient, and μk.
As described above, the fluid stress distribution acting on the surface of a cylinder in a contraction flow can be determined by the product of the dimensionless fluid stress distribution unique to the aspect ratio r0′ of a given cylinder and ik at the position of the cylinder. The dimensionless velocity gradient on the surface and the pressure at the two ends of the cylinder shown in
In the case of a cylinder that is sufficiently slender, i.e., having a sufficiently low aspect ratio (r0′=r0/l<<1), the effect of the two ends on flow around the cylinder disappears, except near the two ends. To determine axisymmetric flow around the cylinder in a case where the effects of the two ends can be neglected, Governing equation (1) is rewritten using the Stokes stream function ϕ(r,z) as follows:
In order to solve Equation (8), the following equation is written.
φ(r,z)=zf(r)
For f(r), Equation (9) is thereby obtained.
Equation (9) is integrated to obtain Equation (10) for f(r).
f(r)=c1r4+c2r2(log r′−½)+c3r2+c4 (10)
Here, c1, c2, c3, and c4 are integral constants. From Equation (8b), the velocity components (vr, vz) are derived as follows.
Here, r′=r/i. The integral constants c1, c2, c3, and c4 are determined for (vr, vz) by the boundary conditions (8c). However, the Stokes approximation is known to Include solutions that increase logarithmically at infinity (Stokes' paradox) (see Nonpatent Documents 11 and 12). Therefore, only c1 is written as 0 from the boundary condition at infinity, and c2 is determined from the boundary condition of the surface of the cylinder. The integral constants are determined as follows from the boundary conditions.
Flow velocity distribution is then determined as follows.
In the flow velocity represented by Equation (13), the first terms on the right side (kr/2, −kz) correspond to a contraction flow in a case where no cylinder is present, i.e., the primary flow; and the rest of the terms on the right side represent deviation from the primary flow. From Equation (13b), the velocity gradient on the surface of the cylinder is obtained in Equation (14).
In
This section will derive the buckling condition of a cylinder in a contraction flow, using fluid stress distribution acting on the surface of the cylinder in the contraction flow as determined in the preceding section. The procedures of buckling and breakage analysis will also be described.
Euler was the first to analyze the buckling condition in the case of inward forces acting on the two ends of a cylinder, and this analysis is known as Euler buckling. Dinh et al. (Nonpatent Document 16) approximated the fluid stress acting on the surface of a cylinder in a contraction flow, and Phelps et al. (see Nonpatent Document 18) used this approximation to derive the buckling condition for a cylinder in a contraction flow, with reference to the results of Euler buckling. Dinh et al. assumed that the velocity gradient at the surface of the cylinder is in proportion to the contraction rate at a position far from the surface of the cylinder (i.e., a place where the cylinder is not present), calling the proportionality coefficient in this case a “resistance coefficient,” but no theoretical derivation was carried out on the resistance coefficient, which remains an uncertain parameter. Thus the buckling condition Phelps et al. derived based on the approximate solution by Dinh et al. also incorporates this same resistance coefficient, which means that it is merely a fitting parameter. In Section (3-2), “Theoretical solution for flow around a cylinder excluding the two end regions,” flow velocity distribution is derived in the case of negligible effect of the two ends of the cylinder, and this actually corresponds to a theoretical derivation of the resistance coefficient in the paper by Dinh et al. The “resistance coefficient” in the paper by Dinh et al. is represented by the product of the velocity gradient of Equation (14) and the viscosity. Therefore, the approximate solution by Dinh et al. can be replaced with the solution obtained in Section (3-2), “Theoretical solution for flow around a cylinder excluding the two end regions,” to eliminate the uncertain parameter.
The derivation of the buckling condition by Phelps et al. has another problem, in addition to the above. In order to apply the Euler buckling problem as is, Phelps et al. replaced the compressive force caused by flow acting on the center of the cylinder, as determined from the approximate solution by Dinh et al., with the inward forces acting on the two ends of the cylinder in the Euler buckling problem to derive the buckling condition. However, no ground for such replacement is shown. Accordingly, doubt remains concerning the reasonableness of the results of Phelps et al.
In this section, the buckling condition of a cylinder in a contraction flow is derived from fluid stress distribution on the surface of the cylinder in the contraction flow, as obtained in the preceding section, instead of the approximate solution by Dinh et al. Instead of directly incorporating the results of a Euler buckling analysis, as done by Phelps et al., the eigenvalue problem derived using fluid stress distribution on the surface of the cylinder, as determined in the preceding section, is solved to derive the buckling condition.
In determining the buckling condition of a cylinder placed in a contraction flow, as shown in
N sin θ1−T cos θ1−(N+dN)sin θ2+(T+dT)cos θ2=0 (15)
Here, the definitions of variables N, dN, T, dT, F, θ1, θ2, θ, and ds are as shown in
Therefore, Equation (15) derives Equation (16).
Both sides of Equation (16) are divided by dz. Because of the balance of forces in the axial (z) direction, the expression dN≈−Fds≈−Fdz holds. Based on these facts, Equation (16) derives Equation (17).
The relationships shown in the following Equation (18) are then used, introducing bending moment M, Young's modulus E, and second moment of area I (=πr04/4).
Equation (17) then derives Equation (19).
Because of the balance of forces in the axial direction inside the cylinder, N in Equation (19) is expressed as follows, as a function of the sum P of fluid pressure and viscous stress acting on the end faces of the cylinder, and viscous stress F acting on the side of the cylinder.
The sum P and the velocity component vz in Equation (20) are calculated using the flow field determined in Section (3), “Fluid stress distribution acting on a cylinder in a contraction flow.”
In the case of a cylinder that is sufficiently slender (r0′<<1) for the velocity gradient on its surface to be approximated by Equation (14), the effects of pressure and viscous stress at the two ends of the cylinder is negligible (P=0), and thus Equation (20) can be expressed as Equation (21).
By substituting Equation (21) into Equation (19), we obtain Equation (22).
The independent variable is converted from z to z′ (=z/l). Equation (22) then derives Equation (23).
The boundary condition for Equation (23) is a condition where the two ends of the cylinder are free ends; in other words, bending moment is zero, and shear force in the x direction at the two ends of the cylinder is zero. For the boundary condition that bending moment is zero, Equation (24) is derived from Equation (18a).
Although shear force in the x direction is generally expressed as a linear combination of T and P, the value of P can be zero for a sufficiently slender cylinder; hence, the shear force equals T. As a result, the boundary condition of zero shear force in the x direction is represented by T=0, and Equation (18b) is used to derive Equation (25).
Differential Equation (23), along with boundary conditions (24) and (25), is an eigenvalue problem where A is an eigenvalue, and the eigenvalue and the eigenfunction δ(z′) are obtained by numerically solving Equations (23), (24), and (25). If λ0 is the minimum eigenvalue, the actual value of λ0 determined by numerical calculation is 5.09.
The buckling condition corresponding to Euler's buckling condition is derived by using Equation (23b) and the minimum eigenvalue of λ0, such that λ≥λ0, and is represented by Equation (26) or (27).
Unlike the buckling condition equation obtained by Phelps et al. using the approximate solution of Dinh et al., Equation (27) does not incorporate uncertain parameters such as a resistance coefficient.
The buckling condition of a cylinder in a contraction flow as determined in Section 3 can be applied to the results of fiber orientation analysis to forecast the position within a flow at which the cylindrical fiber will buckle.
There are two approaches in fiber orientation analysis methods: a method using ellipses or ellipsoids (see Nonpatent Documents 19 and 20), and a method using combinations of particles (see Nonpatent Documents 1 and 3 to 7). Either method can be used.
The method for forecasting the buckling of a fiber in a flow is presented in Section (4-2), “Forecast of fiber buckling using fiber orientation analysis.” However, because buckled fibers do not necessarily break, buckling is only a precondition for breakage. In order to forecast fiber length distribution, it is necessary to either forecast fiber breakage or correlate the breakage condition to the buckling condition. Considering the case of a sufficiently slender fiber (r0′<<1), if μkbu is the value of Ilk when equality holds for the buckling condition (Equation 27), then fiber buckling will begin at a time when μk exceeds μkbu at a place where a fiber is present in the fluid. Breakage is thought to occur when the fiber moves to a place having a higher μk value. The value of μk when breakage occurs is written as μkbr, and the ratio of μkbr to μkbu is written as rbr=μkbr/μkbu. Because it is proportional to the fluid stress μk acting on the surface of the fiber in the contraction flow, rbr (>1) indicates that breakage may occur when the fluid stress acting on the fiber increases to several times its value at the time when buckling begins. Determination of the breakage threshold μkbr or rbr will now be explained.
Concerning an approach for calculating the breakage condition or threshold, a first conceivable method of determining the breakage threshold would involve structural analysis incorporating a cylindrical fiber breakage model. However, to actually perform structural analysis of breakage, a fiber breakage model would be required.
Another conceivable method for determining the breakage condition or threshold involves approximation of rbr by means of a constant, in which rbr is regarded as a fitting parameter and determined by matching it to experimental results. This method requires experimental data on fiber length distribution, but actually, because the value of rbr is thought to be the same with regard to the same fibers, once a value of rbr has been determined by fitting to specific experimental data, that value can be used for the same fibers generally. After the threshold rbr has been obtained, μkbr can be determined from the buckling condition (Equation 27) by replacing μk with μkbr/rbr and replacing the inequality sign with the equality sign, as follows.
As an example of breakage forecasting analysis, the experimental measurement of fiber length distribution by Phelps et al. (see Nonpatent Documents 8 and 9) is traced to verify the extent to which fiber length distribution can be forecast by the present fiber breakage forecasting theory using the buckling condition. Phelps et al. fed a liquid containing cylindrical fibers having a diameter of 17 μm from the upper end of a central protrusion of a device shown in
In addition to measuring fiber length, Phelps et al. compared the fiber length distribution obtained using the fiber breakage forecasting theory proposed by Phelps et al. with their actual measurement data. However, the analysis by Phelps et al. ignores the device's central protrusion and addresses only the disk; and in addition, instead of fiber length distribution at the time of inflow into the device, they take measurement data of fiber length distribution at point A as the inflow condition of fiber length, and then use their analysis to forecast fiber length distribution at points B and C.
However, because flow velocity is slower toward the outside of the disk of the device shown in
In the present analysis, forecasting of fiber length distribution at point A is based on fiber length distribution at a point that is further upstream. Because the papers by Phelps et al. do not mention the dimensions of the device's central protrusion, the present analysis omits the device's central protrusion, as Phelps et al. also did in their analysis, and addresses only the inner portion of the disk. In the flow field analysis, axial symmetry is assumed, and the fluid is assumed to flow downward at a constant velocity from a region having a radius of 5 mm (an approximation because the actual dimensions are unknown) at the center of the upper face of the disk. The inflow volume, viscosity model, and other conditions and parameters stated in the papers by Phelps et al. are used without modification.
With regard to the feeding of fibers, approximately 5000 combined particles mimicking the fibers of
It is not possible to strictly reproduce the fiber length distribution measurement results of Phelps et al. due to uncertainties when tracing some points. Detailed measurement data on fiber length distribution is necessary for verification, but none can be found in the literature; so the measurement data of Phelps et al. Is used in the present analysis. Thus, an object of the present analysis is not to verify the quantitative accuracy of the fiber breakage forecasting theory, but to verify whether it would be generally feasible or completely impossible to reproduce the experimental results on fiber length distribution based on this fiber breakage model.
The specific method for calculation of the fiber length in
In the present analysis, because fiber length at the time of feeding is assumed to be relatively long at 6 mm, it is thought that most fibers are broken before passing through point A. However, because fibers have length distribution at the time of inflow into the device in an actual experiment, short fibers are also thought to exist. Thus, it is thought that some of the fibers do not break, but pass through point A while maintaining the same length as at the time of feeding.
A method of readily determining fluid stress distribution acting on the surface of slender cylindrical fibers in a contraction flow using CFD, etc. is proposed. A theoretical solution is derived for the case where the effects of the two ends of the cylinder are negligible, and based on the resulting fluid stress distribution, the eigenvalue problem is solved to derive a buckling condition. The buckling condition is applied to the results of fiber orientation analysis to propose a method of forecasting the position where fibers in a flow will buckle. The breakage condition is calculated using the buckling condition to propose a method of forecasting fiber breakage in the flow. The experimental results of Phelps et al. for measurement of fiber length distribution are traced, indicating that the present fiber breakage forecasting theory is effective for forecasting fiber length distribution.
As described above, according to the embodiment, the flow field around a cylinder in a contraction flow as addressed by Phelps et al. is calculated by a stricter method, and the eigenvalue problem is then solved with regard to buckling that is caused by the resulting fluid stress distribution around the cylinder, in order to derive the buckling condition of a cylinder placed in a contraction flow. A method of deriving the buckling condition of a fiber moving in a fluid, a method of calculating the breakage condition thereof, and a method of forecasting the fiber length distribution are thereby provided.
Preferred embodiments of the present invention are described with reference to the accompanying drawings. Needless to say, the present invention is not limited by these examples. It will be apparent to those skilled in the art that various variations and modifications can be conceived within the scope described in the claims, and naturally these are understood as falling under the technical scope of the present invention.
| Number | Date | Country | Kind |
|---|---|---|---|
| 2016-152913 | Aug 2016 | JP | national |
| Filing Document | Filing Date | Country | Kind |
|---|---|---|---|
| PCT/JP2017/028124 | 8/2/2017 | WO | 00 |
