Chapter 2: Fiber-Reinforced Composite Materials

2.1 Types of Fiber-Reinforced Composite Materials

2.1.1 CFRP (Carbon Fiber Reinforced Plastic)

A composite material using carbon fiber as reinforcement, widely used in aerospace, automotive, and sporting goods.

2.1.2 GFRP (Glass Fiber Reinforced Plastic)

A composite material using glass fiber as reinforcement, offering excellent cost performance.

2.2 Mechanics of a Single Layer (Lamina)

2.2.1 Stress-Strain Relations

The stress-strain relationship for a unidirectional fiber-reinforced single layer in the principal coordinate system (1-2 axes, fiber direction as 1-axis) is expressed as follows:

where the components of the reduced stiffness matrix \([Q]\) are:

\(E_1, E_2\): longitudinal and transverse moduli, \(\nu_{12}, \nu_{21}\): Poisson's ratios, \(G_{12}\): shear modulus
(Reciprocity: \(\nu_{12}/E_1 = \nu_{21}/E_2\))

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0, <2.0.0

import numpy as np

def lamina_stiffness_matrix(E1, E2, nu12, G12):
    """
    Calculate reduced stiffness matrix [Q] for a single layer

    Parameters:
    -----------
    E1 : float
        Longitudinal modulus [GPa]
    E2 : float
        Transverse modulus [GPa]
    nu12 : float
        Major Poisson's ratio
    G12 : float
        Shear modulus [GPa]

    Returns:
    --------
    Q : ndarray (3x3)
        Reduced stiffness matrix [GPa]
    """
    # Calculate nu21 from reciprocity
    nu21 = nu12 * E2 / E1

    # Q matrix components
    denom = 1 - nu12 * nu21
    Q11 = E1 / denom
    Q22 = E2 / denom
    Q12 = nu12 * E2 / denom
    Q66 = G12

    Q = np.array([
        [Q11, Q12, 0],
        [Q12, Q22, 0],
        [0, 0, Q66]
    ])

    return Q

# CFRP/Epoxy single layer material properties
E1 = 140.0   # GPa
E2 = 10.0    # GPa
nu12 = 0.30
G12 = 5.0    # GPa

Q = lamina_stiffness_matrix(E1, E2, nu12, G12)

print("Reduced stiffness matrix [Q] (GPa):")
print(Q)
print(f"\nQ11 = {Q[0,0]:.2f} GPa")
print(f"Q22 = {Q[1,1]:.2f} GPa")
print(f"Q12 = {Q[0,1]:.2f} GPa")
print(f"Q66 = {Q[2,2]:.2f} GPa")

# Stress calculation example: ε1=0.005, ε2=0, γ12=0
strain = np.array([0.005, 0, 0])
stress = Q @ strain

print(f"\nStrain [ε1, ε2, γ12]: {strain}")
print(f"Stress [σ1, σ2, τ12] (MPa): {stress}")
print(f"Fiber direction stress σ1: {stress[0]:.1f} MPa")

2.2.2 Coordinate Transformation (Off-Axis Loading)

When fibers make an angle \(\theta\) with the loading direction, coordinate transformation is required. The stiffness matrix \([\bar{Q}]\) in the global coordinate system (x-y) is:

The transformation matrix \([T]\) is:

where \(c = \cos\theta\), \(s = \sin\theta\)

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0

import numpy as np
import matplotlib.pyplot as plt

def transformation_matrix(theta_deg):
    """
    Stress-strain transformation matrix [T]

    Parameters:
    -----------
    theta_deg : float
        Fiber orientation angle [degrees]

    Returns:
    --------
    T : ndarray (3x3)
        Transformation matrix
    """
    theta = np.radians(theta_deg)
    c = np.cos(theta)
    s = np.sin(theta)

    T = np.array([
        [c**2, s**2, 2*s*c],
        [s**2, c**2, -2*s*c],
        [-s*c, s*c, c**2 - s**2]
    ])

    return T

def offaxis_stiffness(Q, theta_deg):
    """
    Calculate off-axis stiffness matrix [Q̄]

    Parameters:
    -----------
    Q : ndarray (3x3)
        Stiffness matrix in principal coordinate system
    theta_deg : float
        Fiber orientation angle [degrees]

    Returns:
    --------
    Q_bar : ndarray (3x3)
        Off-axis stiffness matrix
    """
    T = transformation_matrix(theta_deg)
    T_inv = np.linalg.inv(T)

    Q_bar = T_inv @ Q @ T_inv.T

    return Q_bar

# CFRP material properties
E1 = 140.0
E2 = 10.0
nu12 = 0.30
G12 = 5.0

Q = lamina_stiffness_matrix(E1, E2, nu12, G12)

# Range of orientation angles
angles = np.arange(0, 91, 5)
Q_bar_11 = []
Q_bar_22 = []
Q_bar_66 = []
Q_bar_16 = []  # Coupling term

for theta in angles:
    Q_bar = offaxis_stiffness(Q, theta)
    Q_bar_11.append(Q_bar[0, 0])
    Q_bar_22.append(Q_bar[1, 1])
    Q_bar_66.append(Q_bar[2, 2])
    Q_bar_16.append(Q_bar[0, 2])

# Visualization
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Stiffness components
ax1.plot(angles, Q_bar_11, 'b-', linewidth=2, label='Q̄₁₁ (x-direction)')
ax1.plot(angles, Q_bar_22, 'r-', linewidth=2, label='Q̄₂₂ (y-direction)')
ax1.plot(angles, Q_bar_66, 'g-', linewidth=2, label='Q̄₆₆ (shear)')
ax1.set_xlabel('Fiber Orientation Angle θ [degrees]')
ax1.set_ylabel('Stiffness Component [GPa]')
ax1.set_title('Angular Dependence of Off-Axis Stiffness')
ax1.grid(True, alpha=0.3)
ax1.legend()

# Coupling term
ax2.plot(angles, Q_bar_16, 'm-', linewidth=2)
ax2.axhline(y=0, color='k', linestyle='--', alpha=0.3)
ax2.set_xlabel('Fiber Orientation Angle θ [degrees]')
ax2.set_ylabel('Q̄₁₆ [GPa]')
ax2.set_title('Extension-Shear Coupling')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('offaxis_stiffness.png', dpi=300, bbox_inches='tight')
plt.close()

print(f"0 degrees (fiber direction): Q̄11 = {Q_bar_11[0]:.1f} GPa")
print(f"45 degrees: Q̄11 = {Q_bar_11[9]:.1f} GPa, Q̄16 = {Q_bar_16[9]:.1f} GPa")
print(f"90 degrees (perpendicular to fiber): Q̄11 = {Q_bar_11[-1]:.1f} GPa")

2.3 Classical Laminate Theory (CLT)

2.3.1 Assumptions for Laminates

Classical Laminate Theory (CLT) is based on the following assumptions:

• Each layer is perfectly bonded (no interlaminar slip)
• Normal stresses in the thickness direction are negligible (plane stress state)
• Kirchhoff-Love assumption: normals to the midplane remain normal after deformation

Thus, the in-plane strains at distance \(z\) from the midplane are:

\(\epsilon^0\): midplane strains, \(\kappa\): curvatures

2.3.2 A-B-D Matrices

The relationship between resultant forces and moments and midplane strains and curvatures:

Physical meaning of each matrix:

• [A]: Extensional stiffness matrix
• [B]: Bending-extension coupling matrix
• [D]: Bending stiffness matrix

Calculation formulas:

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0, <2.0.0

import numpy as np

class Laminate:
    """Laminate analysis class"""

    def __init__(self, Q, layup, t_ply):
        """
        Parameters:
        -----------
        Q : ndarray (3x3)
            Single layer principal stiffness matrix
        layup : list of float
            Laminate configuration [θ1, θ2, ..., θn] (degrees)
        t_ply : float
            Single layer thickness [mm]
        """
        self.Q = Q
        self.layup = layup
        self.t_ply = t_ply
        self.n_plies = len(layup)
        self.total_thickness = self.n_plies * t_ply

        # Calculate z coordinates (referenced to midplane)
        self.z = np.linspace(
            -self.total_thickness/2,
            self.total_thickness/2,
            self.n_plies + 1
        )

    def compute_ABD(self):
        """Calculate A-B-D matrices"""
        A = np.zeros((3, 3))
        B = np.zeros((3, 3))
        D = np.zeros((3, 3))

        for k in range(self.n_plies):
            # Q̄ matrix for each layer
            Q_bar = offaxis_stiffness(self.Q, self.layup[k])

            # Layer thickness coordinates
            z_k = self.z[k]
            z_k1 = self.z[k + 1]

            # Calculate A, B, D
            A += Q_bar * (z_k1 - z_k)
            B += 0.5 * Q_bar * (z_k1**2 - z_k**2)
            D += (1/3) * Q_bar * (z_k1**3 - z_k**3)

        return A, B, D

    def is_symmetric(self):
        """Determine if laminate configuration is symmetric"""
        n = len(self.layup)
        for i in range(n // 2):
            if self.layup[i] != self.layup[n - 1 - i]:
                return False
        return True

# CFRP material properties
E1 = 140.0
E2 = 10.0
nu12 = 0.30
G12 = 5.0
Q = lamina_stiffness_matrix(E1, E2, nu12, G12)

# Example laminate configurations
layup_symmetric = [0, 45, -45, 90, 90, -45, 45, 0]  # Symmetric laminate
layup_unsymmetric = [0, 45, -45, 90]                 # Unsymmetric laminate

t_ply = 0.125  # mm

# Symmetric laminate
lam_sym = Laminate(Q, layup_symmetric, t_ply)
A_sym, B_sym, D_sym = lam_sym.compute_ABD()

print("=== Symmetric Laminate [0/45/-45/90]s ===")
print(f"Laminate configuration: {layup_symmetric}")
print(f"Total thickness: {lam_sym.total_thickness} mm")
print(f"Symmetric: {lam_sym.is_symmetric()}")
print("\n[A] Matrix (N/mm):")
print(A_sym)
print("\n[B] Matrix (N·mm):")
print(B_sym)
print(f"B matrix norm: {np.linalg.norm(B_sym):.2e} (≈ 0 for symmetric laminates)")
print("\n[D] Matrix (N·mm²):")
print(D_sym)

# Unsymmetric laminate
lam_unsym = Laminate(Q, layup_unsymmetric, t_ply)
A_unsym, B_unsym, D_unsym = lam_unsym.compute_ABD()

print("\n\n=== Unsymmetric Laminate [0/45/-45/90] ===")
print(f"Laminate configuration: {layup_unsymmetric}")
print(f"Symmetric: {lam_unsym.is_symmetric()}")
print("\n[B] Matrix (N·mm):")
print(B_unsym)
print(f"B matrix norm: {np.linalg.norm(B_unsym):.2e} (extension-bending coupling present)")

2.3.3 Quasi-Isotropic Laminates

With specific laminate configurations, it is possible to design laminates with isotropic properties in the plane. The conditions for quasi-isotropy are orientations differing by \(180°/n\) with 3 or more layers:

• 3 layers: [0/60/-60] or [0/60/120]
• 4 layers: [0/45/-45/90] (most common)

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0

import numpy as np
import matplotlib.pyplot as plt

def check_quasi_isotropic(A):
    """
    Check for quasi-isotropy

    Quasi-isotropic conditions:
    - A11 = A22
    - A16 = A26 = 0
    - A66 = (A11 - A12) / 2
    """
    tol = 1e-6

    condition1 = np.abs(A[0,0] - A[1,1]) < tol
    condition2 = np.abs(A[0,2]) < tol and np.abs(A[1,2]) < tol
    condition3 = np.abs(A[2,2] - (A[0,0] - A[0,1])/2) < tol

    is_quasi_iso = condition1 and condition2 and condition3

    return is_quasi_iso, {
        'A11 = A22': condition1,
        'A16 = A26 = 0': condition2,
        'A66 = (A11-A12)/2': condition3
    }

# CFRP material
E1 = 140.0
E2 = 10.0
nu12 = 0.30
G12 = 5.0
Q = lamina_stiffness_matrix(E1, E2, nu12, G12)

# Test various laminate configurations
layups = {
    '[0/45/-45/90]': [0, 45, -45, 90],
    '[0/60/-60]': [0, 60, -60],
    '[0/90]': [0, 90],
    '[45/-45]': [45, -45]
}

t_ply = 0.125

print("Quasi-isotropy determination:\n" + "="*60)
for name, layup in layups.items():
    lam = Laminate(Q, layup, t_ply)
    A, B, D = lam.compute_ABD()

    is_qi, conditions = check_quasi_isotropic(A)

    print(f"\nLaminate configuration: {name}")
    print(f"Quasi-isotropic: {is_qi}")
    for cond, result in conditions.items():
        status = "✓" if result else "✗"
        print(f"  {status} {cond}")

    if is_qi:
        # Calculate equivalent isotropic properties
        E_eq = A[0,0] / lam.total_thickness
        nu_eq = A[0,1] / A[0,0]
        G_eq = A[2,2] / lam.total_thickness

        print(f"  Equivalent Young's modulus: {E_eq:.1f} GPa")
        print(f"  Equivalent Poisson's ratio: {nu_eq:.3f}")
        print(f"  Equivalent shear modulus: {G_eq:.1f} GPa")

2.4 Strength Prediction of Laminates

2.4.1 Maximum Stress Criterion

Stresses in the principal coordinate system of each layer must not exceed allowable values:

2.4.2 Tsai-Wu Failure Criterion

A quadratic form criterion for predicting failure under multiaxial stress states:

Determination of coefficients:

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0

import numpy as np
import matplotlib.pyplot as plt

class TsaiWuCriterion:
    """Tsai-Wu failure criterion"""

    def __init__(self, X_t, X_c, Y_t, Y_c, S):
        """
        Parameters:
        -----------
        X_t, X_c : float
            Longitudinal tensile and compressive strengths [MPa]
        Y_t, Y_c : float
            Transverse tensile and compressive strengths [MPa]
        S : float
            Shear strength [MPa]
        """
        self.X_t = X_t
        self.X_c = X_c
        self.Y_t = Y_t
        self.Y_c = Y_c
        self.S = S

        # Calculate Tsai-Wu coefficients
        self.F1 = 1/X_t - 1/X_c
        self.F2 = 1/Y_t - 1/Y_c
        self.F11 = 1/(X_t * X_c)
        self.F22 = 1/(Y_t * Y_c)
        self.F66 = 1/S**2
        self.F12 = -0.5 * np.sqrt(self.F11 * self.F22)

    def failure_index(self, sigma1, sigma2, tau12):
        """
        Calculate failure index

        Returns:
        --------
        FI : float
            Failure index (FI < 1: safe, FI = 1: failure, FI > 1: failed)
        """
        FI = (self.F1 * sigma1 + self.F2 * sigma2 +
              self.F11 * sigma1**2 + self.F22 * sigma2**2 +
              self.F66 * tau12**2 + 2 * self.F12 * sigma1 * sigma2)

        return FI

    def safety_factor(self, sigma1, sigma2, tau12):
        """Calculate safety factor"""
        FI = self.failure_index(sigma1, sigma2, tau12)
        if FI <= 0:
            return np.inf
        return 1 / np.sqrt(FI)

# CFRP/Epoxy strength properties
X_t = 1500  # MPa
X_c = 1200  # MPa
Y_t = 50    # MPa
Y_c = 200   # MPa
S = 70      # MPa

criterion = TsaiWuCriterion(X_t, X_c, Y_t, Y_c, S)

# Create failure envelope (σ1-σ2 plane)
sigma1_range = np.linspace(-X_c, X_t, 200)
sigma2_range = np.linspace(-Y_c, Y_t, 200)

# Create grid
S1, S2 = np.meshgrid(sigma1_range, sigma2_range)
FI_grid = np.zeros_like(S1)

for i in range(len(sigma2_range)):
    for j in range(len(sigma1_range)):
        FI_grid[i, j] = criterion.failure_index(S1[i, j], S2[i, j], 0)

# Visualization
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# Failure envelope (τ12 = 0)
contour = ax1.contour(S1, S2, FI_grid, levels=[1.0], colors='r', linewidths=2)
ax1.contourf(S1, S2, FI_grid, levels=[0, 1.0], colors=['lightgreen'], alpha=0.3)
ax1.set_xlabel('σ₁ [MPa]')
ax1.set_ylabel('σ₂ [MPa]')
ax1.set_title('Tsai-Wu Failure Envelope (τ₁₂ = 0)')
ax1.grid(True, alpha=0.3)
ax1.axhline(y=0, color='k', linewidth=0.5)
ax1.axvline(x=0, color='k', linewidth=0.5)

# Safety factors for specific stress states
test_cases = [
    (500, 20, 0, "Tension-Tension"),
    (500, -50, 0, "Tension-Compression"),
    (-400, -100, 0, "Compression-Compression"),
    (300, 0, 40, "Tension-Shear")
]

results = []
for sigma1, sigma2, tau12, case_name in test_cases:
    FI = criterion.failure_index(sigma1, sigma2, tau12)
    SF = criterion.safety_factor(sigma1, sigma2, tau12)
    results.append((case_name, sigma1, sigma2, tau12, FI, SF))

    if tau12 == 0:  # Plot only cases with τ12=0
        marker = 'o' if FI < 1 else 'x'
        color = 'blue' if FI < 1 else 'red'
        ax1.plot(sigma1, sigma2, marker, markersize=10, color=color,
                label=f"{case_name} (SF={SF:.2f})")

ax1.legend()

# Results table
ax2.axis('off')
table_data = [["Load State", "σ₁", "σ₂", "τ₁₂", "FI", "SF"]]
for case_name, s1, s2, t12, fi, sf in results:
    table_data.append([
        case_name,
        f"{s1:.0f}",
        f"{s2:.0f}",
        f"{t12:.0f}",
        f"{fi:.3f}",
        f"{sf:.2f}" if sf != np.inf else "∞"
    ])

table = ax2.table(cellText=table_data, cellLoc='center', loc='center',
                  colWidths=[0.25, 0.15, 0.15, 0.15, 0.15, 0.15])
table.auto_set_font_size(False)
table.set_fontsize(10)
table.scale(1, 2)

# Header row styling
for i in range(6):
    table[(0, i)].set_facecolor('#4CAF50')
    table[(0, i)].set_text_props(weight='bold', color='white')

# Data row coloring
for i in range(1, len(table_data)):
    fi_val = float(table_data[i][4])
    color = '#c8e6c9' if fi_val < 1 else '#ffcdd2'
    for j in range(6):
        table[(i, j)].set_facecolor(color)

ax2.set_title('Failure Index and Safety Factor', fontsize=12, weight='bold', pad=20)

plt.tight_layout()
plt.savefig('tsai_wu_failure.png', dpi=300, bbox_inches='tight')
plt.close()

# Console output
print("Strength evaluation using Tsai-Wu failure criterion:")
print("="*70)
for case_name, s1, s2, t12, fi, sf in results:
    status = "Safe" if fi < 1 else "Failed"
    print(f"\n{case_name}: σ1={s1} MPa, σ2={s2} MPa, τ12={t12} MPa")
    print(f"  Failure index FI = {fi:.3f} ({status})")
    print(f"  Safety factor SF = {sf:.2f}" if sf != np.inf else f"  Safety factor SF = ∞")

2.4.3 First Ply Failure (FPF)

Failure of laminates progresses in stages from the load at which the first layer fails (FPF) to final failure (LPF). In FPF analysis, stresses in each layer are calculated to identify which layer satisfies the failure criterion first.

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0, <2.0.0

import numpy as np

def first_ply_failure_analysis(lam, applied_strain, criterion):
    """
    First Ply Failure analysis

    Parameters:
    -----------
    lam : Laminate
        Laminate object
    applied_strain : ndarray (3,)
        Applied strain [εx, εy, γxy]
    criterion : TsaiWuCriterion
        Failure criterion object

    Returns:
    --------
    results : list of dict
        Stresses and failure indices for each layer
    """
    results = []

    for k in range(lam.n_plies):
        # Layer midpoint position
        z_mid = (lam.z[k] + lam.z[k+1]) / 2

        # Calculate stress in global coordinate system
        Q_bar = offaxis_stiffness(lam.Q, lam.layup[k])
        stress_global = Q_bar @ applied_strain

        # Transform to principal coordinate system
        T = transformation_matrix(lam.layup[k])
        stress_local = T @ stress_global

        sigma1, sigma2, tau12 = stress_local

        # Calculate failure index
        FI = criterion.failure_index(sigma1, sigma2, tau12)
        SF = criterion.safety_factor(sigma1, sigma2, tau12)

        results.append({
            'ply': k + 1,
            'angle': lam.layup[k],
            'z': z_mid,
            'stress_global': stress_global,
            'stress_local': stress_local,
            'FI': FI,
            'SF': SF
        })

    return results

# Laminate configuration
E1 = 140.0
E2 = 10.0
nu12 = 0.30
G12 = 5.0
Q = lamina_stiffness_matrix(E1, E2, nu12, G12)

layup = [0, 45, -45, 90]
t_ply = 0.125

lam = Laminate(Q, layup, t_ply)

# Failure criterion
X_t = 1500
X_c = 1200
Y_t = 50
Y_c = 200
S = 70
criterion = TsaiWuCriterion(X_t, X_c, Y_t, Y_c, S)

# Applied strain (x-direction tension)
applied_strain = np.array([0.005, 0, 0])

# Execute FPF analysis
results = first_ply_failure_analysis(lam, applied_strain, criterion)

# Display results
print("First Ply Failure Analysis Results:")
print("="*80)
print(f"Laminate configuration: {layup}")
print(f"Applied strain: εx = {applied_strain[0]:.4f}, εy = {applied_strain[1]:.4f}, "
      f"γxy = {applied_strain[2]:.4f}")
print("\nStress state in each layer:")
print("-"*80)
print(f"{'Ply':>3} {'Angle':>6} {'σ1':>10} {'σ2':>10} {'τ12':>10} {'FI':>8} {'SF':>8}")
print("-"*80)

fpf_ply = None
min_sf = np.inf

for r in results:
    print(f"{r['ply']:3d} {r['angle']:6.0f}° {r['stress_local'][0]:10.1f} "
          f"{r['stress_local'][1]:10.1f} {r['stress_local'][2]:10.1f} "
          f"{r['FI']:8.3f} {r['SF']:8.2f}")

    if r['SF'] < min_sf:
        min_sf = r['SF']
        fpf_ply = r['ply']

print("-"*80)
print(f"\nFirst Ply Failure: Layer {fpf_ply} (angle {layup[fpf_ply-1]}°)")
print(f"Safety factor: {min_sf:.2f}")
print(f"Failure strain: {applied_strain[0] * min_sf:.4f}")

2.5 Summary

In this chapter, we learned about the mechanics of fiber-reinforced composite materials:

• Types and manufacturing methods of CFRP/GFRP
• Stress-strain relationships and coordinate transformations for single layers
• Classical Laminate Theory (CLT) and A-B-D matrices
• Design of symmetric and quasi-isotropic laminates
• Strength prediction using Tsai-Wu failure criterion
• First Ply Failure analysis

In the next chapter, we will learn about particle-reinforced composite materials (MMC, CMC) and laminated composite materials, mastering reinforcement mechanisms and characterization methods.

Exercises

References

1. Jones, R. M., "Mechanics of Composite Materials", 2nd ed., Taylor & Francis, 1999, pp. 190-267
2. Kaw, A. K., "Mechanics of Composite Materials", 2nd ed., CRC Press, 2005, pp. 145-223
3. Tsai, S. W. and Wu, E. M., "A General Theory of Strength for Anisotropic Materials", Journal of Composite Materials, Vol. 5, 1971, pp. 58-80
4. Hyer, M. W., "Stress Analysis of Fiber-Reinforced Composite Materials", DEStech Publications, 2009, pp. 312-389
5. Barbero, E. J., "Introduction to Composite Materials Design", 2nd ed., CRC Press, 2010, pp. 89-156
6. Reddy, J. N., "Mechanics of Laminated Composite Plates and Shells: Theory and Analysis", 2nd ed., CRC Press, 2003, pp. 234-301
7. Daniel, I. M. and Ishai, O., "Engineering Mechanics of Composite Materials", 2nd ed., Oxford University Press, 2006, pp. 147-218
8. Herakovich, C. T., "Mechanics of Fibrous Composites", Wiley, 1998, pp. 178-245