"""
分子动力学模拟 - 简单Lennard-Jones流体
演示原子尺度的多体相互作用模拟
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


class MolecularDynamics:
    """分子动力学模拟器"""
    
    def __init__(self, n_atoms, box_size, temperature, dt=0.001):
        """
        初始化MD系统
        
        参数:
            n_atoms: 原子数量
            box_size: 模拟盒子尺寸 (Lx, Ly, Lz)
            temperature: 初始温度
            dt: 时间步长
        """
        self.n_atoms = n_atoms
        self.box_size = np.array(box_size)
        self.dt = dt
        
        # 初始化位置（简单立方晶格）
        self.positions = self._init_lattice()
        
        # 初始化速度（Maxwell-Boltzmann分布）
        self.velocities = self._init_velocities(temperature)
        
        # 力
        self.forces = np.zeros_like(self.positions)
        
        # 能量记录
        self.energies = []
        self.temperatures = []
        
        # Lennard-Jones参数（无量纲化）
        self.epsilon = 1.0  # 势阱深度
        self.sigma = 1.0    # 原子直径
        self.cutoff = 2.5 * self.sigma  # 截断半径
    
    def _init_lattice(self):
        """初始化简单立方晶格"""
        # 计算晶格点数
        n_per_side = int(np.ceil(self.n_atoms**(1/3)))
        spacing = self.box_size[0] / n_per_side
        
        positions = []
        for i in range(n_per_side):
            for j in range(n_per_side):
                for k in range(n_per_side):
                    if len(positions) < self.n_atoms:
                        pos = np.array([i, j, k]) * spacing + spacing/2
                        positions.append(pos)
        
        return np.array(positions[:self.n_atoms])
    
    def _init_velocities(self, temperature):
        """初始化速度（Maxwell-Boltzmann分布）"""
        # 从正态分布采样
        velocities = np.random.randn(self.n_atoms, 3)
        
        # 移除质心速度
        velocities -= np.mean(velocities, axis=0)
        
        # 调整温度
        current_temp = np.sum(velocities**2) / (3 * self.n_atoms)
        velocities *= np.sqrt(temperature / current_temp)
        
        return velocities
    
    def _apply_pbc(self, positions):
        """应用周期性边界条件"""
        return positions - self.box_size * np.floor(positions / self.box_size)
    
    def _minimum_image(self, r):
        """最小像约定"""
        return r - self.box_size * np.round(r / self.box_size)
    
    def compute_forces(self):
        """计算Lennard-Jones力"""
        self.forces.fill(0)
        potential_energy = 0.0
        
        for i in range(self.n_atoms):
            for j in range(i+1, self.n_atoms):
                # 最小像距离
                r_ij = self._minimum_image(self.positions[i] - self.positions[j])
                r = np.linalg.norm(r_ij)
                
                if r < self.cutoff and r > 0:
                    # LJ势
                    sr6 = (self.sigma / r)**6
                    sr12 = sr6**2
                    
                    # 势能
                    potential_energy += 4 * self.epsilon * (sr12 - sr6)
                    
                    # 力
                    f_mag = 24 * self.epsilon * (2*sr12 - sr6) / r
                    f_vec = f_mag * r_ij / r
                    
                    self.forces[i] += f_vec
                    self.forces[j] -= f_vec
        
        return potential_energy
    
    def velocity_verlet(self):
        """速度Verlet积分"""
        # 半步速度更新
        self.velocities += 0.5 * self.forces * self.dt
        
        # 位置更新
        self.positions += self.velocities * self.dt
        self.positions = self._apply_pbc(self.positions)
        
        # 计算新力
        pe = self.compute_forces()
        
        # 另半步速度更新
        self.velocities += 0.5 * self.forces * self.dt
        
        # 动能
        ke = 0.5 * np.sum(self.velocities**2)
        
        # 温度
        temp = 2 * ke / (3 * self.n_atoms)
        
        return pe, ke, temp
    
    def run(self, n_steps, thermo_interval=100):
        """
        运行MD模拟
        
        参数:
            n_steps: 总步数
            thermo_interval: 热力学量输出间隔
        """
        print(f"运行MD模拟: {self.n_atoms} 原子, {n_steps} 步")
        
        # 初始力计算
        self.compute_forces()
        
        for step in range(n_steps):
            pe, ke, temp = self.velocity_verlet()
            
            if step % thermo_interval == 0:
                total_energy = pe + ke
                self.energies.append([step, pe, ke, total_energy])
                self.temperatures.append([step, temp])
                
                if step % 1000 == 0:
                    print(f"  Step {step:6d}: PE={pe:10.2f}, KE={ke:10.2f}, "
                          f"Total={total_energy:10.2f}, T={temp:.3f}")
        
        print("✓ MD模拟完成")
        return np.array(self.energies), np.array(self.temperatures)
    
    def get_radial_distribution(self, n_bins=100):
        """计算径向分布函数 g(r)"""
        dr = self.box_size[0] / (2 * n_bins)
        r_bins = np.linspace(0, self.box_size[0]/2, n_bins)
        g_r = np.zeros(n_bins)
        
        for i in range(self.n_atoms):
            for j in range(i+1, self.n_atoms):
                r_ij = self._minimum_image(self.positions[i] - self.positions[j])
                r = np.linalg.norm(r_ij)
                
                if r < self.box_size[0]/2:
                    bin_idx = int(r / dr)
                    if bin_idx < n_bins:
                        g_r[bin_idx] += 2  # 每对原子贡献2
        
        # 归一化
        rho = self.n_atoms / np.prod(self.box_size)
        for i in range(n_bins):
            r = r_bins[i]
            shell_volume = 4 * np.pi * r**2 * dr
            if shell_volume > 0:
                g_r[i] /= (self.n_atoms * rho * shell_volume)
        
        return r_bins, g_r


def visualize_md_results(md_sim, energies, temperatures):
    """可视化MD模拟结果"""
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：能量演化
    ax1 = axes[0, 0]
    ax1.plot(energies[:, 0], energies[:, 1], 'b-', label='势能', linewidth=1)
    ax1.plot(energies[:, 0], energies[:, 2], 'r-', label='动能', linewidth=1)
    ax1.plot(energies[:, 0], energies[:, 3], 'g-', label='总能量', linewidth=1.5)
    ax1.set_xlabel('步数')
    ax1.set_ylabel('能量')
    ax1.set_title('能量演化')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 图2：温度演化
    ax2 = axes[0, 1]
    ax2.plot(temperatures[:, 0], temperatures[:, 1], 'm-', linewidth=1)
    ax2.set_xlabel('步数')
    ax2.set_ylabel('温度')
    ax2.set_title('温度演化')
    ax2.grid(True, alpha=0.3)
    
    # 图3：原子位置分布
    ax3 = axes[1, 0]
    ax3.scatter(md_sim.positions[:, 0], md_sim.positions[:, 1], 
               c=md_sim.positions[:, 2], cmap='viridis', s=20)
    ax3.set_xlabel('X')
    ax3.set_ylabel('Y')
    ax3.set_title('原子位置分布 (XY平面)')
    ax3.set_aspect('equal')
    
    # 图4：径向分布函数
    ax4 = axes[1, 1]
    r, g_r = md_sim.get_radial_distribution()
    ax4.plot(r[1:], g_r[1:], 'b-', linewidth=1.5)
    ax4.set_xlabel('r')
    ax4.set_ylabel('g(r)')
    ax4.set_title('径向分布函数')
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('molecular_dynamics_results.png', dpi=150, bbox_inches='tight')
    print("✓ MD结果已保存: molecular_dynamics_results.png")
    plt.close()


# 主程序
if __name__ == "__main__":
    print("=" * 60)
    print("分子动力学模拟")
    print("=" * 60)
    
    # 系统参数
    n_atoms = 256
    box_size = (10.0, 10.0, 10.0)
    temperature = 1.0
    dt = 0.001
    n_steps = 5000
    
    print(f"\n系统参数:")
    print(f"  原子数: {n_atoms}")
    print(f"  盒子尺寸: {box_size}")
    print(f"  初始温度: {temperature}")
    print(f"  时间步长: {dt}")
    print(f"  模拟步数: {n_steps}")
    
    # 创建MD系统
    md = MolecularDynamics(n_atoms, box_size, temperature, dt)
    
    # 运行模拟
    print("\n开始MD模拟...")
    energies, temperatures = md.run(n_steps)
    
    # 输出统计
    print("\n模拟统计:")
    print(f"  平均温度: {np.mean(temperatures[:, 1]):.3f}")
    print(f"  最终总能量: {energies[-1, 3]:.2f}")
    print(f"  能量漂移: {energies[-1, 3] - energies[0, 3]:.4f}")
    
    # 可视化
    print("\n生成可视化结果...")
    visualize_md_results(md, energies, temperatures)
    
    print("\n" + "=" * 60)
    print("✓ 分子动力学模拟完成!")
    print("=" * 60)

7.2 案例二：代表性体积单元(RVE)分析

本案例实现复合材料RVE的均匀化分析，演示微观到宏观的尺度跨越。

"""
代表性体积单元(RVE)分析
复合材料微观结构均匀化
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


class RVEAnalysis:
    """RVE均匀化分析"""
    
    def __init__(self, size, n_elements, fiber_volume_fraction=0.3):
        """
        初始化RVE
        
        参数:
            size: RVE尺寸 (Lx, Ly)
            n_elements: 网格划分 (nx, ny)
            fiber_volume_fraction: 纤维体积分数
        """
        self.Lx, self.Ly = size
        self.nx, self.ny = n_elements
        self.dx = self.Lx / self.nx
        self.dy = self.Ly / self.ny
        
        self.fiber_vf = fiber_volume_fraction
        
        # 材料参数
        self.E_matrix = 3.5e9  # 基体弹性模量 (Pa)
        self.nu_matrix = 0.35   # 基体泊松比
        
        self.E_fiber = 230e9   # 纤维弹性模量 (Pa)
        self.nu_fiber = 0.2     # 纤维泊松比
        
        # 生成微观结构
        self.phase = self._generate_microstructure()
        
        # 计算材料属性场
        self.E_field, self.nu_field = self._compute_material_fields()
    
    def _generate_microstructure(self):
        """生成随机纤维分布"""
        phase = np.zeros((self.ny, self.nx), dtype=int)
        
        # 纤维半径
        fiber_radius = np.sqrt(self.fiber_vf * self.Lx * self.Ly / (np.pi * 20))
        
        # 随机放置纤维
        n_fibers = 20
        np.random.seed(42)
        
        for _ in range(n_fibers):
            cx = np.random.uniform(0, self.Lx)
            cy = np.random.uniform(0, self.Ly)
            
            for i in range(self.ny):
                for j in range(self.nx):
                    x = j * self.dx
                    y = i * self.dy
                    
                    # 周期性边界
                    dx = min(abs(x - cx), self.Lx - abs(x - cx))
                    dy = min(abs(y - cy), self.Ly - abs(y - cy))
                    r = np.sqrt(dx**2 + dy**2)
                    
                    if r < fiber_radius:
                        phase[i, j] = 1  # 纤维
        
        return phase
    
    def _compute_material_fields(self):
        """计算材料属性场"""
        E_field = np.where(self.phase == 1, self.E_fiber, self.E_matrix)
        nu_field = np.where(self.phase == 1, self.nu_fiber, self.nu_matrix)
        return E_field, nu_field
    
    def _compute_stiffness_matrix(self, E, nu):
        """计算平面应力刚度矩阵"""
        # 平面应力本构矩阵
        factor = E / (1 - nu**2)
        D = factor * np.array([
            [1, nu, 0],
            [nu, 1, 0],
            [0, 0, (1-nu)/2]
        ])
        return D
    
    def _get_element_stiffness(self, i, j):
        """获取单元刚度矩阵"""
        E = self.E_field[i, j]
        nu = self.nu_field[i, j]
        D = self._compute_stiffness_matrix(E, nu)
        
        # 简化的单元刚度（假设为双线性四边形）
        # 这里使用简化的刚度矩阵
        k = E / ((1 + nu) * (1 - 2*nu))
        Ke = k * np.array([
            [ 2, -1, -1,  0],
            [-1,  2,  0, -1],
            [-1,  0,  2, -1],
            [ 0, -1, -1,  2]
        ]) * self.dx * self.dy / 4
        
        return Ke
    
    def apply_periodic_bc(self, macro_strain):
        """
        应用周期性边界条件和宏观应变
        
        参数:
            macro_strain: 宏观应变 [εxx, εyy, γxy]
        """
        # 简化的均匀化计算
        # 实际实现需要完整的有限元组装和求解
        
        # 基于混合律估算
        E_eff = self.fiber_vf * self.E_fiber + (1 - self.fiber_vf) * self.E_matrix
        nu_eff = self.fiber_vf * self.nu_fiber + (1 - self.fiber_vf) * self.nu_matrix
        
        # 计算宏观应力
        D_eff = self._compute_stiffness_matrix(E_eff, nu_eff)
        macro_stress = D_eff @ macro_strain
        
        return E_eff, nu_eff, macro_stress
    
    def compute_effective_properties(self):
        """计算等效弹性性质"""
        print("计算RVE等效性质...")
        
        # 施加不同方向的应变
        strain_cases = [
            np.array([0.01, 0, 0]),      # 单轴X
            np.array([0, 0.01, 0]),      # 单轴Y
            np.array([0, 0, 0.01]),      # 剪切
        ]
        
        results = []
        for strain in strain_cases:
            E_eff, nu_eff, stress = self.apply_periodic_bc(strain)
            results.append({
                'strain': strain,
                'stress': stress,
                'E_eff': E_eff,
                'nu_eff': nu_eff
            })
        
        # 提取等效性质
        E_x = results[0]['stress'][0] / results[0]['strain'][0]
        E_y = results[1]['stress'][1] / results[1]['strain'][1]
        G_xy = results[2]['stress'][2] / results[2]['strain'][2]
        nu_xy = -results[0]['stress'][1] / results[0]['stress'][0]
        
        print(f"✓ 等效性质计算完成")
        
        return {
            'E_x': E_x,
            'E_y': E_y,
            'G_xy': G_xy,
            'nu_xy': nu_xy,
            'details': results
        }
    
    def compute_stress_field(self, macro_strain):
        """计算微观应力场"""
        # 简化的应力场计算
        stress_field = np.zeros((self.ny, self.nx, 3))
        
        for i in range(self.ny):
            for j in range(self.nx):
                E = self.E_field[i, j]
                nu = self.nu_field[i, j]
                D = self._compute_stiffness_matrix(E, nu)
                stress_field[i, j] = D @ macro_strain
        
        return stress_field


def visualize_rve(rve, effective_props):
    """可视化RVE分析结果"""
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 图1：微观结构
    ax1 = axes[0, 0]
    im1 = ax1.imshow(rve.phase, extent=[0, rve.Lx, 0, rve.Ly], 
                     cmap='Blues', origin='lower')
    ax1.set_xlabel('X (μm)')
    ax1.set_ylabel('Y (μm)')
    ax1.set_title(f'RVE微观结构 (纤维体积分数={rve.fiber_vf:.2f})')
    
    # 图2：弹性模量分布
    ax2 = axes[0, 1]
    im2 = ax2.imshow(rve.E_field/1e9, extent=[0, rve.Lx, 0, rve.Ly],
                     cmap='viridis', origin='lower')
    ax2.set_xlabel('X (μm)')
    ax2.set_ylabel('Y (μm)')
    ax2.set_title('弹性模量分布 (GPa)')
    plt.colorbar(im2, ax=ax2)
    
    # 图3：应力场
    ax3 = axes[1, 0]
    macro_strain = np.array([0.01, 0, 0])
    stress_field = rve.compute_stress_field(macro_strain)
    von_mises = np.sqrt(stress_field[:,:,0]**2 - stress_field[:,:,0]*stress_field[:,:,1] + 
                        stress_field[:,:,1]**2 + 3*stress_field[:,:,2]**2)
    im3 = ax3.imshow(von_mises/1e6, extent=[0, rve.Lx, 0, rve.Ly],
                     cmap='hot', origin='lower')
    ax3.set_xlabel('X (μm)')
    ax3.set_ylabel('Y (μm)')
    ax3.set_title('Von Mises应力 (MPa)')
    plt.colorbar(im3, ax=ax3)
    
    # 图4：等效性质对比
    ax4 = axes[1, 1]
    
    props = ['E_x', 'E_y', 'G_xy']
    values = [effective_props[p]/1e9 for p in props]
    
    # 理论预测（混合律）
    E_rule_of_mixtures = (rve.fiber_vf * rve.E_fiber + 
                          (1-rve.fiber_vf) * rve.E_matrix) / 1e9
    
    x_pos = np.arange(len(props))
    ax4.bar(x_pos, values, alpha=0.7, label='RVE计算')
    ax4.axhline(y=E_rule_of_mixtures, color='r', linestyle='--', 
                label=f'混合律预测 ({E_rule_of_mixtures:.1f} GPa)')
    
    ax4.set_xticks(x_pos)
    ax4.set_xticklabels(props)
    ax4.set_ylabel('模量 (GPa)')
    ax4.set_title('等效弹性性质')
    ax4.legend()
    ax4.grid(True, alpha=0.3, axis='y')
    
    plt.tight_layout()
    plt.savefig('rve_analysis_results.png', dpi=150, bbox_inches='tight')
    print("✓ RVE分析结果已保存: rve_analysis_results.png")
    plt.close()


# 主程序
if __name__ == "__main__":
    print("=" * 60)
    print("代表性体积单元(RVE)分析")
    print("=" * 60)
    
    # RVE参数
    size = (100e-6, 100e-6)  # 100 μm x 100 μm
    n_elements = (50, 50)
    fiber_vf = 0.3
    
    print(f"\nRVE参数:")
    print(f"  尺寸: {size[0]*1e6:.0f} x {size[1]*1e6:.0f} μm")
    print(f"  网格: {n_elements[0]} x {n_elements[1]}")
    print(f"  纤维体积分数: {fiber_vf:.2f}")
    
    print(f"\n材料参数:")
    print(f"  基体: E={3.5:.1f} GPa, ν={0.35:.2f}")
    print(f"  纤维: E={230:.0f} GPa, ν={0.20:.2f}")
    
    # 创建RVE
    rve = RVEAnalysis(size, n_elements, fiber_vf)
    
    # 计算等效性质
    print("\n开始RVE分析...")
    effective_props = rve.compute_effective_properties()
    
    # 输出结果
    print("\n等效弹性性质:")
    print(f"  E_x = {effective_props['E_x']/1e9:.2f} GPa")
    print(f"  E_y = {effective_props['E_y']/1e9:.2f} GPa")
    print(f"  G_xy = {effective_props['G_xy']/1e9:.2f} GPa")
    print(f"  ν_xy = {effective_props['nu_xy']:.4f}")
    
    # 可视化
    print("\n生成可视化结果...")
    visualize_rve(rve, effective_props)
    
    print("\n" + "=" * 60)
    print("✓ RVE分析完成!")
    print("=" * 60)

7.3 案例三：多尺度有限元计算

本案例实现简化的FE²多尺度计算框架，演示并发多尺度耦合。

"""
多尺度有限元计算 - 简化的FE²框架
演示宏观-微观耦合计算
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import csr_matrix, lil_matrix
from scipy.sparse.linalg import spsolve
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


class MicroModel:
    """微观模型（RVE）"""
    
    def __init__(self, n_points=10):
        """
        初始化微观模型
        
        参数:
            n_points: 微观网格点数
        """
        self.n_points = n_points
        self.xi = np.linspace(-0.5, 0.5, n_points)
        self.yi = np.linspace(-0.5, 0.5, n_points)
        self.dxi = 1.0 / (n_points - 1)
        
        # 微观材料属性（非均匀）
        self.E_micro = self._generate_micro_structure()
    
    def _generate_micro_structure(self):
        """生成微观结构"""
        E = np.ones((self.n_points, self.n_points))
        
        # 添加周期性非均匀性
        for i in range(self.n_points):
            for j in range(self.n_points):
                x, y = self.xi[j], self.yi[i]
                # 周期性变化
                E[i, j] = 1.0 + 0.3 * np.sin(2*np.pi*x) * np.sin(2*np.pi*y)
        
        return E
    
    def solve(self, macro_strain):
        """
        求解微观边界值问题
        
        参数:
            macro_strain: 宏观应变 [εxx, εyy, γxy]
        
        返回:
            homogenized_stress: 均匀化应力
            effective_stiffness: 等效刚度矩阵
        """
        # 简化的均匀化计算
        # 实际应求解周期性边界条件下的微观问题
        
        E_avg = np.mean(self.E_micro)
        nu = 0.3
        
        # 平面应力本构
        factor = E_avg / (1 - nu**2)
        D_eff = factor * np.array([
            [1, nu, 0],
            [nu, 1, 0],
            [0, 0, (1-nu)/2]
        ])
        
        stress = D_eff @ macro_strain
        
        return stress, D_eff


class MacroModel:
    """宏观有限元模型"""
    
    def __init__(self, length, n_elements, micro_models):
        """
        初始化宏观模型
        
        参数:
            length: 结构长度
            n_elements: 单元数
            micro_models: 每个积分点的微观模型列表
        """
        self.L = length
        self.n_elem = n_elements
        self.n_nodes = n_elements + 1
        self.h = length / n_elements
        
        self.micro_models = micro_models
        
        # 节点坐标
        self.x_nodes = np.linspace(0, length, self.n_nodes)
    
    def _get_micro_response(self, elem_id, strain):
        """获取微观响应"""
        micro = self.micro_models[elem_id]
        stress, D_eff = micro.solve(strain)
        return stress, D_eff
    
    def solve_nonlinear(self, load_steps=10):
        """
        非线性求解（简化的一维问题）
        
        参数:
            load_steps: 加载步数
        
        返回:
            displacements: 节点位移
            strains: 单元应变
            stresses: 单元应力
        """
        print("求解宏观问题...")
        
        # 初始化
        u = np.zeros(self.n_nodes)
        strains = np.zeros(self.n_elem)
        stresses = np.zeros(self.n_elem)
        
        # 逐步加载
        max_load = 100e6  # 最大载荷 (Pa)
        
        for step in range(1, load_steps + 1):
            load_factor = step / load_steps
            applied_load = max_load * load_factor
            
            # 组装刚度矩阵（简化的一维）
            K = lil_matrix((self.n_nodes, self.n_nodes))
            F = np.zeros(self.n_nodes)
            
            for e in range(self.n_elem):
                # 获取微观响应
                strain_e = np.array([strains[e], 0, 0])  # 假设单轴
                stress_e, D_eff = self._get_micro_response(e, strain_e)
                
                # 单元刚度
                E_eff = D_eff[0, 0]
                k_e = E_eff * 1.0 / self.h  # 假设单位面积
                
                # 组装
                K[e, e] += k_e
                K[e, e+1] -= k_e
                K[e+1, e] -= k_e
                K[e+1, e+1] += k_e
                
                # 内力
                F[e] -= stress_e[0] * 1.0
                F[e+1] += stress_e[0] * 1.0
            
            # 应用边界条件
            K[0, :] = 0
            K[0, 0] = 1
            F[0] = 0  # 固定端
            
            F[-1] += applied_load * 1.0  # 载荷端
            
            # 求解
            K = K.tocsr()
            du = spsolve(K, F)
            u += du
            
            # 更新应变和应力
            for e in range(self.n_elem):
                strains[e] = (u[e+1] - u[e]) / self.h
                strain_e = np.array([strains[e], 0, 0])
                stress_e, _ = self._get_micro_response(e, strain_e)
                stresses[e] = stress_e[0]
            
            if step % 2 == 0:
                print(f"  Step {step}/{load_steps}: Load={applied_load/1e6:.1f} MPa, "
                      f"Max disp={np.max(np.abs(u))*1e3:.3f} mm")
        
        print("✓ 宏观求解完成")
        return u, strains, stresses


class MultiscaleFE2:
    """FE²多尺度框架"""
    
    def __init__(self, macro_length, n_macro_elem, n_micro_points):
        """
        初始化FE²框架
        
        参数:
            macro_length: 宏观结构长度
            n_macro_elem: 宏观单元数
            n_micro_points: 微观模型网格点数
        """
        self.macro_length = macro_length
        self.n_macro_elem = n_macro_elem
        self.n_micro_points = n_micro_points
        
        # 为每个宏观积分点创建微观模型
        print(f"初始化FE²框架: {n_macro_elem} 宏观单元, 每个含微观模型")
        self.micro_models = [MicroModel(n_micro_points) for _ in range(n_macro_elem)]
        
        # 创建宏观模型
        self.macro_model = MacroModel(macro_length, n_macro_elem, self.micro_models)
    
    def solve(self, load_steps=10):
        """
        求解多尺度问题
        
        参数:
            load_steps: 加载步数
        
        返回:
            results: 结果字典
        """
        print("\n开始FE²多尺度计算...")
        
        # 求解宏观问题
        u, strains, stresses = self.macro_model.solve_nonlinear(load_steps)
        
        # 收集微观结果
        micro_stresses = []
        for i, micro in enumerate(self.micro_models):
            strain = np.array([strains[i], 0, 0])
            stress, _ = micro.solve(strain)
            micro_stresses.append(stress)
        
        results = {
            'displacements': u,
            'strains': strains,
            'stresses': stresses,
            'micro_stresses': np.array(micro_stresses),
            'x_nodes': self.macro_model.x_nodes,
            'x_elements': (self.macro_model.x_nodes[:-1] + self.macro_model.x_nodes[1:]) / 2
        }
        
        return results


def visualize_multiscale_results(results, multiscale):
    """可视化多尺度结果"""
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    x_nodes = results['x_nodes']
    x_elem = results['x_elements']
    
    # 图1：位移分布
    ax1 = axes[0, 0]
    ax1.plot(x_nodes, results['displacements']*1e3, 'b-o', markersize=4)
    ax1.set_xlabel('位置 (m)')
    ax1.set_ylabel('位移 (mm)')
    ax1.set_title('宏观位移分布')
    ax1.grid(True, alpha=0.3)
    
    # 图2：应变分布
    ax2 = axes[0, 1]
    ax2.plot(x_elem, results['strains']*1e3, 'r-s', markersize=4)
    ax2.set_xlabel('位置 (m)')
    ax2.set_ylabel('应变 (×10⁻³)')
    ax2.set_title('宏观应变分布')
    ax2.grid(True, alpha=0.3)
    
    # 图3：应力分布
    ax3 = axes[1, 0]
    ax3.plot(x_elem, results['stresses']/1e6, 'g-^', markersize=4, label='宏观应力')
    ax3.plot(x_elem, results['micro_stresses'][:, 0]/1e6, 'm--', 
             markersize=4, label='微观应力(σxx)')
    ax3.set_xlabel('位置 (m)')
    ax3.set_ylabel('应力 (MPa)')
    ax3.set_title('应力分布对比')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 图4：微观结构示例
    ax4 = axes[1, 1]
    micro = multiscale.micro_models[0]
    im = ax4.imshow(micro.E_micro, extent=[-0.5, 0.5, -0.5, 0.5],
                    cmap='viridis', origin='lower')
    ax4.set_xlabel('ξ')
    ax4.set_ylabel('η')
    ax4.set_title('典型RVE微观结构')
    plt.colorbar(im, ax=ax4, label='相对刚度')
    
    plt.tight_layout()
    plt.savefig('multiscale_fe2_results.png', dpi=150, bbox_inches='tight')
    print("✓ 多尺度结果已保存: multiscale_fe2_results.png")
    plt.close()


# 主程序
if __name__ == "__main__":
    print("=" * 60)
    print("多尺度有限元计算 (FE²)")
    print("=" * 60)
    
    # 参数设置
    macro_length = 1.0  # m
    n_macro_elem = 10
    n_micro_points = 20
    
    print(f"\n模型参数:")
    print(f"  宏观结构长度: {macro_length} m")
    print(f"  宏观单元数: {n_macro_elem}")
    print(f"  微观网格点数: {n_micro_points}x{n_micro_points}")
    
    # 创建FE²框架
    print("\n初始化多尺度框架...")
    multiscale = MultiscaleFE2(macro_length, n_macro_elem, n_micro_points)
    
    # 求解
    print("\n开始多尺度计算...")
    results = multiscale.solve(load_steps=5)
    
    # 输出结果
    print("\n计算结果摘要:")
    print(f"  最大位移: {np.max(np.abs(results['displacements']))*1e3:.4f} mm")
    print(f"  最大应变: {np.max(np.abs(results['strains']))*1e6:.2f} με")
    print(f"  最大应力: {np.max(np.abs(results['stresses']))/1e6:.2f} MPa")
    
    # 可视化
    print("\n生成可视化结果...")
    visualize_multiscale_results(results, multiscale)
    
    print("\n" + "=" * 60)
    print("✓ 多尺度有限元计算完成!")
    print("=" * 60)

---

8. 总结与展望

8.1 本文工作总结

本文系统阐述了多尺度仿真与跨尺度耦合的理论方法和实现技术，主要贡献包括：

理论框架：

• 建立了完整的多尺度仿真框架，涵盖从量子到宏观的六个尺度层次
• 阐明了自下而上、自上而下和双向耦合三种信息传递机制
• 分类讨论了顺序、并发和自适应三种多尺度方法

方法体系：

• 详细介绍了各尺度的建模方法，包括DFT、MD、DDD、CPFEM和FEM
• 深入讨论了均匀化、准连续介质、桥接尺度等跨尺度耦合方法
• 介绍了FE²、晶体塑性、MD-FEM耦合等代表性方法

工程应用：

• 讨论了金属塑性、复合材料、断裂疲劳等领域的多尺度建模策略
• 展示了如何将微观机制与宏观响应关联
• 提供了材料设计和失效预测的多尺度思路

Python实现：

• 实现了分子动力学模拟，演示原子尺度模拟
• 开发了RVE分析程序，展示微观均匀化方法
• 构建了简化的FE²框架，演示并发多尺度耦合

8.2 技术发展趋势

8.2.1 机器学习与多尺度仿真

神经网络势函数：

• 机器学习拟合原子间相互作用
• 兼具DFT精度和MD效率
• 适合大尺度原子模拟

深度算子网络：

• 学习函数到函数的映射
• 实时多尺度分析
• 参数化微观结构优化

代理模型加速：

• 离线训练微观响应模型
• 在线快速预测
• 大幅降低计算成本

8.2.2 自适应与动态多尺度

误差驱动自适应：

• 基于误差估计动态调整尺度
• 计算资源优化配置
• 精度-效率平衡

数字孪生集成：

• 实时数据驱动的多尺度更新
• 在线模型修正
• 预测性维护

异构计算：

• GPU加速原子模拟
• 并行RVE计算
• 云计算平台

8.2.3 跨学科融合

材料基因组工程：

• 高通量多尺度计算
• 材料数据库建设
• 加速材料研发

生物力学应用：

• 骨组织多尺度建模
• 软组织本构
• 医疗器械设计

能源材料：

• 电池材料多尺度设计
• 储氢材料优化
• 太阳能电池

8.3 挑战与建议

8.3.1 当前挑战

计算效率：

• 尺度跨越巨大，计算量惊人
• 需要更高效的算法和硬件
• 降阶模型和代理模型

尺度间一致性：

• 不同尺度模型的匹配
• 边界条件的合理施加
• 信息传递的准确性

验证与确认：

• 实验数据获取困难
• 多尺度模型验证复杂
• 不确定性量化

8.3.2 发展建议

方法论层面：

• 发展自适应多尺度方法
• 建立多尺度验证标准
• 推进开源软件建设

应用层面：

• 聚焦关键工程问题
• 加强产学研合作
• 培养复合型人才

技术层面：

• 融合机器学习技术
• 利用异构计算资源
• 建设材料数据库

8.4 结语

多尺度仿真与跨尺度耦合是现代计算力学的核心发展方向，为理解材料和结构的力学行为提供了强有力的工具。通过将量子力学、分子动力学、介观方法和连续介质力学有机结合，我们能够从原子尺度到宏观尺度全面把握材料的变形、损伤和失效机制。

随着机器学习、自适应算法和异构计算技术的发展，多尺度仿真正朝着更高效、更智能的方向演进。数字孪生技术的兴起更为多尺度仿真提供了实时数据驱动的应用场景。可以预见，多尺度仿真将在新材料设计、结构优化、失效预测等领域发挥越来越重要的作用，为工程实践提供科学依据和技术支撑。

对于从事强度仿真和计算力学的研究人员和工程师而言，掌握多尺度仿真的理论方法和实现技术，将成为应对复杂工程挑战的关键能力。希望本文能够为读者提供系统的知识框架和实用的技术参考，推动多尺度仿真技术在中国的研究与应用。

---

参考文献

1. Liu, W.K., Karpov, E.G., & Park, H.S. (2006). Nano Mechanics and Materials: Theory, Multiscale Methods and Applications. Wiley.

2. Yip, S. (2005). Handbook of Materials Modeling. Springer.

3. Tadmor, E.B., & Miller, R.E. (2011). Modeling Materials: Continuum, Atomistic and Multiscale Techniques. Cambridge University Press.

4. Geers, M.G.D., Kouznetsova, V.G., & Brekelmans, W.A.M. (2010). Multi-scale computational homogenization: Trends and challenges. Journal of Computational and Applied Mathematics, 234(7), 2175-2182.

5. Raabe, D., Roters, F., & Barlat, F. (2004). Continuum Scale Simulation of Engineering Materials: Fundamentals - Microstructures - Process Applications. Wiley-VCH.

6. 杨卫, 赵亚溥. (2009). 多尺度力学问题. 中国科学: 物理学 力学 天文学, 39(12), 1689-1698.

7. 龙连春, 刘应华. (2015). 多尺度有限元方法及其应用. 力学进展, 45, 201502.

8. Feyel, F., & Chaboche, J.L. (2000). FE² multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials. Computer Methods in Applied Mechanics and Engineering, 183(3-4), 309-330.

9. Tadmor, E.B., Ortiz, M., & Phillips, R. (1996). Quasicontinuum analysis of defects in solids. Philosophical Magazine A, 73(6), 1529-1563.

10. Behler, J., & Parrinello, M. (2007). Generalized neural-network representation of high-dimensional potential-energy surfaces. Physical Review Letters, 98(14), 146401.

---

附录

附录A：常用势函数形式

Lennard-Jones势：
$$U(r)=4ϵ\left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right]$$

EAM势：
$$U=\sum _i{F}_{\alpha }\left(\sum _{j\ne i}{\rho }_{\beta }(r_{ij})\right)+\frac{1}{2}\sum _{i\ne j}{\varphi }_{\alpha \beta }(r_{ij})$$

Morse势：
$$U(r)=D_e{\left[1-e^{-a(r-r_e)}\right]}^2$$

附录B：周期性边界条件实现

最小像约定：
$$r_{ij}^{mic}=r_{ij}-L\cdot \text{round}\left(\frac{r_{ij}}{L}\right)$$

应变施加：
$${\mathbf{r}}^{\prime }=(\mathbf{I}+\mathbf{\epsilon })\cdot \mathbf{r}$$

附录C：均匀化方法对比

方法	精度	计算成本	适用场景
Voigt	低	极低	快速估算
Reuss	低	极低	快速估算
Hashin-Shtrikman	中	低	界限估计
自洽方法	中	中	多相材料
渐近均匀化	高	高	周期性结构
FE²	最高	最高	精确分析

附录D：开源多尺度软件

• LAMMPS: 大规模原子/分子并行模拟器
• OVITO: 原子结构可视化分析
• FEniCS: 有限元计算平台
• FEPX: 晶体塑性有限元
• DAMASK: 多尺度材料力学分析