An inverse application of unit cells for extracting fibre properties from effective properties of composites

2.1. Complex function solution of the GPE problem

For the conventional plane strain problem in the cross-section plane transverse to the fibre direction, and the anticlastic problem for the shear deformation along fibres out of the cross-section plane, there have been elegant solutions separately presented in the theory of elasticity employing complex functions [16,17]. The GPE problem represents a slight generalisation by merging both into the same formulation while allowing non-vanishing constant direct strain in the fibre direction [15].

Introduce the coordinate system in such a way that the cross-section of the composite lies in the x-y plane while the fibres run along the z-axis. The coordinate plane x-y coincides with the complex plane on which any point is defined by a complex variable ζ = x + iy. The solution to any GPE problem can be given in terms of three analytic complex functions, φ(ζ), ψ(ζ), and ω(ζ), also called potentials since they are all analytic and hence bear a path-independent characteristic, through which the solution can always be expressed as:

σ x + σ y = 4 Re φ ′ ( ζ ) σ y − σ x + 2 i τ xy = 2 ( ζ ¯ φ ′ ′ ( ζ ) + ψ ′ ( ζ ) ) 2 G T ( u + iv ) = κ φ ( ζ ) − ζ φ ′ ( ζ ) ¯ − ψ ( ζ ) ¯ − 2 G T ν L ε 0 ζ

(1a)

where σ x , σ y , and τ xy are the stresses, and u and v are the displacements in the x-y plane, and ε 0 is the direct strain in the fibre direction which is constant over the cross-section plane [16]; and

τ xz − i τ yz = ω ′ ( ζ ) G L w = Re ω ( ζ ) + G L ε 0

(1b)

where τ xz and τ yz are the stresses and w is the displacement out of the x-y plane [17]. As a GPE problem, the remaining stress component can always be obtained as follows:

σ z = ν L ( σ x + σ y ) + E L ε 0 .

(2)

In equations (1) and (2), E_L, E_T, G_L, G_T, ν_L, and ν_T are the longitudinal and transverse Young’s moduli, shear moduli, and Poisson’s ratios, respectively, where the transverse properties satisfy,

G T = E T 2 ( 1 + ν T )

(3)

as a requirement from transverse isotropy which most UD composites feature;

κ = − 1 + 4 ( 1 − ν L 2 E T / E L ) ( 1 + ν T )

(4)

is a material constant; i = − 1 is the imaginary unity; Re represents the real part of a complex function while the corresponding imaginary part is denoted as Im as will be employed later. In equation (1), the length of the composite in the z-direction has been assumed to be unity. A prime on a function represents its derivative with respect to the complex variable ζ. A bar on the top of a function means its complex conjugate.

From equation (1), it is straightforward to obtain the individual components of the stresses and displacements and expressed in terms of the three complex potentials as follows:

σ x = Re ( 2 φ ′ ( ζ ) − ζ ¯ φ ′ ′ ( ζ ) − ψ ′ ( ζ ) ) σ y = Re ( 2 φ ′ ( ζ ) + ζ ¯ φ ′ ′ ( ζ ) + ψ ′ ( ζ ) ) τ xy = Im ( ζ ¯ φ ′ ′ ( ζ ) + ψ ′ ( ζ ) )

(5a)

τ xz = Re ω ′ ( ζ ) τ yz = − Im ω ′ ( ζ )

(5b)

in addition to the longitudinal direct stress σ z as given in equation (2). The corresponding displacements are:

u = 1 2 G T Re ( κ φ ( ζ ) − ζ φ ′ ( ζ ) ¯ − ψ ( ζ ) ¯ ) − ν L ε 0 x v = 1 2 G T Im ( κ φ ( ζ ) − ζ φ ′ ( ζ ) ¯ − ψ ( ζ ) ¯ ) − ν L ε 0 y

(6a)

w = 1 G L Re ω ( ζ ) + ε 0 .

(6b)

It is now clear that the solution to the GPE problem relies on three complex potentials, with which all governing equations of elasticity are satisfied. All they need to meet are the boundary conditions, including those on the internal boundaries and the external ones.

2.2. Series form of the complex potentials

To model the behaviour of the composites using a unit cell, whether square one or hexagonal one, as shown in Figure 1, there is one fibre involved in each unit cell surrounded by matrix. Thus, there are two regions in the domain, corresponding to the fibre and the matrix, respectively. For simplicity, the problem can always be scaled appropriately, so that the fibre region in the domain is defined as a unit circle centred at the origin of the x-y plane, that is, ζ = ρ = e i θ , with θ being the polar angle. Complex functions have to be defined separately within these two distinct regions. While the same notation will be adopted, each will be given a subscript of either f or m representing fibre and matrix, respectively. Solution in a form of series will be pursued below.

OPEN IN VIEWER

Figure 1.

(a) A square unit cell, and (b) a hexagonal unit cell.

The complex potentials for the fibre, that is, defined within the unit circle, are φ_f, ψ_f, and ω_f, expressed in terms of Taylor series as:

φ f = ∑ k = 1 ∞ a k ζ k , ψ f = ∑ k = 1 ∞ b k ζ k , and ω f = ∑ k = 1 ∞ c k ζ k

(7)

where the zeroth-order terms have been dropped as the rigid body displacements are constrained at the origin, that is,

u | ζ = 0 = v | ζ = 0 = w | ζ = 0 = 0 or φ f ( 0 ) = ψ f ( 0 ) = ω f ( 0 ) = 0 .

(8)

In the matrix region, that is, within the unit cell but outside the unit circle, the complex potentials for the matrix, φ_m, ψ_m, and ω_m, can be expressed in terms of Laurent series as:

φ m = ∑ k = 0 ∞ A k ζ k + ∑ k = 1 ∞ D k ζ − k , ψ m = ∑ k = 0 ∞ B k ζ k + ∑ k = 1 ∞ F k ζ − k , and ω m = ∑ k = 0 ∞ C k ζ k + ∑ k = 1 ∞ H k ζ − k .

(9)

The rigid body displacements should not be constrained in the matrix directly. Instead, their constraints will be accomplished through the continuity conditions between the fibre and the matrix. For this reason, the zeroth-order terms have been retained in equation (9).

It is well established that Taylor series in a solid domain and Laurent series in a hollow domain as involved in the current problem are both analytic. The validity of the solution can therefore be taken for granted once the coefficients in series (7) and (9) have been correctly determined. Apparently, coefficients a_k, b_k, and c_k in equation (7) are associated with the fibre, and A_k, B_k, C_k, D_k, F_k, and H_k in equation (9) with the matrix. They should be determined by the continuity conditions between fibre and matrix, and the boundary conditions for the unit cell concerned. It will be shown below that those associated with the matrix can be expressed in terms of those associated with the fibre after imposing the continuity conditions between the fibre and the matrix across the interface.

It should be pointed out that, unlike the unit cells analysed using FEM where traction boundary conditions are natural boundary conditions and do not need any attention from users, both displacement boundary conditions and the traction boundary conditions will have to be satisfied simultaneously as the present solution is no longer based on a variational principle.

2.3. Continuity across fibre–matrix interface

It is assumed that perfect bonding is present between the fibre and the matrix. Thus, both the displacements and traction components should be continuous at interface ζ = ρ . Denote the traction components in the x-y plane as:

X = σ x cos β + τ xy sin β Y = τ xy cos β + σ y sin β

(10)

where β is the angle between the outward normal to the unit circle and the x-axis. They are related to the complex potential functions in general as follows [16]:

i ∫ ( X + iY ) ds = φ ( ρ ) + ρ φ ′ ( ρ ) ¯ + ψ ( ρ ) ¯

(11)

where s is the arc length along the interface. The traction component in the z-direction can be given as [12]:

Z = τ xz cos β + τ yz sin β = ρ ω ′ ( ρ ) + ρ ω ′ ( ρ ) ¯ .

(12)

The traction continuity conditions between fibre and matrix across the interface can be obtained as:

φ m ( ρ ) + ρ φ ′ m ( ρ ) ¯ + ψ m ( ρ ) ¯ = φ f ( ρ ) + ρ φ ′ f ( ρ ) ¯ + ψ f ( ρ ) ¯ ρ ω ′ m ( ρ ) + ρ ω ′ m ( ρ ) ¯ = ρ ω ′ f ( ρ ) + ρ ω ′ f ( ρ ) ¯ .

(13a)

Expressing displacement (6) in the fibre and the matrix, respectively, the displacement continuity conditions across the interface are given as:

κ m φ m ( ρ ) − ρ φ ′ m ( ρ ) ¯ − ψ m ( ρ ) ¯ = η ( κ f φ f ( ρ ) − ρ φ ′ f ( ρ ) ¯ − ψ f ( ρ ) ¯ ) + 2 G Tm ( ν Lm − ν Lf ) ε 0 ρ ω m ( ρ ) + ω m ( ρ ) ¯ = λ ( ω f ( ρ ) + ω f ( ρ ) ¯ )

(13b)

where κ_f and κ_m are the same expression as given in equation (4) but in the fibre and the matrix, respectively

η = G Tm G Tf and λ = G Lm G Lf .

(14)

The term associated with ε 0 also appears in displacement w. However, it does not show in the second equation of equation (13b) because it appears on both sides of the equation with equal magnitude as it corresponds to the same z-coordinate and therefore is duly cancelled.

Substitute equations (7) and (9) into equation (13), bearing in mind the following equalities from the unit circle:

ρ = e i θ , ρ k ¯ = e − ik θ = ( e ik θ ) − 1 = ( ρ k ) − 1 = ρ − k and ρ − k ¯ = e − ik θ ¯ = e ik θ = ρ k

(15)

the coefficients A_k, B_k, C_k, D_k, F_k, and H_k in equation (9) can be linearly expressed in terms of a_k, b_k, and c_k in equation (7), after comparing the coefficients to the terms of the same order on both sides of the equations:

A 0 = 2 ( 1 − η ) 1 + κ m a ¯ 2 , A 1 = 1 1 + κ m ( a ¯ 1 + a 1 − η ( a ¯ 1 − κ f a 1 ) ) + 2 G Tm ( ν Lm − ν Lf ) 1 + κ m ε 0

A k = 1 + η κ f 1 + κ m a k ( k = 2 , 3 , … … . )

B 0 = 2 ( κ m + η − 1 − η κ f ) 1 + κ m a 2

B k = 1 1 + κ m ( ( k + 2 ) ( κ m + η − 1 − η κ f ) a k + 2 + ( κ m + η ) b k ) ( k = 1 , 2 , 3 , . . . . . . . )

(16)

D k = 1 − η 1 + κ m ( ( k + 2 ) a ¯ k + 2 + b ¯ k ) ( k = 1 , 2 , 3 , . . . . . . . )

F 1 = 1 1 + κ m ( ( κ m − 1 + η − η κ f ) ( a ¯ 1 + a 1 ) ) − 4 G Tm ( ν Lm − ν Lf ) 1 + κ m ε 0

F k = ( κ m − η κ f + k ( k − 2 ) ( 1 − η ) ) a ¯ k + ( k − 2 ) ( 1 − η ) b ¯ k − 2 1 + κ m ( k = 2 , 3 , . . . . . . . )

C 0 + C ¯ 0 = 0 hence C 0 = 0

C k = ( 1 + λ ) c k 2 and H k = − ( 1 − λ ) c ¯ k 2 ( k = 1 , 2 , . . . . . . . ) .

In the above, C 0 + C ¯ 0 = 0 only means Re C 0 = 0 . It does not place any constraint on the ImC₀. From equation (6b), it is clear that displacement w is only related to the real part of ω. The imaginary part of ω is of no physical significance at all. As a result, it can take any value arbitrarily without affecting the solution. For convenience, take ImC₀ and hence C 0 = 0 .

Substituting the expressions of A_k, B_k, C_k, D_k, F_k and H_k as obtained in equation (16) into series (9), those potentials for the matrix can be expressed in terms of a_k, b_k, and c_k as follows:

φ m = 2 G Tm ( ν Lm − ν Lf ) 1 + κ m ε 0 ζ + 1 1 + κ m ( a ¯ 1 + a 1 − η ( a ¯ 1 − κ f a 1 ) ) ζ + 1 + η κ f 1 + κ m a 2 ζ 2 + 2 ( 1 − η ) 1 + κ m a ¯ 2 + 1 + η κ f 1 + κ m ∑ k = 3 ∞ a k ζ k + 1 − η 1 + κ m ∑ k = 3 ∞ k a ¯ k ζ − ( k − 2 ) + 1 − η 1 + κ m ∑ k = 1 ∞ b ¯ k ζ − k

ψ m = − 4 G Tm ( ν Lm − ν Lf ) 1 + κ m ζ − 1 ε 0 + 1 1 + κ m ( κ m − 1 + η − η κ f ) ζ − 1 ( a ¯ 1 + a 1 ) + 2 1 + κ m ( κ m + η − 1 − η κ f ) a 2 + 1 1 + κ m ∑ k = 1 ∞ ζ k ( ( k + 2 ) ( κ m + η − 1 − η κ f ) a ¯ k + 2 + ( κ m + η ) b ¯ k ) + 1 1 + κ m ∑ k = 2 ∞ ζ − k ( κ m − η κ f + k ( k − 2 ) ( 1 − η ) ) a ¯ k + ( k − 2 ) ( 1 − η ) b ¯ k − 2

ω m = 1 2 ∑ k = 1 ∞ ( ( λ + 1 ) c k ζ k + ( λ − 1 ) c ¯ k ζ − k ) .

(17)

In general, between the two phases of different materials, the continuity conditions in displacement and traction will be satisfied as long as series (7) and (9) are truncated in a consistent manner. This can be easily demonstrated if one substitutes equations (7) and (17) into equation (13).

2.4. Relative displacement boundary condition and periodic traction boundary conditions

After satisfying the continuity conditions between the fibre and matrix, coefficients a_k, b_k, and c_k become the only set of unknowns to be determined. The conditions available are the boundary conditions for the unit cell. According to the study by Li and Sitnikova [5,6], they are the relative displacement boundary conditions and the periodic traction boundary conditions. Because the traction boundary conditions involve traction on different parts of the boundary in the same equation, the expression of traction as given in equation (11) is no longer convenient to apply. The periodic traction boundary conditions can be given as [6]:

[ σ x Δ τ xy Δ τ xz Δ τ xy Δ σ y Δ τ yz Δ τ xz Δ τ yz Δ σ z Δ ] { cos β sin β 0 } = [ σ x τ xy τ xz τ xy σ y τ yz τ xz τ yz σ z ] { cos β sin β 0 } .

(18a)

In a GPE problem, the z component of the outward normal to the boundary always vanishes. The outward normal at paired points on the opposite side of the unit cell is the same but opposite in direction, as reflected in equation (18a). Stresses σ_x, σ_y, σ_z, τ_yz, τ_xz and τ_xy, and σ_xΔ, σ_yΔ, σ_zΔ, τ_yzΔ, τ_xzΔ and τ_xyΔ are at paired points, (x, y) and (x_Δ, y_Δ), on the boundary.

The relative displacement boundary conditions for the unit cell are [5,6]:

{ u Δ v Δ w Δ } − { u v w } = [ ε x 0 0 0 γ xy 0 ε y 0 0 γ xz 0 γ yz 0 ε 0 ] { Δ x Δ y Δ z } = [ ε x 0 0 0 γ xy 0 ε y 0 0 γ xz 0 γ yz 0 ε 0 ] { cos β sin β 0 } Δ = [ ε x 0 0 γ xy 0 ε y 0 γ xz 0 γ yz 0 ] { cos β sin β } Δ

(18b)

where u, v and w, and u_Δ, v_Δ and w_Δ are the displacements at the paired points, Δ is the distance between the paired points, that is the amount of translation involved in the formulation of the unit cell, with Δ _x , Δ _y and Δ _z as its components, and ε x 0 , ε y 0 , ε 0 , γ yz 0 γ xz 0 and γ xy 0 are the average strains in the unit cell. Since the two paired points are always on the same cross-section of the material, they share the common z-coordinate and hence Δ _z  = 0. As a result, ε 0 will not appear in the actual relative displacement boundary conditions.

To prescribe the periodic traction boundary conditions, all stress components can be expressed in terms of the complex potentials as follows:

σ x Δ = Re ( 2 φ ′ m ( ζ Δ ) − ζ ¯ Δ φ ′ ′ m ( ζ Δ ) − ψ ′ m ( ζ Δ ) ) σ x = Re ( 2 φ ′ m ( ζ ) − ζ ¯ φ ′ ′ m ( ζ ) − ψ ′ m ( ζ ) ) σ y Δ = Re ( 2 φ ′ m ( ζ Δ ) + ζ ¯ Δ φ ′ ′ m ( ζ Δ ) + ψ ′ m ( ζ Δ ) ) σ y = Re ( 2 φ ′ m ( ζ ) + ζ ¯ φ ′ ′ m ( ζ ) + ψ ′ m ( ζ ) ) τ xy Δ = Im ( ζ ¯ Δ φ ′ ′ m ( ζ Δ ) + ψ ′ m ( ζ Δ ) ) τ xy = Im ( ζ ¯ φ ′ ′ m ( ζ ) + ψ ′ m ( ζ ) )

τ xz Δ = Re ω ′ m ( ζ Δ ) τ xz = Re ω ′ m ( ζ ) τ yz Δ = − Im ω ′ m ( ζ Δ ) τ yz = − Im ω ′ m ( ζ )

(19)

where all complex potentials involved are those associated with the matrix as the boundary of the unit cell lies completely in the matrix.

The following relationships have been used repeatedly in the derivation:

ζ = | ζ | ρ = | ζ | e i θ , ζ ¯ = | ζ | ρ − 1 = | ζ | e − i θ , ζ k = | ζ | k ρ k = | ζ | k e ik θ , ζ k ¯ = | ζ | k ρ − k = | ζ | k e − ik θ , ζ − k = | ζ | − k ρ − k = | ζ | − k e − ik θ , and ζ − k ¯ = | ζ | − k ρ k = | ζ | − k e ik θ .

(20)

In the present GPE problem, the part within the plane of the cross-section and the part out of this plane are fully decoupled. They can be addressed separately for both clarity and computational efficiency. The part in the cross-section plane will be considered first below. The terms involved can be given in the form of series containing unknown coefficients a_k, b_k, and c_k and their conjugates only as follows:

φ ′ m ( ζ ) = 1 1 + κ m ( 2 G Tm ( ν Lm − ν Lf ) ε 0 + ( a ¯ 1 + a 1 − η ( a ¯ 1 − κ f a 1 ) ) + 2 ( 1 + η κ f ) a 2 ζ + ∑ k = 3 ∞ k ( 1 + η κ f ) a k ζ k − 1 + ∑ k = 3 ∞ − k ( k − 2 ) ( 1 − η ) a ¯ k ζ − ( k − 1 ) + ∑ k = 1 ∞ − k ( 1 − η ) b ¯ k ζ − ( k + 1 ) )

ζ ¯ φ ′ ′ m ( ζ ) = 1 1 + κ m ( 2 ( 1 + η κ f ) a 2 | ζ | ρ − 1 + ∑ k = 3 ∞ k ( k − 1 ) ( 1 + η κ f ) a k | ζ | k − 1 ρ k − 3 + ∑ k = 3 ∞ k ( k − 1 ) ( k − 2 ) ( 1 − η ) a ¯ k | ζ | − ( k − 1 ) ρ − ( k + 1 ) + ∑ k = 1 ∞ k ( k + 1 ) ( 1 − η ) b ¯ k | ζ | − ( k + 1 ) ρ − ( k + 3 ) )

ψ ′ m ( ζ ) = 1 1 + κ m ( 4 G Tm ( ν Lm − ν Lf ) ζ − 2 ε 0 − ( κ m − 1 + η − η κ f ) ζ − 2 ( a ¯ 1 + a 1 ) − 2 ( κ m − η κ f ) ζ − 3 a ¯ 2 + ∑ k = 3 ∞ k ( k − 2 ) ( κ m + η − 1 − η κ f ) ζ k − 3 a k − ∑ k = 3 ∞ k ( κ m − η κ f + k ( k − 2 ) ( 1 − η ) ) ζ − ( k + 1 ) a ¯ k + ∑ k = 1 ∞ k ( κ m + η ) ζ k − 1 b k − ∑ k = 1 ∞ k ( k + 2 ) ( 1 − η ) ζ − ( k + 3 ) b ¯ k ) .

(21)

Replacing ζ with ζ_Δ, one obtains φ ′ m ( ζ Δ ) , ζ ¯ Δ φ ′ ′ m ( ζ Δ ) , and ψ ′ m ( ζ Δ ) accordingly, where ζ and ζ_Δ are the complex coordinates of the two paired points (x, y) and (x_Δ, y_Δ). Then, the periodic traction boundary conditions can be expressed as:

Re ( 2 φ ′ m ( ζ Δ ) − ζ ¯ Δ φ ′ ′ m ( ζ Δ ) − ψ ′ m ( ζ Δ ) ) cos β + Im ( ζ ¯ Δ φ ′ ′ m ( ζ Δ ) + ψ ′ m ( ζ Δ ) ) sin β = Re ( 2 φ ′ m ( ζ ) − ζ ¯ φ ′ ′ m ( ζ ) − ψ ′ m ( ζ ) ) cos β + Im ( ζ ¯ φ ′ ′ m ( ζ ) + ψ ′ m ( ζ ) ) sin β Im ( ζ ¯ Δ φ ′ ′ m ( ζ Δ ) + ψ ′ m ( ζ Δ ) ) cos β + Re ( 2 φ ′ m ( ζ Δ ) + ζ ¯ Δ φ ′ ′ m ( ζ Δ ) + ψ ′ m ( ζ Δ ) ) sin β = Im ( ζ ¯ φ ′ ′ m ( ζ ) + ψ ′ m ( ζ ) ) cos β + Re ( 2 φ ′ m ( ζ ) + ζ ¯ φ ′ ′ m ( ζ ) + ψ ′ m ( ζ ) ) sin β .

(22a)

Since they are paired, with the origin of the coordinate system being placed at the centre of the unit cell, the polar radii of these two points are always the same, that is,

| ζ Δ | = | ζ | .

(23)

Substituting equation (21) and their counterparts associated with ζ_Δ into equation (22a), the periodic traction boundary conditions in the cross-section plane become:

( 2 Re Q 1 − Re Q 2 ) cos β + Im Q 2 sin β = − Re R cos β + Im R sin β ( 2 Re Q 1 + Re Q 2 ) sin β + Im Q 2 cos β = Re R sin β + Im R cos β

(22b)

where

Q 1 = ( 2 ( 1 + η κ f ) | ζ | ( ρ Δ − ρ ) a 2 + ∑ k = 3 ∞ k ( 1 + η κ f ) | ζ | k − 1 ( ρ Δ k − 1 − ρ k − 1 ) a k − ∑ k = 3 ∞ k ( k − 2 ) ( 1 − η ) | ζ | − ( k − 1 ) ( ρ Δ − ( k − 1 ) − ρ − ( k − 1 ) ) a ¯ k − ∑ k = 1 ∞ k ( 1 − η ) | ζ | − ( k + 1 ) ( ρ Δ − ( k + 1 ) − ρ − ( k + 1 ) ) b ¯ k )

(24a)

Q 2 = ( − ( κ m − 1 + η − η κ f ) | ζ | − 2 ( ρ Δ − 2 − ρ − 2 ) a 1 − ( κ m − 1 + η − η κ f ) | ζ | − 2 ( ρ Δ − 2 − ρ − 2 ) a ¯ 1 + 2 ( 1 + η κ f ) | ζ | ( ρ Δ − 1 − ρ − 1 ) a 2 − 2 ( κ m − η κ f ) | ζ | − 3 ( ρ Δ − 3 − ρ − 3 ) a ¯ 2 + ∑ k = 3 ∞ k ( ( k − 1 ) ( 1 + η κ f ) | ζ | k − 1 + ( k − 2 ) ( κ m + η − 1 − η κ f ) | ζ | k − 3 ) ( ρ Δ k − 3 − ρ k − 3 ) a k + ∑ k = 3 ∞ k ( ( k − 1 ) ( k − 2 ) ( 1 − η ) | ζ | − ( k − 1 ) − ( κ m − η κ f + k ( k − 2 ) ( 1 − η ) ) | ζ | − ( k + 1 ) ) ( ρ Δ − ( k + 1 ) − ρ − ( k + 1 ) ) a ¯ k + ∑ k = 1 ∞ k ( κ m + η ) | ζ | k − 1 ( ρ Δ k − 1 − ρ k − 1 ) b k + ∑ k = 1 ∞ k ( 1 − η ) ( ( k + 1 ) | ζ | − ( k + 1 ) − ( k + 2 ) | ζ | − ( k + 3 ) ) ( ρ Δ − ( k + 3 ) − ρ − ( k + 3 ) ) b ¯ k )

(24b)

R = − 4 G Tm ( ν Lm − ν Lf ) | ζ | − 2 ( ρ Δ − 2 − ρ − 2 ) ε 0 .

(24c)

Equation (23) is the periodic traction boundary conditions for the unit cell. They are simultaneous equations for the unknowns a_k, b_k, and c_k and their conjugates, with expressions Q₁, Q₂ and R involved given in equation (24), respectively.

Now, attention can be moved to the relative displacement boundary conditions. According to equation (6a), displacements in the cross-section plane can be expressed in terms of the following terms:

φ m ( ζ ) = 1 1 + κ m ( 2 G Tm ( ν Lm − ν Lf ) ε 0 | ζ | ρ + ( a 1 + a ¯ 1 + η κ f a 1 − η a ¯ 1 ) | ζ | ρ + ( 1 + η κ f ) a 2 | ζ | 2 ρ 2 + 2 ( 1 − η ) a ¯ 2 + ∑ k = 3 ∞ ( 1 + η κ f ) a k | ζ | k ρ k + ∑ k = 3 ∞ k ( 1 − η ) a ¯ k | ζ | − ( k − 2 ) ρ − ( k − 2 ) + ∑ k = 1 ∞ ( 1 − η ) b ¯ k | ζ | − k ρ − k )

(25a)

ζ φ ¯ ′ m ( ζ ) = 1 1 + κ m ( 2 G Tm ( ν Lm − ν Lf ) | ζ | ρ ε 0 + ( a 1 + a ¯ 1 + η κ f a ¯ 1 − η a 1 ) | ζ | ρ + 2 ( 1 + η κ f ) | ζ | 2 a ¯ 2 + ∑ k = 3 ∞ k ( 1 + η κ f ) | ζ | k ρ − ( k − 2 ) a ¯ k − ∑ k = 3 ∞ k ( k − 2 ) ( 1 − η ) | ζ | − ( k − 2 ) ρ k a k − ∑ k = 1 ∞ k ( 1 − η ) | ζ | − k ρ k + 2 b k )

(25b)

ψ ¯ m ( ζ ) = 1 1 + κ m ( − 4 G Tm ( ν Lm − ν Lf ) | ζ | − 1 ρ ε 0 + ( κ m − 1 + η − η κ f ) | ζ | − 1 ρ ( a ¯ 1 + a 1 ) + 2 ( κ m + η − 1 − η κ f ) a ¯ 2 + ( κ m − η κ f ) | ζ | − 2 ρ 2 a 2 + ∑ k = 3 ∞ k ( κ m + η − 1 − η κ f ) | ζ | k − 2 ρ − ( k − 2 ) a ¯ k + ∑ k = 3 ∞ ( κ m − η κ f + k ( k − 2 ) ( 1 − η ) ) | ζ | − k ρ k a k + ∑ k = 1 ∞ ( κ m + η ) | ζ | k ρ − k b ¯ k + ∑ k = 1 ∞ k ( 1 − η ) | ζ | − ( k + 2 ) ρ k + 2 b k ) .

(25c)

Replacing ζ with ζ_Δ, one obtains φ m ( ζ Δ ) , ζ Δ φ ¯ ′ m ( ζ Δ ) , and ψ ¯ m ( ζ Δ ) accordingly. Then, the relative displacement boundary conditions in the cross-section plane can be obtained as: