Investigating mathematical methods for high-throughput prediction of the critical buckling load of non-uniform wool fibers

Abstract

The sensation of prickle from textile garments is directly related to the force that a fiber protruding from the fabric surface can exert on the skin without buckling – its critical buckling load (CBL). Finite element modeling (FEM) has previously been used in the literature to predict CBLs for a set of 25 fibers with different along-fiber morphology. With a view to high-throughput analysis of fibers, we investigated two analytical methods that were potentially faster and less computationally intensive than FEM and applied them to calculate CBLs for the same set of fibers. In addition, we tested a numerical integration and gradient search (NIGS) method that we developed by adapting a previously published, non-FEM, numerical approach.

The analytical methods that we tested were either inadequately formulated or prone to instability. Our NIGS method was more reliable that the analytical methods (but slower to compute), and its results appeared more accurate than the published FEM results, based on an inconsistency metric that we developed.

The published FEM results and the NIGS predictions agreed within 5% for 60% of the fibers, and within 10% for 72% of the fibers (with differences ranging from 0.4% to 19.1%) and generally showed qualitative agreement on the response of CBL to fiber shape, with some notable exceptions.

The response of CBL to dimensional variation was complex. This, and the inconsistency between methods, highlights the need for caution when analyzing complicated biological structures, such as wool, and the value of verifying the reliability of any predictions from any approach.

Keywords

textile prickle critical buckling load wool fiber finite element modeling numerical modeling Euler–Bernoulli beam theory flexural rigidity

Single wool fibers are highly heterogeneous both in composition and shape; effectively, they are a multi-component cellular-based biomaterial. An adequate understanding of how the industrially relevant properties of single fibers are created from underlying composition requires consideration of structure at a range of scales from nanometers to micrometers,¹ and this has been the motivation for significant past modeling efforts.² Liu and Bryson³ addressed the composite nature of wool fibers by modeling the effective stiffness of a fiber as a function of the Young’s moduli of its cuticle, orthocortex and meso/paracortex, weighted by the relative abundances of these components, but did not extend their approach to investigate emergent effects of along-fiber variation in composition. Conversely, the effect of shape on buckling was investigated by He⁴ and He and Wang,⁵ who modeled smooth sinusoidal variations in diameter along the length of a fiber, but assumed its composition to be constant. For irregular shapes, explicit solutions are not necessarily available, and thus finite element modeling (FEM) methods were used. Munro and Carnaby⁶ also used FEM to connect ratios of idealized macrofibrils (the structural basis of orthocortex and paracortex) to single fiber curvature. Our long-term interest is also in describing how along-fiber variation in morphology, cellular and sub-cellular composition scales up to affect whole-fiber mechanical properties from biological data. We have chosen “prickle” as a property for initial investigation as it is a relatively simple, well understood and physically measurable single-fiber textile property.

The sensation of prickle from textile garments or surfaces is largely governed by the buckling behavior of individual fibers whose ends protrude onto the skin. A fiber will resist a degree of end-on force from the skin without buckling. The greater this exerted force, the greater the sensation of prickle, once a nerve sensitivity threshold has been reached.⁷ A critical buckling load (CBL) of greater than 75 mg has been shown to stimulate pain nerve receptors in the skin, and sufficient stimulation within a given area of skin results in a prickle sensation. Thus the CBL of a fiber (the minimum force that will cause buckling) is a major factor in the tactile acceptability of textiles, especially garments worn next to the skin (Figure 1).

Figure 1.

Fibers protruding from the surface of a textile garment are subject to an end-on force from the skin of the wearer, which generates a resultant force of the same magnitude exerted from the fiber onto the skin. If fiber can resist this force without buckling, the wearer of the garment will experience a sensation of prickle if the force is over a nerve sensitivity threshold. Discomfort increases with the magnitude of the force. If the fiber buckles it causes very little prickle.

The CBL of a thin column or rod can be described by the Euler–Bernoulli beam theory, for example, see Timoshenko and Young.⁸ This theory was applied in mathematical studies of wool fiber stiffness by Veitch and Naylor,⁹ who modeled the wool fiber as a materially homogenous, circular rod with uniform cross-section, for which an explicit, analytical solution is readily obtainable. More recently, Asad et al.¹⁰ applied the theory to study a range of animal and plant fibers, using similar assumptions of homogeneity and a uniform, circular cross-section of fibers.

A natural extension of the previous approaches is to model how both variation in fiber morphology and composition along the length of the fiber might generate single-fiber stiffness, affecting CBL and, thus, fiber prickle. This is necessary for future studies that investigate how actual variations in real wool fibers affect CBL. FEM is widely used for studying the mechanical properties of irregular structures. For biological structures such as bone, the performance of FEM may depend on the number and shapes of mesh elements used to represent the structure.¹¹ Thus, accuracy may require manual input into mesh design rather than relying on a purely automated approach to mesh generation. This consideration is equally applicable in the context of complex wool fibers. Therefore, from the future aspiration of automated, high-throughput analysis of fibers from real fiber and sub-fiber structural data, we investigated alternative methodologies for finding the CBL of irregular fibers.

Methods

Equations for critical buckling load

The buckling of fibers is analogous to the buckling of columns in engineering applications, despite the large difference in scale. The governing equation for buckling of columns not subject to axial restraints can be obtained from the moment-equilibrium equation for a free body and was given by Coşkun and Atay¹² as a fourth-order differential equation (DE)

\frac{d^{2}}{{dx}^{2}} (E (x) I (x) \frac{d^{2} w}{{dx}^{2}}) + P \frac{d^{2} w}{{dx}^{2}} = 0

(1)

where

0 \leq x \leq L

is the distance from one end of the fiber (which is of length L),

w \equiv w (x)

is the displacement of the fiber at each point

x,

E and I (respectively the longitudinal stiffness and moment of inertia) are, in principle, functions of x, and P is the buckling load. Values of P that satisfy equation (1) will form a discrete set of numbers corresponding to different buckling modes. The CBL is the smallest, strictly positive solution to P, and corresponds to the first mode of buckling.

We can express equation (1) in a non-dimensional form as

\frac{d^{2}}{{dx}^{2}} (ζ (x) \frac{d^{2} w}{{dx}^{2}}) + β \frac{d^{2} w}{{dx}^{2}} = 0

(2)

where x has been transformed by

x \to x / L

, so that

x \in [0, 1],

and ζ is the normalized flexural rigidity

ζ (x) \equiv E (x) I (x) / ζ_{0} ζ_{0} = E_{0} I_{0}

(3)

where

E_{0}

and I₀ are reference values of E and I evaluated at

x = 0

, and β is the normalized buckling load

β \equiv {PL}^{2} / ζ_{0}

(4)

Wool fibers vary in diameter and cross-sectional shape along the short distances relevant to our prickle scenario (Figure 2). In our model we assume that the fiber has an elliptical cross-section along its length, and the moment of inertia of the fiber is described as

I (x) = \frac{A (x)^{2}}{4 π ε (x)}

(5)

where

A (x)

is the cross-sectional area of the fiber at distance x along its length, and we are using the definition for ellipticity of

ε (x) = r_{1} / r_{2}

, where r₁ and r₂ are the lengths of the major and minor axes, respectively.

Figure 2.

A scanning electron microscope image of wool fibers from a Perendale sheep showing the typical irregularities that occur in what would otherwise be a fiber of elliptical cross-section (image by the authors using standard scanning electron microscope methods).

Figure 3 depicts the assumptions made for the analysis of prickle. The end of the fiber that is anchored in the garment can be considered to be clamped (at $x = 0)$ . The other end can be considered pinned or hinged (at $x = 1$ ); that is, fixed into the skin and able to rotate but not slide. The projection of the buckled fiber onto the x-axis will be less than unity, but this difference is assumed negligible in the formulation. These assumptions give rise to two classical boundary conditions¹³

Clampedend : w = 0 and \frac{dw}{dx} = 0

(6)

Pinnedend : w = 0 and \frac{d^{2} w}{{dx}^{2}} = 0

(7)

Figure 3.

Representation of the buckling problem for a textile fabric against skin, where the fiber length has been normalized to the unit length, and the fiber is subject to a normalized force, β, by the skin, with buckling occurring if this force is greater than or equal to the normalized critical buckling load, $β_{cr}$ , of the fiber. The fiber is assumed to be firmly anchored in the textile fabric, with no rotation or sideways movement of the fiber at the fabric surface ( $x = 0$ ). The end of the fiber against the skin (at x = 1) is assumed to be pinned or hinged, so that there is no sideways motion, but rotation is possible. This gives rise to the boundary conditions described in equations (6) and (7). The projection of the buckled fiber onto the x-axis will be less than unity, but this difference is assumed to be negligible in the formulation.

Integrating equation (2) twice and using the boundary conditions gives

ζ (x) \frac{d^{2} w}{{dx}^{2}} + β w - a (1 - x) = 0, a = ζ (0) (\frac{d^{2} w}{{dx}^{2}})_{| x = 0}

(8)

However, many approaches to solving the buckling problem use the form in equation (2).

Solution approaches alternative to FEM

Exact, closed-form solutions for w and β are not generally available, although solutions to equation (2) have been calculated for a few special cases.^14–19 A number of analytical approaches have been used to obtain approximate solutions, including a functional perturbation method involving convolutions of ζ with Green’s functions,²⁰ application of the variational iteration method (VIM),^12,21,22 homotopy analysis²³ or by assuming a power series solution to w and generating linear relations to identify coefficients.^13,15

Numerical approaches other than FEM have been used to investigate column buckling, such as the boundary element²⁴ and differential quadrature²⁵ methods. Alternative approaches involve the numerical solution of equation (2) (or equivalent first-order formulations involving tangent angles and displacements) with an iterative search to find unknown initial values and buckling load.^26–28 Rychlewska²⁹ used a piecewise, exponential approximation of ζ, generating an eigenvalue problem incorporating continuity conditions across segments.

Analytical approaches have generally been developed using relatively simple functional forms of ζ, which are more applicable to structures of engineering interest than complex biological structures for which these approaches may become intractable. However, they are useful for testing the performance of other methods that yield approximate solutions and can be relatively fast to compute compared to numerical approaches that do not suffer the same tractability issues, but often entail a balance between precision of representation and calculation versus computation time and/or memory.

Investigating alternative approaches to FEM

We evaluated three non-FEM methods by analyzing the same set of fiber morphologies from the FEM studies of He⁴ and He and Wang⁵ and then comparing our results with those in their publication. For each of the fibers considered by He⁴ and He and Wang,⁵ elasticity was specified as $E =$ 3000 MPa (constant along the fiber) and fiber length specified as $L = 1.0$ mm. The diameter, d, of each fiber was varied sinusoidally along the length of the fiber as $d = d_{0} (1 - r \sin (2 m π x + θ))$ , where $d_{0} = 30$ µm for all fibers, and fibers varied in the combination of values for r (which took values of 0, 0.1, 0.3 or 0.5), m (which took values of 1, 5 or 10 mm⁻¹) and the offset parameter θ (which took values of 0, $π / 2$ , π, or $3 π / 2$ ), and $0 \leq x \leq 1.0 mm$ . Here, r controlled the amplitude of the diameter variation along the fiber, m controlled the number of sinusoidal modulations along the fiber and θ governed the positions along the fiber where the diameter was the thickest and thinnest. CBLs were reported as masses (mg) rather than forces, in these studies, for this conversion we used the value $g = 9.81$ ms⁻². The normalized flexural rigidity is obtained from equations (3) and (5)

ζ = (1 + r \sin (2 m π x + θ))^{4}, ζ_{0} = \frac{π d_{0}^{4}}{64} E

(9)

In equation (9) we have made the transformation, $m \to m \times 1.0 mm$ and $x \to x \times 1.0 {mm}^{- 1}$ .

To obtain the actual CBL (P) from the normalized CBL (β) in equation (4) and express it in terms of mg to facilitate comparisons between studies, multiplication by a scale factor of 12.2 mg is necessary

P = \frac{π d_{0}^{4}}{64} \frac{E}{L^{2} g} β = 12.2 β

(10)

The challenge for each method was to substitute equation (9) into equation (2) and solve for the smallest positive value of β. We investigated this for the same combinations of r, m and θ considered in the FEM studies.^4,5

The first of the three methods we investigated was the VIM developed by He³⁰ to solve the general DE, $Lu (t) + Nu (t) = g (t,)$ where L and N are linear and non-linear operators and g is an inhomogeneous term. The VIM provides the iterative solution $u_{n + 1} (t) = u_{n} (t) + \int_{t_{o}}^{t} λ {{Lu}_{n} (s) + N {\tilde{u}}_{n} (s) - g (s)} ds$ , where λ is a general Lagrange multiplier identified via variational theory, u_n is the nth approximation to the solution and ${\tilde{u}}_{n}$ denotes the restricted variation $δ {\tilde{u}}_{n} = 0$ . Sufficient iterations would yield an approximation that converges to the true value. The VIM has proved useful in the analysis of certain buckling problems where flexural rigidity has varied along the length of a column.^12,22 Coşkun and Atay¹² presented the following VIM correction functional for general Euler buckling represented by equation (2)

w_{n + 1} (x) = w_{n} (x) + \int_{0}^{x} λ (ξ, x) (\frac{d^{4} w (ξ)}{d ξ^{4}} + \frac{2}{ζ (ξ)} \frac{d ζ (ξ)}{d ξ} \frac{d^{3} w (ξ)}{d ξ^{3}} + \frac{1}{ζ (ξ)} [\frac{d^{2} ζ (ξ)}{d ξ^{2}} + β_{0}] \frac{d^{2} w (ξ)}{d ξ^{2}}) d ξ

(11)

λ (ξ, x) = \frac{1}{6} (ξ^{3} - 3 x ξ^{2} + 3 x^{2} ξ - x^{3})

(12)

with the initial approximation

w_{0} = {Ax}^{3} + {Bx}^{2} + Cx + D

. Imposing the boundary conditions on the solution generates an eigenvalue problem whose characteristic equation is a polynomial in

β_{0}

, with the smallest positive, real root being the CBL.¹² The challenge is solving the integral, which can quickly become intractable for complicated functions of ζ, such as equation (9). Coşkun and Atay¹² overcame this problem for exponentially varying functions of ζ by representing

1 / ζ

using a Taylor series expansion. For our functional form of ζ (equation (2)), the Taylor series for

ζ^{″} / ζ

ζ' / ζ

and

1 / ζ

all have very small radii of convergence unless an excessive number of terms is used. Thus in addition to the Taylor series, we also represented these functions using ordinary polynomial functions fitted to the curves numerically for each combination of parameter values. We used the symbolic package Maxima (http://maxima.sourceforge.net) to perform the Taylor expansion and numerical polynomial approximations, as well as the calculations in equation (11) and solution of the resulting eigenvalue problem.

The second method investigated (which we will refer to as the power series method or PSM) was presented by Huang and Luo,¹³ who assumed the solution to equation (2) to be an Nth-order power series expansion

w (x) = \sum_{i = 0}^{N} \frac{a_{i} x^{i}}{i!}, 0 \leq x \leq 1

(13)

and thus turned the problem into one of finding the

N + 1

coefficients,

a_{i}, i = 0, \dots, N .

From equation (2) one can obtain a number of relationships between the coefficients

where

C_{i}^{m} \equiv \frac{i!}{(i - m)!}

(15)

Substituting equation (13) into the boundary condition equations gives four more relationships, which together with equation (14) can be written as

T_{1} a + β T_{2} a = 0

(16)

$T_{1}$ and $T_{2}$ are $(N + 1) \times (N + 1)$ matrices, and $a = (a_{0}, a_{1}, \dots, a_{N})^{T}$ is a vector of the coefficients in equation (13). Solving equation (16) gives rise to a generalized eigen problem in which the smallest positive eigenvalue corresponds to the CBL, and the corresponding eigenvector gives an arbitrary multiple of the displacement of the buckled column/fiber.

Equation (9) can be exactly represented as a finite sum of sine and cosine terms, which we used in the PSM. With a view to analyzing real fibers whose flexural rigidity is available from data measured at discrete points, we also investigated the PSM using both Fast Fourier Transform (FFT) representations (PSM FFT) and polynomial representations (PSM Poly) of ζ to compare with the PSM using the exact representation (PSM Exact). The FFT and polynomial representations were obtained by calculating values for ζ at discrete points of x, for different values of r, m and θ and using, respectively, the fft and polyfit commands in MATLAB. All three PSM methods were implemented in MATLAB (The Mathworks Inc., Natick, MA) and used the eig function to solve the generalized eigen problem. Following the recommendation of Huang and Luo,¹³ we used a polynomial of order 14 to represent the buckling deformation, and thus the dimensions of the matrices involved were 15 × 15.

The third method tested was motivated by the approach of Lee et al.,²⁷ who numerically integrated equation (2) from $x = 0$ to 1, which would require values for β and the initial values of w and its first, second and third derivatives. Lee et al.²⁷ supplied two initial values from the boundary conditions at $x = 0$ , and treated β and apparently the remaining two initial values as parameters that could be varied using the determinant search method and the Regular Falsi method to achieve the boundary conditions at $x = 1$ , thereby identifying values for the CBL.²⁶

We adapted this version by making use of the boundary conditions at both ends of the fiber to determine the values of w and its first, second and third derivatives at $ξ = 0$ , which for a clamped/pinned fiber are

w (0) = 0, w^{'} (0) = 0, w^{″} (0) = - ζ (0), w^{″'} (0) = ζ (0) + ζ^{'} (0)

(17)

Equation (2) is readily expressed as a system of first-order ordinary differential equations (ODEs) with initial values given by equation (17), leaving just β as the unknown. We used the MATLAB function fminsearch (which implements a Levenberg Marquard gradient search) to find the CBL, β, using the MATLAB integration routine ode113 to perform the numerical integration. In each case we used a starting value of 0.1 for β (corresponding to 1.22 mg for the non-normalized problem). To avoid convergence on $β = 0$ , we used the following objective function

f = \frac{1}{β} (\frac{w (1)^{2}}{max (w^{2})} + \frac{w^{″} (1)^{2}}{max (w^{″ 2})}), β \neq 0

(18)

We refer to this numerical integration and gradient search method as “NIGS”. We note that because equation (2) is linear, equation (17) (and indeed any solution to w) can be scaled by a constant and still be valid for the same value of β.

Checking prediction consistency for the investigated approaches

Experimental testing of predictions would generally be required to validate theoretical predictions. If independent methods agree in their prediction of the CBL, β, then this may on its own suggest that the predictions are reliable. If methods differ, then determining which (if any) result is correct can be challenging without recourse to experiments. However, it may be possible to detect anomalies by checking the consistency of the predicted β with the corresponding solution for the deformation, w.

We calculated a dimensionless inconsistency score, γ, to use as an indicator of how well the solutions satisfied both equation (8), and the boundary conditions, defining it as

γ = \sqrt{\frac{\int_{0}^{1} (ζ (x) w^{″} + β w - a (1 - x))^{2} dx}{β^{2} \int_{0}^{1} w^{2} dx}} + \sqrt{\frac{w (1)^{2} + w^{″} (1)^{2}}{w^{″} (0)^{2} + w^{″'} (0)^{2}}}, a = ζ (0) w ″_{| x = 0}

(19)

The first term of equation (19) was obtained by taking the root mean square (RMS) of the left-hand side of equation (8) and dividing by the RMS of $β w$ to obtain a quantity invariant under scaling of w or β. The second term checks for deviation from the boundary conditions at the pinned end ( $x = 1$ ), and is normalized by the summed square of the clamped-end boundary conditions to make it invariant under scaling of w.

Ideally, predictions for β and w should yield a value of zero for γ when substituted into equation (19). However, we do not suggest that γ quantifies the uncertainty in the corresponding estimate of β. Further, γ does not detect whether or not β is the CBL, since γ should be zero for all buckling modes.

We calculated γ for from the predicted β and w for each method that we implemented. When a method expressed the solution to w as a polynomial, equation (18) was calculated by global adaptive quadrature using the integral function in MATLAB. When the solution to w was a set of numerical values corresponding to discrete values of x, ζ was correspondingly discretized and equation (19) was calculated by the trapezoidal method, using the trapz function in MATLAB. Since the values of $β_{cr}$ provided by He and Wang⁵ from their FEM studies were not accompanied by expressions for the corresponding deformation, we estimated w by integrating equation (2) using the FEM-predicted $β_{cr}$ , starting with initial conditions specified in equation (1).

Results

The uniform case was the only case in which the VIM performed well. In the cases where

r \neq 0

, the Taylor’s expansions of

1 / ζ

required an excessive number of terms, which ruled out its use for the VIM. For example, for the case

r = 0.1

m = 1

and

θ = 0

, the Taylor’s expansion around

x = 0.5

required 180 terms for a good approximation (converging within 2% of ζ), and this quickly exhausted the memory limits of a typical desktop computer (the degree of the polynomial solution would increase by 182 with each iteration). Therefore, the VIM was tested using a numerically fitted polynomial of degree 15 for the case when

m = 1

. However, the VIM did not converge to a solution, but fluctuated widely before the size of the symbolic problem overwhelmed the memory limits of the computer (over 20 iterations). Thus, the VIM was discarded as an alternative. Results for all other methods tested were recorded to one decimal place for consistency with He and Wang⁵ and are presented along with their FEM values in Table 1.

Table 1.

Predicted critical buckling loads ( $β_{cr}$ in mg) and discrepancy scores of predictions (γ in mg) for radially symmetric wool fibers with a length of 1.0 mm, constant elasticity of 3000 MPa and with sinusoidal variation of diameters along their length. The diameter along each fiber is described by $30 (1 - r \sin (2 m π x + θ))$ µm, where r governs the degree of variation (Level), m in (mm⁻¹) is the number of oscillations along the fiber length (Freq.), θ (in radians) governs the position of the peaks in diameter and $x \in [0, 1]$ is the normalized length of the fiber. Finite element modeling (FEM0 predictions presented by He⁴ and He and Wang⁵ are compared with predictions of four of the methods considered in this paper: the numerical integration and gradient search (NIGS) and the polynomial solution method in which the dimensional variation of the fiber was represented using either an exact equation (PSM Exact), a Fast Fourier Transform (FFT) approximation (PSM FFT) or polynomial approximation (PSM Poly). The discrepancy score expresses the mismatch between predicted $β_{cr}$ and predicted fiber deformation (a lower score indicating greater consistency). In order to calculate discrepancy scores for the FEM predictions, we estimated the deformation corresponding to each FEM value for $β_{cr}$ , by numerically integrating equation (2) with initial values given by equation (17)

	Variations in fiber			Predicted critical buckling load, $β_{cr}$ (mg), and discrepancy score, γ (mg)
	Level	Freq.	Phase	FEM		NIGS		PSM Exact		PSM FFT		PSM Poly
Case	(%)	(mm⁻¹)	(rad)	$β_{cr}$	γ	$β_{cr}$	γ	$β_{cr}$	γ	$β_{cr}$	γ	$β_{cr}$	γ
1	0	1	0	240.0	0.04	245.6	2 × 10⁻¹⁰	245.6	2 × 10⁻⁹	245.6	2 × 10⁻⁹	245.6	2 × 10⁻⁹
2	0.1	1	0	235.5	0.17	265.4	2 × 10⁻¹⁰	265.4	1 × 10⁻³	265.4	1 × 10⁻³	265.4	2 × 10⁻³
	0.3			217.2	0.02	225.9	7 × 10⁻¹¹	227.2	0.01	226.3	0.11	233.3	0.11
	0.5			118.3	0.02	123.0	1 × 10⁻¹⁰	135.5	0.21	134.5	0.21	167.2	0.21
3	0.1	1	0	235.5	0.17	265.4	2 × 10⁻¹⁰	265.4	1 × 10⁻³	265.4	1 × 10⁻³	265.4	2 × 10⁻³
		5		218.2	0.01	227.1	4 × 10⁻¹¹	238.9	0.21	238.9	0.21	234.7	0.91
		10		221.6	0.00	230.6	4 × 10⁻¹¹	248.7	0.21	248.7	0.21	867.7	0.81
4	0.3	1	0	217.2	0.02	225.9	7 × 10⁻¹¹	227.2	0.01	226.3	0.11	233.3	0.11
		5		138.4	2 × 10^–3	143.7	4 × 10⁻¹¹	230.3	0.51	229.6	0.51	169.7	1.6
		10		144.9	1 × 10^–3	150.1	9 × 10⁻¹¹	290.3	0.75	290.3	0.76	202.0	2.8
5	0.5	1	0	118.3	0.02	123.0	1 × 10⁻¹⁰	135.5	0.21	134.5	0.21	167.2	0.21
		5		56.6	1 × 10^–3	58.3	8 × 10⁻¹¹	254.9	0.91	255.3	0.91	134.6	5.9
		10		61.8	8 × 10^–5	62.1	1 × 10⁻¹⁰	370.4	1.4	370.4	1.4	155.4	8.6
6	0.1	1	0	235.5	0.17	265.4	2 × 10⁻¹⁰	265.4	1 × 10⁻³	265.4	1 × 10⁻³	265.4	2 × 10⁻³
			π/2	223.3	0.19	250.6	2 × 10⁻¹⁰	250.6	2 × 10⁻⁴	250.6	2 × 10⁻⁴	250.6	4 × 10⁻⁴
			π	201.5	0.02	209.9	8 × 10⁻¹¹	209.9	1 × 10⁻³	209.9	1 × 10⁻³	209.9	7 × 10⁻³
			3π/2	240.4	0.05	232.5	2 × 10⁻¹⁰	232.5	2 × 10⁻³	232.5	2 × 10⁻³	232.5	4 × 10⁻³
7	0.3	1	0	217.2	0.02	225.9	7 × 10⁻¹¹	227.2	0.01	226.3	0.11	233.3	0.11
			π/2	179.4	0.43	221.6	4 × 10⁻¹⁰	221.6	9 × 10⁻³	221.6	9 × 10⁻³	221.2	0.01
			π	127.7	0.01	133.0	3 × 10⁻¹¹	132.9	0.11	133.2	0.01	134.2	0.11
			3π/2	212.0	0.20	186.8	3 × 10⁻¹⁰	195.0	0.21	192.0	0.11	192.9	0.11
8	0.5	1	0	118.3	0.02	123.0	1 × 10⁻¹⁰	135.5	0.21	134.5	0.21	167.2	0.21
			π/2	108.7	0.62	132.9	6 × 10⁻¹⁰	132.2	0.21	132.2	0.21	129.7	0.49
			π	75.1	3 × 10⁻³	78.2	4 × 10⁻¹¹	82.3	0.31	82.3	0.31	85.3	0.33
			3π/2	126.9	0.25	113.3	7 × 10⁻¹⁰	138.3	0.28	138.3	0.28	143.3	0.23
A	0.1	1/2	π/2	223.3	0.07	243.1	2 × 10⁻¹¹	243.1	1 × 10⁻⁴	242.7	0.51	243.1	1 × 10⁻⁵
	0.3			167.0	0.04	191.5	3 × 10⁻¹¹	191.5	1 × 10⁻³	201.3	11	191.5	2 × 10⁻³
	0.5			94.4	0.01	103.8	2 × 10⁻¹¹	103.6	0.11	317.4	2.4	103.7	0.01
B	0.1	1/2	3π/2	223.6	0.11	232.8	5 × 10⁻¹⁰	232.8	1 × 10⁻⁴	232.9	0.21	232.8	4 × 10⁻⁶
	0.3			184.3	0.50	174.2	1 × 10⁻⁹	174.2	6 × 10⁻³	174.2	0.41	174.2	7 × 10⁻⁴
	0.5			100.1	0.55	98.5	9 × 10⁻⁹	98.6	0.01	99.2	1.1	98.5	0.01

For PSM FFT, 51 discrete values of ζ were used to obtain its Fourier series approximation. For PSM Poly a polynomial of degree 15 was used to represent ζ. A polynomial of this magnitude was not sufficient for accurately approximating ζ when m was 5 or 10, and using a higher order polynomial introduced spurious oscillations or approximated poorly at the ends of the fiber.

The results for PSM FFT and PSM Exact (which used the exact functional form of ζ) were in good agreement (within 0.1 mg) for only the instances satisfying { $r \leq 0.1}$ , { $m = 10$ } or { $r = 0.5$ , $m = 1$ , $θ \neq 0$ }, and for the combinations { $r = 0.3$ , $m = 1$ , $θ = π / 2$ } and { $r = 0.3$ , $m = 1 / 2$ , $θ = 3 π / 2$ }. With the exception of the latter combination, PSM FFT and PSM Exact generally differed for cases with $m = 1 / 2$ . PSM Poly and PSM Exact agreed well (within 0.1 mg) for only instances satisfying { $r \leq 0.1$ , $m = 1$ } or { $m = 1 / 2}$ . PSM Poly and PSM FFT agreed only for instances where they also agreed with PSM Exact. Based on the agreement with PSM Exact, PSM FFT performed better than PSM Poly on fibers with an integer number of sinusoidal modulations, while PSM Poly performed better than PSM FFT for fibers with only half a sinusoidal modulation.

PSM Exact, PSM FFT and PSM Poly showed very large variations in γ, ranging from a lower value of 2 $\times 10^{- 9}$ for all methods (for the uniform case) to upper values of 1.3, 11 and 8.6, respectively. This suggested that some predictions from PSM Exact, PSM FFT and PSM Poly should be approached with caution. Indeed, further investigation revealed an instability in these methods. This was detected by varying the degree of the polynomial approximating the deformation, which is given in equation (13). Small increases in the order of equation (13) could give very large fluctuations in the solution in some instances, with no sign of convergence. This is because the matrix T₁ in equation (16) could be near singular in these cases, leading to an unstable generalized eigen problem.

In order to ensure the NIGS was converging to an actual root and achieved the boundary conditions for $x = 1$ , we aimed to achieve the condition $β_{cr} f_{end} \leq 1 \times 10^{- 9}$ for each case, where $f_{end}$ is the converged value of equation (18). Although ode113 is an automatic multistep method, we specified initial and maximum step sizes of $InitialStep = 10^{- 20}$ and $MaxStep = 10^{- 5}$ to achieve the criterion. With these specifications, the NIGS demonstrated good consistency between predicted $β_{cr}$ and corresponding w, with $γ \leq 9 \times 10^{- 9}$ for all cases. This indicated that all NIGS predictions corresponded to actual buckling modes. Given that searches were started from a near-zero value of β in each case, we would suggest that the NIGS converged to $β_{cr}$ in each case.

The FEM results always differed by more than 0.1 mg from the results for any of the other methods for any of the model fibers. This was notable for the uniform case where all the tested methods gave the theoretically derived value while the FEM results differed by 5.6 mg. Agreement to either within 1 mg or to within 1% occurred only in one instance – for the NIGS result for { $r = 0.5$ , $m = 10$ , $θ = 0$ }. FEM and the NIGS agreed within 5% for 15 of the 25 unique combinations and within 10% for 18 combinations. The disagreement between FEM and the NIGS predictions varied from 0.3 to 42 mg (0.4–19%). Surprisingly, some of the larger disagreements occurred for fibers whose diameters had relatively little modulation, while the best agreement occurred for the fiber with the greatest degree of modulation.

However, on a qualitative level, the FEM and NIGS results generally showed agreement on the direction of the response of CBL to changes in the shape of the fiber. However, for the fiber described by { $r = 0.1$ , $m = 1$ , $θ = 0$ }, the FEM results suggested that CBL would decrease, while the NIGS predicted it would increase. Conversely, for { $r = 0.1$ , $m = 10$ , $θ = 3 π / 2$ }, the FEM results suggested an increase in CBL while the NIGS predicted a decrease.

To estimate w, in order to calculate γ for each finite element model, we integrated equation (2) using ode113 with $InitialStep = 10^{- 20}$ and $MaxStep = 10^{- 5}$ , as for the NIGS. The inconsistency values for the FEM results varied from $γ = 8 \times 10^{- 5}$ to $0.62$ , and were generally smaller for smaller disagreement between FEM and the NIGS. These higher inconsistency scores indicated that the FEM results did not correspond as closely to $β_{cr}$ as the NIGS results.

Computation times were measured using the tic and toc functions in MATLAB. To reduce computation time for the NIGS, we performed the NIGS with $MaxStep = 10^{- 3}$ to find a preliminary solution for β that then provided a starting value for the NIGS with $MaxStep = 10^{- 5}$ . This still gave a relatively long computation time of $159 \pm 8$ s, in contrast with computation times of $0.013 \pm 0.004$ , $0.12 \pm 0.03$ and $0.5 \pm 0.5 s$ for PSM Exact, FFT and Poly, respectively.

However, we found that if we used the NIGS with $MaxStep = 10^{- 4}$ after the preliminary search, this reduced the computation time to $18 \pm 1$ s across all fiber combinations, but increased γ only marginally (with values between $2 \times 10^{- 11}$ and $10^{- 8}$ ) but, importantly, differences in $β_{cr}$ predictions were less that with $10^{- 8}$ mg. Using the NIGS with $MaxStep = 10^{- 3}$ reduced computation times to $2.5 \pm 0.1$ s, but gave higher γ values (between between $2 \times 10^{- 11}$ and $9 \times 10^{- 8}$ ) and affected $β_{cr}$ predictions by less than 0.05 mg across all fiber combinations. Thus, predictions for $β_{cr}$ were less sensitive to integration tolerances than the predicted deformations.

Discussion

For the fibers considered, the buckling equation can be solved explicitly for only the simplest fiber (uniform along its length). Thus, in the absence of experimental verification (which is outside the scope of this paper), the inconsistency score, γ, described by equation (19) may provide a useful gauge for evaluation of the goodness of predictions. The first term in equation (19) tests for consistency of equation (2) without reference to boundary conditions, and would be negligible if w is obtained from $β_{cr}$ by numerical integration, irrespective of the accuracy of $β_{cr}$ . The second term tests that boundary conditions at $x = 1$ are met, and would be negligible for the PSMs, even if they yield the wrong $β_{cr}$ . Thus, both terms are required for robustly testing predictions.

Based on the γ scores, the NIGS method appeared to give more reliable predictions than the other methods considered, provided stringent error tolerances were observed. This came at the cost of having a much longer computation time (apart from the VIM). Relaxing the tolerances reduced the computation time without having a significant effect on the predicted CBL and, thus, is an option if faster prediction is required. This could have the consequence of increasing the inconsistency score (equation (19)) due to less accurate prediction of fiber deformation, but this would be a concern only if the accuracy of the results were in question.

The PSM Exact method, while giving good agreement with the NIGS method in many cases, was found to be unstable in others, and should be used with caution: it appears adequate for fibers/columns with small deviations from a uniform cylinder or for large deviations if the rate of change along the length is relatively low. This is certainly a scenario that is expected in keratin fibers.

The generally better agreement of PSM FFT with PSM Exact than PSM Poly for integer values of m reflects that a Fourier series representation is likely to be more robust than a polynomial representation when dealing with biological structures like wool, which often show seasonal (sinusoidal) variations in diameter. Conversely, the better agreement of PSM Poly with PSM Exact for half-sine-wave variations was expected as such a shape is more readily approximated by a polynomial, and since FFT performs better when the sampling window corresponds to integer cycles of periodic data. However, the PSMs in general appear to be unreliable because they can yield an unstable eigen problem, giving some erroneous results.

One would expect that if the FEM CBL predictions of He and Wang⁵ were correct, then our numerical integration of equation (2) would yield correct values for the corresponding fiber deformations, and thus the calculated γ scores would be negligible. On this basis, the FEM predictions appear to be less accurate than the NIGS predictions. It is unclear why and may be related to the choice of mesh shape and density. He⁴ reported using a mesh with 84 elements for the cross-section of the fiber and 120 layers along the fiber for all simulations performed. While this may have been too sparse, especially for cases when five or 10 undulations were modeled along the length of the fiber, the best agreement between the FEM and NIGS predictions were a fiber with 10 undulations of the largest diameter variation considered.

The inability of the VIM to converge was puzzling, but appears to be related to the choice of Lagrange multiplier. Applying variational iteration theory³⁰ to equation (8), λ must satisfy both the relationship

\frac{d^{4} λ}{d ξ^{4}} - 2 \frac{d^{3}}{d ξ^{3}} (\frac{ζ^{'} λ}{ζ}) + \frac{d^{2}}{d ξ^{2}} (λ \frac{ζ^{″} + β_{0}}{ζ}) = 0

(20)

and the boundary conditions arising from the restricted variation. The expression for λ (equation (11)) given by Coşkun and Atay¹² to cover general Euler buckling would imply from equation (20) that

2 \frac{d^{3}}{d ξ^{3}} (\frac{ζ' λ}{ζ}) = \frac{d^{2}}{d ξ^{2}} (λ \frac{ζ^{″} + β_{0}}{ζ})

, which would appear to restrict the functional forms of ζ for which the VIM is valid. Solving equation (20) to obtain a value of λ for an arbitrary functional form of ζ appears to be intractable, and thus should be done on a case-by-case basis. We note that applying the variational iteration theory to the form given in equation (8) does not yield a tractable solution.

The boundary conditions in equation (6) assume that a fiber protruding from a textile garment can be treated as being clamped by surrounding fibers at the garment end, that is, that it cannot move). Similarly, the boundary conditions in equation (7) assume that a fiber pushing against the skin can be treated as pinned, that is, that the end of the fiber is perhaps lodged in a pore and will not slide along the skin but can swivel. The accuracy of these assumptions has obvious implications for the accuracy of the predictions.

Conclusion

Our simulation results indicate how sinusoidal modulation of the fiber diameter can cause significant effects on its CBL, compared with a uniform fiber having the same mean diameter. This general conclusion was no surprise and therefore reassuring, because our aim was to compare the fine details of results from alternative methods. CBL was affected by the magnitude of the modulations, the number of modulations occurring along the fiber length and the position along the fiber where peak diameter occurred. When changes in diameter were relatively small and gradual along the length of the fiber, the resistance to end-on forces can be increased. However, for larger variations in diameter, increasing the magnitude of variation while retaining a gradual change along the fiber resulted in a weakening of the fiber. Interestingly, when the diameter at the clamped end equaled the mean fiber diameter, increasing the number of diameter modulations along the fiber showed a complex, diameter-dependent response. When the peak diameter was 10% or 30% of the mean diameter, the CBL decreased when the number of modulations in the fiber increased from one to five, but increased when the number of modulations in the fiber increased from five to 10. When the peak diameter was 50% of the mean, CBL increased monotonically with the number of modulations in the fiber.

This complex behavior even in highly simplified abstract structures (simple wave-form diameter changes) underscores the need for reliable computational methods to investigate CBLs when analyzing complicated biological structures, such as wool. The significantly poor performance of published methods (e.g. PSM and VIM) that had been “validated” against exact solutions for simple column geometries indicates the need for caution when using CBL predictions for irregular geometries. Our inconsistency analyses for the predictions of different methods indicate that the NIGS approach gave more reliable predictions than previously published results that had been obtained by FEM, and that we had originally intended to use as a standard.

We stress that this does not reflect on the FEM approach itself, but rather underscores that other approaches can be equally valid and that there is always a need to investigate the reliability of predictions, irrespective of the method used.

Wool fibers, in addition to being morphologically irregular, vary considerably in chemical and structural composition (and thus elasticity) along their lengths. Consideration of morphological and compositional irregularities together presents a more demanding computational problem as a future research challenge.

Footnotes

Acknowledgements

We would like to thank Peter Brorens, Andy Bray and Warren Meade for discussion at various stages of this work.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the New Zealand Ministry of Business, Innovation and Employment (contract number C10X0710) and the AgResearch Core funded Integrated Wool Research Science Programme.

References

Hearle

JWS

. A total model for the structural mechanics of wool. Wool Technol Sheep Breed 2003; 51: 95–117.

Chen

. Modelling and predicting textile behaviour, Cambridge: Woodside Publishing Limited, 2009.

Liu

Bryson

. A three-component model of the wool fibre - effects of morphology, elasticity and intermediate filament arrangement on fibre stiffness. J Text Inst 2002; 93: 121–131.

. The mechanical behaviour of irregular fibre materials, Ph.D. Thesis,: School of Engineering and Technology, Deakin University, 2002.

Wang

. Mechanical behavior of irregular fibers Part III: the flexural buckling behavior. Text Res J 2002; 72: 573–578.

Munro

Carnaby

. Wool fibre crimp, Part I: the effects of micro-fibrillar geometry. J Text Inst 1999; 90: 123–136.

Naylor

GRS

. The role of coarse fibres in fabric prickle using blended acrylic fibres of different diameters. Wool Technol Sheep Breed 1992; 40: 14–18.

Timoshenko

Young

. Elements of strength of materials, 5th ed. New York: Van Nostrand Reinhold Company, 1968.

Veitch

Naylor

GRS

. The mechanics of fibre buckling in relation to fabric-evoked prickle. Wool Technol Sheep Breed 1992; 40: 31–34.

10.

Asad

RAM

Zheng

Y-H

et al.

Characterization of prickle tactile discomfort properties of different textile single fibers using an axial fiber-compression-bending analyzer (FICBA). Text Res J 2015; 85: 512–523.

11.

Burkhart

Andrews

Dunning

. Finite element modeling mesh quality, energy balance and validation methods: A review with recommendations associated with the modeling of bone tissue. J Biomech 2013; 46: 1477–1488.

12.

Coşkun

Atay

. Determination of critical buckling load for elastic columns of constant and variable cross-sections using variational iteration method. Comput Math Appl 2009; 58: 2260–2266.

13.

Huang

Luo

Q-Z

. A simple method to determine the critical buckling loads for axially inhomogeneous beams with elastic restraint. Comput Math Appl 2011; 61: 2510–2517.

14.

Elishakoff

Rollot

. New closed-form solutions for buckling of a variable stiffness column by Mathematica®. J Sound Vib 1999; 224: 172–182.

15.

Huang

X-F

. An analytic approach for exactly determining critical loads of buckling of nonuniform columns. Int J Struct Stabil Dynam 2012; 12: 1250027-1–1250027-1.

16.

Cao

. Stability analysis of bars with varying cross-section. Int J Solid Struct 1995; 32: 3217–3228.

17.

Cao

. Stability analysis of a bar with multi-segments of varying cross-section. Comput Struct 1994; 53: 1085–1089.

18.

. Exact solutions for buckling of non-uniform columns under axial concentrated and distributed loading. Eur J Mech A/Solids 2001; 20: 485–500.

19.

Wang

Reddy

. Exact solutions for buckling of structural members. CRC Series in Computational Mechanics and Applied Analysis, Boca Raton, FL: CRC Press LLC, 2005.

20.

Totry

Altus

Proskura

. Buckling of non-uniform beams by a direct functional perturbation method. Probabil Eng Mech 2007; 22: 88–99.

21.

Afsharfard

Farshidianfar

. Finding the buckling load of non-uniform columns using the iteration perturbation method. Theor Appl Mech Lett 2014; 4: 041011–041011.

22.

Pinarbasi

. Lateral torsional buckling of rectangular beams using variational iteration method. Sci Res Essays 2011; 6: 1445–1457.

23.

Başbük

Atay

Erilmazza

. Numerical solutions for buckling loads of elastic columns with constant cross-section. Trans IoT Cloud Comput 2 2014; 2: 1–13.

24.

Katsikadelis

Tsiatas

. Non-linear dynamic analysis of beams with variable stiffness. J Sound Vib 2004; 270: 847–863.

25.

Yilmaz

Girgin

Evran

. Buckling analyses of axially functionally graded nonuniform columns with elastic restraint using a localized differential quadrature method. Math Probl Eng https://dx-doi-org.web.bisu.edu.cn/10.1155/2013/793062.

26.

Lee

Carr

Lee

et al.

Elasticas and buckling loads of shear deformable tapered columns. Int J Struct Stabil Dynam 2005; 5: 317–335.

27.

Lee

Carr

Lee

et al.

Buckling loads of columns with constant volume. J Sound Vib 2006; 294: 381–387.

28.

Lee

. Elastica and buckling load of simple tapered columns with constant volume. Int J Solid Struct 2000; 37: 2507–2518.

29.

Rychlewska

. A new approach for buckling analysis of axially functionally graded beams. J Appl Math Comput Mech 2015; 14: 95–102.

30.

J-H

. Variational iteration method - a kind of non-linear analytic technique: some examples. Int J Non Lin Mech 1999; 34: 699–708.