Chapter 5 Python Practice

5.1 Complete Implementation of Classical Lamination Theory

5.1.1 Object-Oriented Design

We design a composite material analysis tool based on classes to ensure reusability and extensibility.

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

import numpy as np
import matplotlib.pyplot as plt
from dataclasses import dataclass
from typing import List, Tuple, Optional

@dataclass
class Material:
    """Class to hold single-layer material properties"""
    name: str
    E1: float  # Longitudinal elastic modulus [GPa]
    E2: float  # Transverse elastic modulus [GPa]
    nu12: float  # Major Poisson's ratio
    G12: float  # Shear modulus [GPa]
    Xt: float  # Longitudinal tensile strength [MPa]
    Xc: float  # Longitudinal compressive strength [MPa]
    Yt: float  # Transverse tensile strength [MPa]
    Yc: float  # Transverse compressive strength [MPa]
    S: float   # Shear strength [MPa]

    def __post_init__(self):
        """Verify reciprocal theorem"""
        self.nu21 = self.nu12 * self.E2 / self.E1

    def Q_matrix(self) -> np.ndarray:
        """Calculate reduced stiffness matrix [Q]"""
        denom = 1 - self.nu12 * self.nu21
        Q11 = self.E1 / denom
        Q22 = self.E2 / denom
        Q12 = self.nu12 * self.E2 / denom
        Q66 = self.G12

        return np.array([
            [Q11, Q12, 0],
            [Q12, Q22, 0],
            [0, 0, Q66]
        ]) * 1000  # GPa → MPa

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

    def __init__(self, material: Material, layup: List[float], t_ply: float):
        """
        Parameters:
        -----------
        material : Material
            Single-layer material
        layup : List[float]
            Stacking sequence [θ1, θ2, ..., θn] (degrees)
        t_ply : float
            Single-layer thickness [mm]
        """
        self.material = material
        self.layup = np.array(layup)
        self.t_ply = t_ply
        self.n_plies = len(layup)
        self.total_thickness = self.n_plies * t_ply

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

        # Calculate A, B, D matrices
        self.A, self.B, self.D = self._compute_ABD()

    @staticmethod
    def transformation_matrix(theta: float) -> np.ndarray:
        """Coordinate transformation matrix [T]"""
        theta_rad = np.radians(theta)
        c = np.cos(theta_rad)
        s = np.sin(theta_rad)

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

    def Q_bar(self, theta: float) -> np.ndarray:
        """Off-axis stiffness matrix [Q̄]"""
        Q = self.material.Q_matrix()
        T = self.transformation_matrix(theta)
        T_inv = np.linalg.inv(T)

        return T_inv @ Q @ T_inv.T

    def _compute_ABD(self) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
        """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_bar = self.Q_bar(self.layup[k])
            z_k = self.z[k]
            z_k1 = self.z[k + 1]

            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 ABD_matrix(self) -> np.ndarray:
        """Complete 6×6 ABD matrix"""
        return np.block([
            [self.A, self.B],
            [self.B, self.D]
        ])

    def is_symmetric(self) -> bool:
        """Determine whether stacking sequence 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

    def effective_properties(self) -> dict:
        """Calculate equivalent in-plane properties"""
        # Compliance matrix
        a = np.linalg.inv(self.A)

        # Equivalent Young's moduli
        Ex = 1 / (a[0, 0] * self.total_thickness)
        Ey = 1 / (a[1, 1] * self.total_thickness)

        # Equivalent Poisson's ratios
        nu_xy = -a[0, 1] / a[0, 0]
        nu_yx = -a[1, 0] / a[1, 1]

        # Equivalent shear modulus
        Gxy = 1 / (a[2, 2] * self.total_thickness)

        return {
            'Ex': Ex / 1000,  # MPa → GPa
            'Ey': Ey / 1000,
            'nu_xy': nu_xy,
            'nu_yx': nu_yx,
            'Gxy': Gxy / 1000
        }

    def print_summary(self):
        """Display laminate information"""
        print("="*70)
        print(f"Laminate Summary: {self.material.name}")
        print("="*70)
        print(f"Stacking Sequence: {self.layup}")
        print(f"Number of Plies: {self.n_plies}")
        print(f"Single Ply Thickness: {self.t_ply} mm")
        print(f"Total Thickness: {self.total_thickness} mm")
        print(f"Symmetric Laminate: {self.is_symmetric()}")

        print("\n[A] Matrix (N/mm):")
        print(self.A)

        print("\n[B] Matrix (N):")
        print(self.B)
        print(f"B-matrix Norm: {np.linalg.norm(self.B):.2e}")

        print("\n[D] Matrix (N·mm):")
        print(self.D)

        props = self.effective_properties()
        print("\nEquivalent In-Plane Properties:")
        print(f"  Ex = {props['Ex']:.1f} GPa")
        print(f"  Ey = {props['Ey']:.1f} GPa")
        print(f"  νxy = {props['nu_xy']:.3f}")
        print(f"  Gxy = {props['Gxy']:.1f} GPa")
        print("="*70)

# Usage example
# CFRP material definition
cfrp = Material(
    name="T300/Epoxy",
    E1=140.0, E2=10.0, nu12=0.30, G12=5.0,
    Xt=1500, Xc=1200, Yt=50, Yc=200, S=70
)

# Stacking sequences
layup_symmetric = [0, 45, -45, 90, 90, -45, 45, 0]
layup_quasi_iso = [0, 45, -45, 90]

# Create laminates
lam_sym = Laminate(cfrp, layup_symmetric, t_ply=0.125)
lam_qi = Laminate(cfrp, layup_quasi_iso, t_ply=0.125)

# Display summaries
lam_sym.print_summary()
print("\n")
lam_qi.print_summary()

5.1.2 Stress and Strain Analysis

Calculate the stress in each layer from applied loads and compare with failure criteria.

class FailureCriterion:
    """Base class for failure criteria"""

    def __init__(self, material: Material):
        self.material = material

    def failure_index(self, sigma1: float, sigma2: float, tau12: float) -> float:
        """Calculate failure index (implemented in derived classes)"""
        raise NotImplementedError

class TsaiWuCriterion(FailureCriterion):
    """Tsai-Wu failure criterion"""

    def __init__(self, material: Material):
        super().__init__(material)

        # Tsai-Wu coefficients
        self.F1 = 1/material.Xt - 1/material.Xc
        self.F2 = 1/material.Yt - 1/material.Yc
        self.F11 = 1/(material.Xt * material.Xc)
        self.F22 = 1/(material.Yt * material.Yc)
        self.F66 = 1/material.S**2
        self.F12 = -0.5 * np.sqrt(self.F11 * self.F22)

    def failure_index(self, sigma1: float, sigma2: float, tau12: float) -> float:
        """Tsai-Wu failure index"""
        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

class LaminateAnalysis:
    """Laminate load analysis class"""

    def __init__(self, laminate: Laminate, criterion: FailureCriterion):
        self.laminate = laminate
        self.criterion = criterion

    def analyze_loading(self, Nx: float, Ny: float, Nxy: float,
                        Mx: float = 0, My: float = 0, Mxy: float = 0) -> List[dict]:
        """
        Stress analysis of each layer under loading conditions

        Parameters:
        -----------
        Nx, Ny, Nxy : float
            Resultant forces [N/mm]
        Mx, My, Mxy : float
            Resultant moments [N·mm/mm]

        Returns:
        --------
        results : List[dict]
            Stress and failure index for each layer
        """
        # Inverse of ABD matrix
        ABD_inv = np.linalg.inv(self.laminate.ABD_matrix())

        # Load vector
        load = np.array([Nx, Ny, Nxy, Mx, My, Mxy])

        # Mid-plane strains and curvatures
        strain_curvature = ABD_inv @ load
        epsilon0 = strain_curvature[:3]
        kappa = strain_curvature[3:]

        results = []

        for k in range(self.laminate.n_plies):
            # Layer mid-plane position
            z_mid = (self.laminate.z[k] + self.laminate.z[k + 1]) / 2

            # Strain in global coordinate system
            epsilon_global = epsilon0 + z_mid * kappa

            # Stress in global coordinate system
            Q_bar = self.laminate.Q_bar(self.laminate.layup[k])
            stress_global = Q_bar @ epsilon_global

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

            sigma1, sigma2, tau12 = stress_local

            # Failure index
            FI = self.criterion.failure_index(sigma1, sigma2, tau12)
            SF = 1 / np.sqrt(FI) if FI > 0 else np.inf

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

        return results

    def first_ply_failure(self, Nx: float, Ny: float, Nxy: float) -> Tuple[int, float]:
        """
        Determine First Ply Failure load

        Returns:
        --------
        fpf_ply : int
            Number of first layer to fail
        fpf_load : float
            FPF load multiplier
        """
        # Analysis under unit load
        results = self.analyze_loading(Nx, Ny, Nxy)

        # Find minimum safety factor
        min_sf = min(r['SF'] for r in results)
        fpf_ply = min((r for r in results), key=lambda r: r['SF'])['ply']

        return fpf_ply, min_sf

# Usage example
cfrp = Material(
    name="T300/Epoxy",
    E1=140.0, E2=10.0, nu12=0.30, G12=5.0,
    Xt=1500, Xc=1200, Yt=50, Yc=200, S=70
)

layup = [0, 45, -45, 90]
lam = Laminate(cfrp, layup, t_ply=0.125)

# Tsai-Wu criterion
criterion = TsaiWuCriterion(cfrp)

# Analysis object
analysis = LaminateAnalysis(lam, criterion)

# Loading conditions (uniaxial tension)
Nx = 100  # N/mm
Ny = 0
Nxy = 0

# Stress analysis
results = analysis.analyze_loading(Nx, Ny, Nxy)

# Display results
print("Laminate Stress Analysis Results:")
print("="*80)
print(f"Loading: Nx = {Nx} N/mm, Ny = {Ny} N/mm, Nxy = {Nxy} N/mm")
print("-"*80)
print(f"{'Ply':>3} {'Angle':>6} {'σ1':>10} {'σ2':>10} {'τ12':>10} {'FI':>8} {'SF':>8}")
print("-"*80)

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}")

# FPF
fpf_ply, fpf_sf = analysis.first_ply_failure(Nx, Ny, Nxy)
print("-"*80)
print(f"First Ply Failure: Ply {fpf_ply} (Angle {layup[fpf_ply-1]}°)")
print(f"Safety Factor: {fpf_sf:.2f}")
print(f"Failure Load: Nx = {Nx * fpf_sf:.1f} N/mm")

5.2 Optimal Stacking Design

5.2.1 Genetic Algorithm (GA)

Genetic algorithms are effective for optimization of discrete variables (fiber orientation angles).

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

import random
from typing import List, Callable
import numpy as np

class GeneticAlgorithm:
    """Genetic algorithm for stacking design optimization"""

    def __init__(self, n_plies: int, angle_options: List[float],
                 objective_func: Callable, symmetric: bool = True,
                 pop_size: int = 50, generations: int = 100):
        """
        Parameters:
        -----------
        n_plies : int
            Number of plies
        angle_options : List[float]
            Available angles [degrees]
        objective_func : Callable
            Objective function (layup → score, smaller is better)
        symmetric : bool
            Whether to enforce symmetric laminate
        """
        self.n_plies = n_plies
        self.angle_options = angle_options
        self.objective_func = objective_func
        self.symmetric = symmetric
        self.pop_size = pop_size
        self.generations = generations

        # For symmetric laminates, treat only half as genes
        self.gene_length = n_plies // 2 if symmetric else n_plies

    def create_individual(self) -> List[float]:
        """Randomly generate individual (stacking sequence)"""
        genes = [random.choice(self.angle_options) for _ in range(self.gene_length)]

        if self.symmetric:
            # Make symmetric
            return genes + genes[::-1]
        else:
            return genes

    def fitness(self, individual: List[float]) -> float:
        """Fitness (inverse of objective function)"""
        score = self.objective_func(individual)
        return 1 / (1 + score)  # Higher fitness for smaller score

    def selection(self, population: List[List[float]]) -> List[List[float]]:
        """Tournament selection"""
        tournament_size = 3
        selected = []

        for _ in range(len(population)):
            tournament = random.sample(population, tournament_size)
            winner = max(tournament, key=self.fitness)
            selected.append(winner)

        return selected

    def crossover(self, parent1: List[float], parent2: List[float]) -> List[float]:
        """Single-point crossover"""
        point = random.randint(1, self.gene_length - 1)

        if self.symmetric:
            # Crossover on half of genes
            genes1 = parent1[:self.gene_length]
            genes2 = parent2[:self.gene_length]
            child_genes = genes1[:point] + genes2[point:]
            return child_genes + child_genes[::-1]
        else:
            return parent1[:point] + parent2[point:]

    def mutate(self, individual: List[float], mutation_rate: float = 0.1) -> List[float]:
        """Mutation"""
        if self.symmetric:
            genes = individual[:self.gene_length]
            mutated_genes = []

            for gene in genes:
                if random.random() < mutation_rate:
                    mutated_genes.append(random.choice(self.angle_options))
                else:
                    mutated_genes.append(gene)

            return mutated_genes + mutated_genes[::-1]
        else:
            return [
                random.choice(self.angle_options) if random.random() < mutation_rate else gene
                for gene in individual
            ]

    def optimize(self) -> Tuple[List[float], float]:
        """Execute optimization"""
        # Initial population
        population = [self.create_individual() for _ in range(self.pop_size)]

        best_history = []

        for gen in range(self.generations):
            # Fitness evaluation
            fitnesses = [self.fitness(ind) for ind in population]
            best_idx = np.argmax(fitnesses)
            best_individual = population[best_idx]
            best_score = self.objective_func(best_individual)

            best_history.append(best_score)

            if gen % 10 == 0:
                print(f"Generation {gen}: Best Score = {best_score:.4f}, "
                      f"Layup = {best_individual}")

            # Selection
            selected = self.selection(population)

            # Next generation
            next_population = [best_individual]  # Elitism

            while len(next_population) < self.pop_size:
                parent1, parent2 = random.sample(selected, 2)
                child = self.crossover(parent1, parent2)
                child = self.mutate(child)
                next_population.append(child)

            population = next_population

        # Best individual from final generation
        fitnesses = [self.fitness(ind) for ind in population]
        best_idx = np.argmax(fitnesses)
        best_individual = population[best_idx]
        best_score = self.objective_func(best_individual)

        return best_individual, best_score

# Define optimization problem
cfrp = Material(
    name="T300/Epoxy",
    E1=140.0, E2=10.0, nu12=0.30, G12=5.0,
    Xt=1500, Xc=1200, Yt=50, Yc=200, S=70
)

t_ply = 0.125

# Objective function: Minimize difference between Ex and Ey (approach quasi-isotropic)
def objective_quasi_isotropic(layup):
    lam = Laminate(cfrp, layup, t_ply)
    props = lam.effective_properties()
    # Relative difference between Ex and Ey
    diff = abs(props['Ex'] - props['Ey']) / props['Ex']
    return diff

# Execute GA
angle_options = [0, 45, -45, 90]
n_plies = 8

ga = GeneticAlgorithm(
    n_plies=n_plies,
    angle_options=angle_options,
    objective_func=objective_quasi_isotropic,
    symmetric=True,
    pop_size=50,
    generations=100
)

best_layup, best_score = ga.optimize()

print("\n" + "="*70)
print("Optimal Stacking Design Results:")
print("="*70)
print(f"Optimal Layup: {best_layup}")
print(f"Objective Function Value (Ex-Ey difference): {best_score:.4f}")

# Detailed analysis of optimal layup
lam_opt = Laminate(cfrp, best_layup, t_ply)
lam_opt.print_summary()

5.2.2 Multi-Objective Optimization

Simultaneously optimize multiple objectives such as strength, stiffness, and weight.

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0

from scipy.optimize import differential_evolution
import matplotlib.pyplot as plt

def multi_objective_function(layup_continuous, material, t_ply,
                              target_Nx, target_Ny):
    """
    Multi-objective function: Minimize weight & Ensure strength

    Parameters:
    -----------
    layup_continuous : array
        Stacking sequence as continuous variables [0-3] → [0, 45, -45, 90]

    Returns:
    --------
    objectives : tuple
        (weight, inverse safety factor)
    """
    # Convert continuous variables to discrete angles
    angle_map = {0: 0, 1: 45, 2: -45, 3: 90}
    layup = [angle_map[int(round(x))] for x in layup_continuous]

    # Create laminate
    lam = Laminate(material, layup, t_ply)

    # Weight (proportional to thickness)
    mass = lam.total_thickness

    # Safety factor (Tsai-Wu)
    criterion = TsaiWuCriterion(material)
    analysis = LaminateAnalysis(lam, criterion)

    try:
        results = analysis.analyze_loading(target_Nx, target_Ny, 0)
        min_sf = min(r['SF'] for r in results)
    except:
        min_sf = 0.1  # Penalty on error

    # Objectives: Minimize weight, maximize safety factor (minimize inverse)
    return mass, 1 / min_sf

# Pareto frontier exploration (simplified: scalarization method)
cfrp = Material(
    name="T300/Epoxy",
    E1=140.0, E2=10.0, nu12=0.30, G12=5.0,
    Xt=1500, Xc=1200, Yt=50, Yc=200, S=70
)

t_ply = 0.125
target_Nx = 150  # N/mm
target_Ny = 50   # N/mm
n_plies = 12

# Weighted scalarization method
weights = np.linspace(0, 1, 11)
pareto_solutions = []

for w in weights:
    def scalarized_objective(x):
        mass, inv_sf = multi_objective_function(x, cfrp, t_ply, target_Nx, target_Ny)
        # Normalize and weighted sum
        return w * mass / 2.0 + (1 - w) * inv_sf * 10

    # Optimization (differential_evolution)
    bounds = [(0, 3)] * n_plies  # Continuous values 0-3
    result = differential_evolution(
        scalarized_objective,
        bounds,
        maxiter=50,
        seed=123,
        atol=0.1,
        tol=0.1
    )

    # Optimal solution
    angle_map = {0: 0, 1: 45, 2: -45, 3: 90}
    best_layup = [angle_map[int(round(x))] for x in result.x]
    mass, inv_sf = multi_objective_function(result.x, cfrp, t_ply, target_Nx, target_Ny)

    pareto_solutions.append({
        'weight': w,
        'layup': best_layup,
        'mass': mass,
        'safety_factor': 1 / inv_sf
    })

    print(f"Weight w={w:.1f}: Mass={mass:.3f} mm, SF={1/inv_sf:.2f}, Layup={best_layup}")

# Visualize Pareto front
masses = [sol['mass'] for sol in pareto_solutions]
sfs = [sol['safety_factor'] for sol in pareto_solutions]

plt.figure(figsize=(10, 6))
plt.plot(masses, sfs, 'bo-', linewidth=2, markersize=8)
plt.xlabel('Thickness [mm]')
plt.ylabel('Safety Factor')
plt.title('Pareto Frontier (Weight vs Safety Factor)')
plt.grid(True, alpha=0.3)

# Label each point
for i, sol in enumerate(pareto_solutions[::2]):  # Display every other
    plt.annotate(f"w={sol['weight']:.1f}",
                 (sol['mass'], sol['safety_factor']),
                 textcoords="offset points", xytext=(5,5), fontsize=8)

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

5.3 Finite Element Method Preprocessing

5.3.1 Mesh Generation

In finite element analysis of composite materials, each layer is treated as separate elements or as integration points of shell elements.

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

import numpy as np

class CompositeMesh:
    """Mesh generation for composite material plates"""

    def __init__(self, length: float, width: float,
                 nx: int, ny: int, laminate: Laminate):
        """
        Parameters:
        -----------
        length, width : float
            Plate dimensions [mm]
        nx, ny : int
            Element divisions
        laminate : Laminate
            Laminate object
        """
        self.length = length
        self.width = width
        self.nx = nx
        self.ny = ny
        self.laminate = laminate

        self.nodes = []
        self.elements = []

        self._generate_mesh()

    def _generate_mesh(self):
        """Generate nodes and elements"""
        # Node generation
        dx = self.length / self.nx
        dy = self.width / self.ny

        node_id = 1
        for j in range(self.ny + 1):
            for i in range(self.nx + 1):
                x = i * dx
                y = j * dy
                self.nodes.append({
                    'id': node_id,
                    'x': x,
                    'y': y,
                    'z': 0
                })
                node_id += 1

        # Element generation (4-node shell elements)
        elem_id = 1
        for j in range(self.ny):
            for i in range(self.nx):
                n1 = j * (self.nx + 1) + i + 1
                n2 = n1 + 1
                n3 = n1 + (self.nx + 1) + 1
                n4 = n1 + (self.nx + 1)

                self.elements.append({
                    'id': elem_id,
                    'nodes': [n1, n2, n3, n4]
                })
                elem_id += 1

    def export_abaqus_inp(self, filename: str):
        """Export Abaqus input file"""
        with open(filename, 'w') as f:
            f.write("*HEADING\n")
            f.write("Composite Laminate Mesh\n")

            # Nodes
            f.write("*NODE\n")
            for node in self.nodes:
                f.write(f"{node['id']}, {node['x']:.4f}, {node['y']:.4f}, {node['z']:.4f}\n")

            # Elements
            f.write("*ELEMENT, TYPE=S4R, ELSET=PLATE\n")
            for elem in self.elements:
                nodes_str = ", ".join(map(str, elem['nodes']))
                f.write(f"{elem['id']}, {nodes_str}\n")

            # Shell section
            f.write("*SHELL SECTION, ELSET=PLATE, COMPOSITE\n")
            for k, angle in enumerate(self.laminate.layup):
                # Thickness, integration points, material name, angle
                f.write(f"{self.laminate.t_ply}, 3, MAT1, {angle}\n")

            # Material properties
            mat = self.laminate.material
            f.write("*MATERIAL, NAME=MAT1\n")
            f.write("*ELASTIC, TYPE=LAMINA\n")
            f.write(f"{mat.E1*1000}, {mat.E2*1000}, {mat.nu12}, "
                    f"{mat.G12*1000}, {mat.G12*1000}, {mat.G12*1000}\n")

            # Boundary conditions (simply supported)
            f.write("*BOUNDARY\n")
            # Left edge (x=0): UX=0
            for node in self.nodes:
                if abs(node['x']) < 1e-6:
                    f.write(f"{node['id']}, 1\n")

            # Bottom edge (y=0): UY=0
            for node in self.nodes:
                if abs(node['y']) < 1e-6:
                    f.write(f"{node['id']}, 2\n")

            # Load step
            f.write("*STEP\n")
            f.write("*STATIC\n")

            # Distributed load (pressure on top surface)
            f.write("*DLOAD\n")
            for elem in self.elements:
                f.write(f"{elem['id']}, P, 0.1\n")  # 0.1 MPa

            f.write("*OUTPUT, FIELD\n")
            f.write("*NODE OUTPUT\n")
            f.write("U, RF\n")
            f.write("*ELEMENT OUTPUT\n")
            f.write("S, E\n")
            f.write("*END STEP\n")

        print(f"Abaqus input file exported: {filename}")

# Mesh generation example
cfrp = Material(
    name="T300/Epoxy",
    E1=140.0, E2=10.0, nu12=0.30, G12=5.0,
    Xt=1500, Xc=1200, Yt=50, Yc=200, S=70
)

layup = [0, 45, -45, 90, 90, -45, 45, 0]
lam = Laminate(cfrp, layup, t_ply=0.125)

# Rectangular plate mesh
mesh = CompositeMesh(length=100, width=100, nx=10, ny=10, laminate=lam)

# Export Abaqus input file
mesh.export_abaqus_inp("composite_plate.inp")

print(f"Generated Mesh Information:")
print(f"  Number of Nodes: {len(mesh.nodes)}")
print(f"  Number of Elements: {len(mesh.elements)}")
print(f"  Stacking Sequence: {layup}")
print(f"  Total Thickness: {lam.total_thickness} mm")

5.3.2 Post-Processing and Data Visualization

Automate reading and visualization of FEA results.

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

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
from matplotlib.collections import PatchCollection

def visualize_stress_distribution(mesh: CompositeMesh, stress_values: np.ndarray,
                                    component: str = 'Sxx', cmap: str = 'jet'):
    """
    Visualize stress distribution on mesh

    Parameters:
    -----------
    mesh : CompositeMesh
        Mesh object
    stress_values : ndarray
        Stress value for each element [MPa]
    component : str
        Stress component name
    cmap : str
        Colormap
    """
    fig, ax = plt.subplots(figsize=(10, 8))

    patches = []
    colors = []

    for elem, stress in zip(mesh.elements, stress_values):
        # Get coordinates of 4 nodes of element
        node_ids = elem['nodes']
        coords = np.array([[mesh.nodes[nid-1]['x'], mesh.nodes[nid-1]['y']]
                           for nid in node_ids])

        # Create rectangular patch
        x_min, y_min = coords.min(axis=0)
        width = coords[:, 0].max() - x_min
        height = coords[:, 1].max() - y_min

        rect = Rectangle((x_min, y_min), width, height)
        patches.append(rect)
        colors.append(stress)

    # Patch collection
    p = PatchCollection(patches, cmap=cmap, edgecolors='black', linewidths=0.5)
    p.set_array(np.array(colors))

    ax.add_collection(p)

    # Colorbar
    cbar = plt.colorbar(p, ax=ax)
    cbar.set_label(f'{component} [MPa]', fontsize=12)

    ax.set_xlim(0, mesh.length)
    ax.set_ylim(0, mesh.width)
    ax.set_aspect('equal')
    ax.set_xlabel('X [mm]', fontsize=12)
    ax.set_ylabel('Y [mm]', fontsize=12)
    ax.set_title(f'Stress Distribution: {component}', fontsize=14, weight='bold')

    plt.tight_layout()
    plt.savefig(f'stress_{component}.png', dpi=300, bbox_inches='tight')
    plt.close()

# Generate simulated stress data
n_elements = len(mesh.elements)

# Simulate distribution with higher stress at plate center
stress_Sxx = []
for elem in mesh.elements:
    node_ids = elem['nodes']
    x_center = np.mean([mesh.nodes[nid-1]['x'] for nid in node_ids])
    y_center = np.mean([mesh.nodes[nid-1]['y'] for nid in node_ids])

    # Stress based on distance from plate center (50, 50)
    dist = np.sqrt((x_center - 50)**2 + (y_center - 50)**2)
    stress = 100 * np.exp(-dist / 30)  # Gaussian distribution

    stress_Sxx.append(stress)

stress_Sxx = np.array(stress_Sxx)

# Visualization
visualize_stress_distribution(mesh, stress_Sxx, component='Sxx', cmap='jet')
print("Stress distribution diagram exported: stress_Sxx.png")

5.4 Summary

In this chapter, we learned practical implementation of composite material analysis using Python:

• Complete implementation of Classical Lamination Theory (CLT) with object-oriented design
• Stress analysis and First Ply Failure prediction
• Optimal stacking design using genetic algorithms
• Multi-objective optimization and Pareto frontier
• Finite element method preprocessing (mesh generation, Abaqus input)
• Result visualization and post-processing

By combining these techniques, professional-level composite material design and analysis becomes possible. For advanced study, we recommend engaging with damage mechanics, probabilistic design, and multiscale analysis.

Exercises

References

1. Reddy, J. N., "Mechanics of Laminated Composite Plates and Shells: Theory and Analysis", 2nd ed., CRC Press, 2003, pp. 456-534
2. Kaw, A. K., "Mechanics of Composite Materials", 2nd ed., CRC Press, 2005, pp. 312-389
3. Goldberg, D. E., "Genetic Algorithms in Search, Optimization, and Machine Learning", Addison-Wesley, 1989, pp. 1-89
4. Deb, K., "Multi-Objective Optimization Using Evolutionary Algorithms", Wiley, 2001, pp. 234-312
5. Liu, B., Haftka, R. T., and Akgun, M. A., "Two-level Composite Wing Structural Optimization Using Response Surfaces", Structural and Multidisciplinary Optimization, Vol. 20, 2000, pp. 87-96
6. Simulia, "Abaqus Analysis User's Guide: Composite Materials", Dassault Systèmes, 2020, pp. 23.1.1-23.6.8
7. Bathe, K. J., "Finite Element Procedures", Prentice Hall, 1996, pp. 634-712
8. Hunter, J. D., "Matplotlib: A 2D Graphics Environment", Computing in Science & Engineering, Vol. 9, 2007, pp. 90-95