材料力学：数据驱动预测位移

一、导入库

import numpy as np
import torch
import torch.nn as nn
import pandas as pd
import matplotlib.pyplot as plt 
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

二、构建神经网络 

class Net(nn.Module):

    def __init__(self, n_input, n_output, n_layer, n_nodes):
        super(Net, self).__init__()
        self.n_layer = n_layer
        
        self.Input = nn.Linear(n_input, n_nodes)   # linear layer
        nn.init.xavier_uniform_(self.Input.weight) # wigths and bias initiation
        nn.init.normal_(self.Input.bias)

        self.Output = nn.Linear(n_nodes, n_output)
        nn.init.xavier_uniform_(self.Output.weight)
        nn.init.normal_(self.Output.bias)
        
        self.Hidden = nn.ModuleList() # hidden layer list
        for i in range(n_layer):
            self.Hidden.append(nn.Linear(n_nodes, n_nodes))
        for layer in self.Hidden:
            nn.init.xavier_uniform_(layer.weight)
            nn.init.normal_(layer.bias)
        

    def forward(self, x):
        y = torch.tanh(self.Input(x)) # tanh activation function
        for layer in self.Hidden:
            y = torch.tanh(layer(y))
        y = self.Output(y)
        return y

三、定义求导函数

def derivative(x, Net):
    
    w = Net(x)  #输出 w 的形状是 (batch_size, n_output)
    dw_x = []
    dw_y = []
    #偏导数 dw_xy 是一个元组，dw_xy[0] 表示输出 w 对输入 x 的导数。由于输入 x 是二维的，所以 dw_xy[0][:,0] 和 dw_xy[0][:,1] 分别对应 x 的第一个分量和第二个分量的导数。
    for i in range(w.size()[1]):
        
        #print(Net(x),func(x).view(-1,1),w)
        dw_xy = torch.autograd.grad(w[:,i], x, torch.ones_like(w[:,i]), retain_graph=True,  create_graph=True,allow_unused=True)
        dw_x.append(dw_xy[0][:,0].view(-1,1)) #dw_x.append(dw_xy[0][:,0].view(-1,1))：将第 i 个输出分量相对于 x 第一个分量的导数添加到 dw_x 列表中。
        dw_y.append(dw_xy[0][:,1].view(-1,1))

    return w, dw_x, dw_y

四、定义均方误差函数

def err(X, Y):
    
    return torch.mean(torch.mean((X-Y)**2))

五、准备数据

class Dataset(torch.utils.data.Dataset):

    def __init__(self, X, Y):
        
        self.X = X
        self.Y = Y
        
    def __len__(self):
        
        return len(self.X)

    def __getitem__(self, index):
        
        x = self.X[index]
        y = self.Y[index]

        return x, y
data = pd.read_csv('rec-2.csv')
X_train = data.iloc[:, 5:7].to_numpy()
U = data.iloc[:,11:13].to_numpy()
LE = data.iloc[:,13:16].to_numpy()
sig = data.iloc[:,16:19].to_numpy()
U1 = U[:,0].reshape(-1, 1)
U2 = U[:,1].reshape(-1, 1)
eps11 = LE[:,0].reshape(-1, 1)
eps22 = LE[:,1].reshape(-1, 1)
eps12 = LE[:,2].reshape(-1, 1)
sig11 = sig[:,0].reshape(-1, 1)
sig22 = sig[:,1].reshape(-1, 1)
sig12 = sig[:,2].reshape(-1, 1)
X_train = torch.tensor(X_train, dtype=torch.float32)
X_train = X_train/10
U = torch.tensor(U, dtype=torch.float32)
U = U / 10
# Prepare training data
E, mu = 70, 0.3
Net_u = Net(2, 2, 10, 10)

六、训练模型

nepoches = 6000

learning_rate = 1.0e-3
optimizer = torch.optim.Adam(Net_u.parameters(), lr=learning_rate)
training_set = Dataset(X_train, U)
training_generator = torch.utils.data.DataLoader(training_set, batch_size= 128)
training_loss = []
sig_loss = []
for epoch in range(nepoches):
    
    for X_batch, Y_batch in training_generator:
        
        Y_pred = Net_u(X_batch)
    
        
        loss = err(Y_pred, Y_batch)
        loss.backward()
        
        optimizer.step()
        optimizer.zero_grad()
    if (epoch+1) % 100 == 0:
        print(f'epoch:{epoch+1}: loss:{loss:.4e} ')
        
    if (epoch+1) % 10 ==0:
        with torch.no_grad():
            Y_pred = Net_u(X_batch)
            loss = err(Y_pred, Y_batch)
            training_loss.append(loss)
            #使用导数计算得到的 du_x、du_y 和 dv_x、dv_y 计算信号 sig_x、sig_y 和 sig_xy。这里使用了一些物理公式（可能涉及材料科学或流体力学的背景）
        X_train.requires_grad = True
        U_pred, dux, duy = derivative(X_train, Net_u)
        du_x, du_y = dux[0], duy[0]
        dv_x, dv_y = dux[1], duy[1]
        
        sig_x = (du_x + mu*dv_y)*E/(1 - mu**2) 
        sig_y = (dv_y + mu*du_x)*E/(1 - mu**2)
        sig_xy = (dv_x + du_y)*E/(1 + mu)/2.
        
        sig11_loss = err(sig_x, torch.tensor(sig11))
        sig22_loss = err(sig_y, torch.tensor(sig22))
        sig12_loss = err(sig_xy, torch.tensor(sig12))
        sig_loss.append([sig11_loss,sig22_loss,sig12_loss])

七、相关性分析 

xx = [[xx[0].detach().numpy(),xx[1].detach().numpy(),xx[2].detach().numpy()] for xx in sig_loss]


# In[ ]:


pd.DataFrame(xx).to_clipboard()


# In[ ]:


X_train.requires_grad = True
U_pred, dux, duy = derivative(X_train, Net_u)

X = X_train[:,0].detach().numpy()*10
Y = X_train[:,1].detach().numpy()*10
U_pred = U_pred.detach().numpy()*10

u1 = U_pred[:,0].reshape(-1,1)
u2 = U_pred[:,1].reshape(-1,1)


# In[ ]:


corr_matrix1 = np.corrcoef(U1.reshape(1, -1), u1.reshape(1, -1))
print(corr_matrix1)

八、可视化

# Then, "ALWAYS use sans-serif fonts"
plt.rcParams['font.family'] = 'Times New Roman'
fig, ax = plt.subplots(figsize=(5.8, 4.8)) 
surf = ax.scatter(X, Y, c = u1, vmin=0, vmax=0.14, cmap=cm.rainbow)
#cbar = fig.colorbar(ax)
cb = fig.colorbar(surf)
cb.ax.locator_params(nbins=7)
cb.ax.tick_params(labelsize=16)
#cb.set_label(label =r'$\sigma_xx (MPa)$', fontsize=16)
#cb.set_label(fontsize=16)
ax.axis('equal')
ax.set_xlabel('X Position (mm)', fontsize=18)
ax.set_ylabel('Y Position (mm)', fontsize=18)
for tick in ax.get_xticklabels():
    #tick.set_fontname('Times New Roman')
    tick.set_fontsize(16)
for tick in ax.get_yticklabels():
    #tick.set_fontname('Times New Roman')
    tick.set_fontsize(16)
#plt.savefig('Flat-U-NN5-10.png', dpi=600, transparent=True)
plt.show()

九、总结

这段代码实现了一个基于神经网络的数据驱动模型，用于预测位移场和对应的应力分布。以下是代码的概述：

1. 模型架构

• 定义了一个全连接的前馈神经网络类（Net），其主要结构为：
  • 输入：二维坐标（x, y）。
  • 输出：位移分量（u 表示水平位移，v 表示垂直位移）。
  • 隐藏层：可配置的隐藏层和节点数量，权重使用 Xavier 均匀分布初始化。
  • 激活函数：使用 tanh 激活函数为非线性变换。

2. 偏导数计算

• derivative 函数计算预测位移（u, v）对输入（x, y）的偏导数，这对基于应变-位移关系计算应力是必要的。

3. 损失函数

• 使用均方误差（MSE）作为损失函数（err），计算预测位移与真实位移之间的误差。
• 此外，代码还计算应力分量（sig11，sig22，sig12）的损失，基于位移梯度和弹性理论中的应力-应变关系。

4. 数据准备

• 训练数据从 CSV 文件（rec-2.csv）中加载。输入数据是 x 和 y 坐标，而目标变量包括位移和应力分量。
• 数据集在训练前进行了归一化，并转换为 PyTorch 张量格式。

5. 训练过程

• 定义了自定义的 Dataset 类，用于处理训练的输入和输出数据对。
• 训练过程采用小批量数据进行迭代，使用 Adam 优化器优化网络权重。
• 每训练 100 轮时，打印训练损失；每训练 10 轮时，计算与应力相关的量并存储。

6. 后处理

• 训练完成后，生成整个数据集的位移预测，并计算预测位移与真实位移之间的相关性矩阵（U1 和 u1）。
• 同时，代码从位移的偏导数计算应力场（sig_x, sig_y, sig_xy），并通过散点图和颜色映射进行可视化。
• 还计算了预测剪切应力和实际剪切应力之间的相关性矩阵（sig12 和 sig_xy）。

7. 数据输出

• 代码将预测的位移和应力结果合并到一个输出数据集中，并将其保存为 CSV 文件（Flat-NN-10-10.csv）。

主要目标：

• 训练神经网络，基于二维空间数据（x, y）模拟位移场。
• 通过位移导数，使用应变-位移关系计算应力场。
• 可视化位移和应力的分布情况。
• 评估模型性能，通过比较预测和实际的位移及应力分量。

这个过程用于基于位移场和应力-应变关系，模拟材料和结构的力学行为。