文章复现-PINN-薄板弯曲挠度预测-SciANN代码

参考文献:

[1] 唐明健, 唐和生. 基于物理信息的深度学习求解矩形薄板力学正反问题[J]. Chinese Journal of Computational Mechanics/Jisuan Lixue Xuebao, 2022, 39(1).

[2] Haghighat E, Juanes R. Sciann: A keras/tensorflow wrapper for scientific computations and physics-informed deep learning using artificial neural networks[J]. Computer Methods in Applied Mechanics and Engineering, 2021, 373: 113552.

1. 薄板弯曲理论

1.1 矩形薄板应变与挠度关系

 \varepsilon_x=-(\partial^2\ w)/(\partial x^2\ )\ z\\\varepsilon_y=-(\partial^2\ w)/(\partial y^2\ )\ z\\\gamma=-2\ (\partial^2\ w)/(\partial x\partial y)\ z\\

1.2 本构方程

\sigma_x=E/(1-\upsilon^2\ )\ (\varepsilon_x+\upsilon\varepsilon_y\ )\\ \sigma_y=E/(1-\upsilon^2\ )\ (\varepsilon_y+\upsilon\varepsilon_x\ )\\\tau=E/2(1+\upsilon)\ \ \gamma\\

1.3 薄板弹性曲面微分方程

D\nabla^4\ w=q\\

其中D=(E\delta^3)/(12(1-\mu^2\ ))

2. SciANN代码

2.1 导入第三方库

import numpy as np
import matplotlib.pyplot as plt 
import sciann as sn 
import time
from numpy import pi
from sciann.utils.math import diff, sin,cos

2.2 四边简支矩形薄板的纳维解(徐芝纶 弹性力学下 第五版)构建数据集

2.2.1 物性参数

E=2.1*10**5 #弹性模量
v=0.3       #泊松比
q=50        #均布荷载
a=1
b=1
d=1
D=E*d**3/(12.0*(1-v**2))

2.2.2 划分采样点(20X20)

x_data, y_data = np.meshgrid(
    np.linspace(0, 1, 20), 
    np.linspace(0, 1, 20)
)

x_data = x_data.reshape(-1,1)
y_data = y_data.reshape(-1,1)
xy = np.concatenate([x_data,y_data],axis=1)

2.2.3 计算理论解,构建数据集

w1 = []
fx1 = []
fy1 = []
fxy1 = []
x_real=np.linspace(0,1,20)
y_real=np.linspace(0,1,20)

for y in y_real:
    for x in x_real:
        w2 = 4*q / pi**6 / D *np.sin(pi * x)*np.sin(pi * y)          ##位移
        fx2 = 24*q/pi**4/d**2*(1+v)*np.sin(pi * x)*np.sin(pi * y)    ##x方向正应力
        fy2 = 24*q/pi**4/d**2*(1+v)*np.sin(pi * x)*np.sin(pi * y)    ##y方向正应力
        fxy2 = -24*q/pi**4/d**2*(1-v)*np.cos(pi * x)*np.cos(pi * y)  ##切应力
        w1.append(w2)
        fx1.append(fx2)
        fy1.append(fy2)
        fxy1.append(fxy2)

w1= np.array(w1).reshape(-1,1)
fx1=np.array(fx1).reshape(-1,1)
fy1= np.array(fy1).reshape(-1,1)
fxy1= np.array(fxy1).reshape(-1,1)

##归一化
w3=(w1-w1.mean())/w1.std()
fx3=(fx1-fx1.mean())/fx1.std()
fy3=(fy1-fy1.mean())/fy1.std()
fxy3=(fxy1-fxy1.mean())/fxy1.std()

##组合到一起
train_1 = np.concatenate([w3,fx3],axis=1)
train_2 = np.concatenate([train_1,fy3],axis=1)
train_3 = np.concatenate([train_2,fxy3],axis=1)
train_data=np.concatenate([xy,train_3],axis=1)

2.3 构建神经网络

2.3.1 输入

x = sn.Variable('x',dtype='float64')
y = sn.Variable('y',dtype='float64')

2.3.2 输出

设置3个隐藏层,每个隐藏层40个神经单元,激活函数为双曲正切函数。

sigam_x1 = sn.Functional('sigam_x', [x,y], 3*[40], 'tanh')
sigam_y1 = sn.Functional('sigam_y', [x,y], 3*[40], 'tanh')
sigam_xy1 = sn.Functional('sigam_xy', [x,y], 3*[40], 'tanh')

2.4 损失函数

2.4.1 归一化复原

def invstd(X, mu, sig):
    return X*sig + mu

w = invstd(w_1, w1.mean(), w1.std()) 
sigam_x = invstd(sigam_x1, fx1.mean(), fx1.std())
sigam_y = invstd(sigam_y1, fy1.mean(), fy1.std())
sigam_xy = invstd(sigam_xy1, fxy1.mean(), fxy1.std())

WBAR = np.min([fx1.std(),fy1.std(),fxy1.std()])

2.4.2 求导设置

z=d/2.0
w_x1 = diff(w,x)
w_y1 = diff(w,y)

w_x2 = diff(w,x,order=2)
w_y2 = diff(w,y,order=2)
w_xy = diff(w_x1,y)

u_x=- z *w_x2
u_y=- z *w_y2
u_xy = -2.0*z*w_xy

w_x4=diff(w, x, order=4)
w_y4=diff(w, y, order=4)
w_x2y2 = diff(w_x2, y ,order=2)

2.4.3 构建损失函数

将损失函数在基于数据驱动的神经网络模型基础上添加薄板的弹性曲面微分方程,为了表征尽可能多的物理信息,又将应力应变的本构关系式纳入到损失函数中。

L1 = sn.Data(w_1)  ##挠度
L2 = sn.Data(sigam_x1)   ##x方向正应力
L3 = sn.Data(sigam_y1)   ##x方向正应力
L4 = sn.Data(sigam_xy1)  ##切应力
L5 = sn.Tie(-w_x4,w_y4+2.0*w_x2y2-q/D)  ##弹性曲面微分方程

##本构方程
L6 = sn.Tie(sigam_x/WBAR,E/(1.0-v**2)*(u_x+v*u_y)/WBAR)
L7 = sn.Tie(sigam_y/WBAR,E/(1.0-v**2)*(u_y+v*u_x)/WBAR)
L8 = sn.Tie(sigam_xy/WBAR,E/(2.0*(1+v))*u_xy/WBAR)

2.5 优化模型

model = sn.SciModel(
    inputs=[x, y],
    targets=[L1, L2, L3, L4, L5,L6,L7,L8],
    loss_func="mse",
    plot_to_file='Ban.png'
)

2.6 训练

data_L1 = train_data[:,2]
data_L2 = train_data[:,3]
data_L3 = train_data[:,4]
data_L4 = train_data[:,5]
data_L5 = 'zeros'
data_L6 = 'zeros'
data_L7 = 'zeros'
data_L8 = 'zeros'
target_data = [data_L1, data_L2, data_L3, data_L4,data_L5,data_L6,data_L7,data_L8]

input_data = [train_data[:,0], train_data[:,1]]

t = time.time()
history = model.train(
    x_true=input_data,
    y_true=target_data,
    epochs=2500,
    batch_size=32,
    shuffle=True,
    learning_rate=0.001,
    reduce_lr_after=100,
    stop_loss_value=1e-6,
    verbose=0
)
t = time.time() - t
  • 21
    点赞
  • 23
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
好的,以下是用 PINN 求解二阶 Diffusion 方程的 PyTorch 代码: ```python import torch import numpy as np import matplotlib.pyplot as plt # 定义求解区域的左右边界 x_left = 0.0 x_right = 1.0 # 定义时间和空间步长 delta_t = 0.001 delta_x = 0.05 # 定义神经网络的输入维度和隐藏层维度 input_dim = 2 hidden_dim = 10 # 定义神经网络的输出维度 output_dim = 1 # 定义 Diffusion 方程的初始条件和边界条件 def u_init(x): return np.sin(np.pi * x) def u_left(t): return 0.0 def u_right(t): return 0.0 # 定义神经网络模型 class PINN(torch.nn.Module): def __init__(self): super(PINN, self).__init__() self.fc1 = torch.nn.Linear(input_dim, hidden_dim) self.fc2 = torch.nn.Linear(hidden_dim, hidden_dim) self.fc3 = torch.nn.Linear(hidden_dim, output_dim) def forward(self, x, t): # 拼接时间和空间坐标作为神经网络的输入 xt = torch.cat([x, t], dim=1) h1 = torch.tanh(self.fc1(xt)) h2 = torch.tanh(self.fc2(h1)) out = self.fc3(h2) return out # 定义 PINN 的损失函数 def pinn_loss(model, x_left, x_right, delta_t, delta_x): # 定义空间和时间的网格 x = torch.linspace(x_left, x_right, int((x_right - x_left) / delta_x) + 1).unsqueeze(1) t = torch.linspace(0, 1, int(1 / delta_t) + 1).unsqueeze(1) # 计算模型在内部网格点上的预测值 x_internal = x[1:-1] t_internal = t[1:-1] u_internal = model(x_internal, t_internal) # 计算模型在左右边界上的预测值 u_left_boundary = model(x_left, t_internal) u_right_boundary = model(x_right, t_internal) # 计算初始条件和边界条件的损失 init_loss = torch.mean((model(x, torch.zeros_like(x)) - u_init(x)) ** 2) left_boundary_loss = torch.mean((u_left_boundary - u_left(t_internal)) ** 2) right_boundary_loss = torch.mean((u_right_boundary - u_right(t_internal)) ** 2) # 计算 Diffusion 方程的残差 u_x = torch.autograd.grad(u_internal, x_internal, grad_outputs=torch.ones_like(u_internal), create_graph=True)[0] u_xx = torch.autograd.grad(u_x, x_internal, grad_outputs=torch.ones_like(u_x), create_graph=True)[0] f = u_xx - (1 / delta_t) * (u_internal - u_init(x_internal)) # 计算残差的损失 residual_loss = torch.mean(f ** 2) # 计算总的损失 loss = init_loss + left_boundary_loss + right_boundary_loss + residual_loss return loss # 初始化神经网络模型和优化器 model = PINN() optimizer = torch.optim.Adam(model.parameters(), lr=0.001) # 训练神经网络模型 for i in range(10000): optimizer.zero_grad() loss = pinn_loss(model, x_left, x_right, delta_t, delta_x) loss.backward() optimizer.step() if i % 1000 == 0: print("Iteration {}, Loss = {}".format(i, loss.item())) # 使用训练好的神经网络模型进行预测 x_test = torch.linspace(x_left, x_right, 101).unsqueeze(1) t_test = torch.linspace(0, 1, 1001).unsqueeze(1) u_test = model(x_test, t_test).detach().numpy() # 绘制预测结果 X, T = np.meshgrid(x_test, t_test) plt.pcolormesh(X, T, u_test, cmap='coolwarm') plt.colorbar() plt.xlabel("x") plt.ylabel("t") plt.title("PINN Solution to Diffusion Equation") plt.show() ``` 这段代码中,首先定义了求解区域的左右边界、时间和空间步长等参数,然后定义了 Diffusion 方程的初始条件和边界条件,以及神经网络模型和损失函数。在训练神经网络模型时,使用了 Adam 优化器,最后使用训练好的模型进行预测,并绘制了预测结果。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值