2D Plane stress Problem, linear elasticity

[ruchitkini]

Dear @lululxvi ,

I am trying to solve a plane stress problem of Perforated Strip Plate as described in the following paper: https://link.springer.com/article/10.1007/s00466-022-02252-0

Below is my code:

"""
PINN Model for linear elasticity analysis of a Perforated strip plate using plane stress condition as in the paper
" A physics-informed neural network technique based on a modified loss function for computational 2D and 3D
solid mechanics " : https://doi.org/10.1007/s00466-022-02252-0
"""
import deepxde as dde

# Define Parameters
E = 5
nu = 0.3
lmbd = 0.333
mu = 0.385

# Define the geometry
rect = dde.geometry.Rectangle([0.0, 0.0], [1.0, 1.0])
circle = dde.geometry.Disk([0, 0], 0.2)
geom = dde.geometry.CSGDifference(rect, circle)


# Define the boundary functions
def boundary_left(x, on_boundary):
    return on_boundary and dde.utils.isclose(x[0], 0.0)


def boundary_right(x, on_boundary):
    return on_boundary and dde.utils.isclose(x[0], 1.0)


def boundary_bottom(x, on_boundary):
    return on_boundary and dde.utils.isclose(x[1], 0.0)


# Define boundary conditions
ux_left_bc = dde.icbc.DirichletBC(geom, lambda x: 0.0, boundary_left, component=0)
uy_bottom_bc = dde.icbc.DirichletBC(geom, lambda x: 0.0, boundary_bottom, component=1)
sx_right_bc = dde.icbc.DirichletBC(geom, lambda x: 1/E, boundary_right, component=2)


# Define governing equations
def pde(x, f):
    E_xx = dde.grad.jacobian(f, x, i=0, j=0)
    E_yy = dde.grad.jacobian(f, x, i=1, j=1)
    E_xy = 0.5 * (dde.grad.jacobian(f, x, i=0, j=1) + dde.grad.jacobian(f, x, i=1, j=0))

    S_xx = E_xx * (2 * mu + lmbd) + E_yy * lmbd
    S_yy = E_yy * (2 * mu + lmbd) + E_xx * lmbd
    S_xy = E_xy * 2 * mu

    Sxx_x = dde.grad.jacobian(f, x, i=2, j=0)
    Syy_y = dde.grad.jacobian(f, x, i=3, j=1)
    Sxy_x = dde.grad.jacobian(f, x, i=4, j=0)
    Sxy_y = dde.grad.jacobian(f, x, i=4, j=1)

    momentum_x = Sxx_x + Sxy_y
    momentum_y = Sxy_x + Syy_y

    stress_x = S_xx - f[:, 2:3]
    stress_y = S_yy - f[:, 3:4]
    stress_xy = S_xy - f[:, 4:5]

    return [momentum_x, momentum_y, stress_x, stress_y, stress_xy]


# Define data
data = dde.data.PDE(geom, pde, [ux_left_bc, uy_bottom_bc, sx_right_bc],
                    num_domain=25000, num_boundary=15000, num_test=30000)

# Define the Neural Network
layers = [2, [80] * 5, [80] * 5, [80] * 5, [80] * 5, [80] * 5, 5]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.PFNN(layers, activation, initializer)

# Define the model, optimiser and learning rate
model = dde.Model(data, net)
model.compile("adam", lr=0.0001)

# Save the model during training.
checkpointer = dde.callbacks.ModelCheckpoint("model/model", verbose=1, save_better_only=True)

losshistory, train_state = model.train(iterations=50000, callbacks=[checkpointer])
dde.saveplot(losshistory, train_state, issave=True, isplot=True)

Although the I get the loss in order of $10^{-8}$, the results are not at all correct as in the paper. I cannot find any reason for this.
It would be great if I can get some help from you. Is there any problem in my code? Or my boundary condition definition is not correct?

[whyzjuhit]

Could you please show some results for comparison to see the difference?

[ruchitkini]

I think I have also tried providing zero stress conditions in the y direction on the top boundary. But still no good results. What would be boundary conditions on the arc according to you? In the paper, it is not mentioned what boundary conditions have been used on the arc boundary.

[whyzjuhit]

I guess the arc stress could be extracted from FEM results. And you could also check if y direction is zero stress condition according to FEM simulation.

[ruchitkini]

Okay. I will try this out.

Just a general question. For NN training do I need to provide some extra stress conditions are better or a faster solution? Or the solution will ultimately converge if I provide many iterations. For example: In this problem, solving it on FEM software, Only the displacement boundary condition on the left and bottom boundaries are provided and the stress condition is only on the right boundary. Will the NN train the model with these 3 conditions?

[whyzjuhit]

I‘m working on this issue recently. According to my results, the main characteristic keeps well, for example, if you do a tensile test, in the tensile direction, you may get a reasonable result, however, in the other direction, may be the strain is not the same order of magnitude, to get a reasonable result, you may need extra constrain. I haven't found out why this happens, maybe I didn't set the NN well, maybe it's the mechanical problems. So if you want to train the model with these 3 conditions, in my advice, you may firstly know if it is solvable in mechanical perspective. As far as I know, extra information won't get a faster solution, on the contrary, if the points picked are too many or the pde contains second derivative, it becomes extremley slow.

[ruchitkini]

Hey, did you happen to get any progress on this issue? In one of the issues, I learned that the traction conditions should also be mentioned for other boundaries(#263). Also, there the stress conditions are implemented as OperatorBC, by specifying the stress equation. I have a doubt, doesn't the Dirichlet boundary condition calculate the same thing when the specific component is mentioned or I am wrong? @lululxvi , could you please help with this problem. How to formulate the stress conditions for the above mentioned problem.

@whyzjuhit, I also checked the code for the paper: https://github.com/JinshuaiBai/LSWR_loss_function_PINN/blob/main/03_perforated_strip_plate/PINN.py
If this could be helpful for solving the problem.

[ruchitkini]

Hello, I finally got the correct result by using hard constraints and applying suitable traction conditions as shown below.

I am leaving the code below if someone needs it.

import deepxde as dde
import numpy as np
import torch
from deepxde import config, optimizers

dde.config.set_default_float('float64')

# Define Parameters
E = 5
nu = 0.3
lmbd = nu / ((1 + nu) * (1 - nu))
mu = 0.5 / (1 + nu)

# Define the geometry
rect = dde.geometry.Rectangle([0, 0], [1.0, 1.0])
circle = dde.geometry.Disk([0, 0], 0.2)
geom = dde.geometry.CSGDifference(rect, circle)


# Define the boundary functions
def boundary_left(x, on_boundary):
    return on_boundary and dde.utils.isclose(x[0], 0.0)


def boundary_right(x, on_boundary):
    return on_boundary and dde.utils.isclose(x[0], 1.0)


def boundary_top(x, on_boundary):
    return on_boundary and dde.utils.isclose(x[1], 1.0)


def boundary_bottom(x, on_boundary):
    return on_boundary and dde.utils.isclose(x[1], 0.0)


def boundary_arc(x, on_boundary):
    return on_boundary and (x[0] < 0.2 and x[1] < 0.2)


def traction_arc11(x, f, _):
    theta = torch.atan2(x[:, 1:2], x[:, 0:1])
    s11 = (2 * mu + lmbd) * (dde.grad.jacobian(f, x, i=0, j=0)) + lmbd * (dde.grad.jacobian(f, x, i=1, j=1))
    s12 = mu * (dde.grad.jacobian(f, x, i=0, j=1) + dde.grad.jacobian(f, x, i=1, j=0))
    value11 = s11 * torch.cos(theta) + s12 * torch.sin(theta)

    return value11


def traction_arc22(x, f, _):
    theta = torch.atan2(x[:, 1:2], x[:, 0:1])
    s22 = (2 * mu + lmbd) * (dde.grad.jacobian(f, x, i=1, j=1)) + lmbd * (dde.grad.jacobian(f, x, i=0, j=0))
    s12 = mu * (dde.grad.jacobian(f, x, i=0, j=1) + dde.grad.jacobian(f, x, i=1, j=0))
    value22 = s12 * torch.cos(theta) + s22 * torch.sin(theta)

    return value22


# Define boundary conditions
s11_right_bc = dde.icbc.DirichletBC(geom, lambda x: 1.0/E, boundary_right, component=2)
s22_top_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_top, component=3)
s12_right_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_right, component=4)
s12_left_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_left, component=4)
s12_top_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_top, component=4)
s12_bottom_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_bottom, component=4)

s1_arc_bc = dde.icbc.OperatorBC(geom, traction_arc11, boundary_arc)
s2_arc_bc = dde.icbc.OperatorBC(geom, traction_arc22, boundary_arc)


# Define governing equations
def pde(x, f):
    E_xx = dde.grad.jacobian(f, x, i=0, j=0)
    E_yy = dde.grad.jacobian(f, x, i=1, j=1)
    E_xy = 0.5 * (dde.grad.jacobian(f, x, i=0, j=1) + dde.grad.jacobian(f, x, i=1, j=0))

    S_xx = E_xx * (2 * mu + lmbd) + E_yy * lmbd
    S_yy = E_yy * (2 * mu + lmbd) + E_xx * lmbd
    S_xy = E_xy * 2 * mu

    Sxx_x = dde.grad.jacobian(f, x, i=2, j=0)
    Syy_y = dde.grad.jacobian(f, x, i=3, j=1)
    Sxy_x = dde.grad.jacobian(f, x, i=4, j=0)
    Sxy_y = dde.grad.jacobian(f, x, i=4, j=1)

    momentum_x = Sxx_x + Sxy_y
    momentum_y = Sxy_x + Syy_y

    stress_x = S_xx - f[:, 2:3]
    stress_y = S_yy - f[:, 3:4]
    stress_xy = S_xy - f[:, 4:5]

    return [momentum_x, momentum_y, stress_x, stress_y, stress_xy]


# Define data
data = dde.data.PDE(
    geom,
    pde,
    [s11_right_bc, s22_top_bc, s12_right_bc, s12_left_bc, s12_top_bc, s12_bottom_bc,
     s1_arc_bc, s2_arc_bc],
    num_domain=4000,
    num_boundary=1000,
    num_test=20000)

# Define the Neural Network
layers = [2, [20] * 5, [20] * 5, [20] * 5, [20] * 5, [20] * 5, 5]
# layers = [2] + [60] * 5 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.PFNN(layers, activation, initializer)


def modify_output(X, f):
    x, y = X[:, 0:1], X[:, 1:2]
    ux, uy, sx, sy, sxy = f[:, 0:1], f[:, 1:2], f[:, 2:3], f[:, 3:4], f[:, 4:5]
    ux_new = x * ux
    uy_new = y * uy
    sx_new = sx
    sy_new = sy
    sxy_new = sxy
    return torch.cat((ux_new, uy_new, sx_new, sy_new, sxy_new), dim=1)


net.apply_output_transform(modify_output)

# Define the model, optimiser and learning rate
model = dde.Model(data, net)
model.compile("adam", lr=0.001, loss_weights=[1, 1, 1, 1, 1, 10, 10, 10, 10, 10, 10, 10, 10])

# Save the model during training.
checkpointer = dde.callbacks.ModelCheckpoint("model2/model", verbose=1, save_better_only=True)

# Train the model with Adam and then with BFGS
model.train(iterations=5000, display_every=100, callbacks=[checkpointer])
dde.optimizers.config.set_LBFGS_options(maxiter=25000)
model.compile("L-BFGS")
losshistory, train_state = model.train(callbacks=[checkpointer])

dde.saveplot(losshistory, train_state, issave=True, isplot=True)

RESULTS:

I want to ask @lululxvi , if there is any other way to implement traction conditions on the arc?
Also, for the earlier question, for plane surface, the traction given through OperatorBC and DirichletBC would be the same right?

Thanks.

[MinfeiLiang]

Thanks! I think there are mistakes with your calculation of the lame constant and shear modulus. Could you take a look?

[ruchitkini]

I did not multiply the lame's parameters by Young's modulus as I am scaling the equations by 1/E, same is done for traction value. Is this what you are talking about?

[MinfeiLiang]

Ok, I get what you mean with the scaling of 1/E. But the lmbd should equal to Enu / ((1 + nu) * (1 - 2nu)), right?

[ruchitkini]

Yes

[aililo-123]

Hi, may I ask why my calculations are completely wrong after changing Young's modulus of elasticity to 2.1e9 using your code?

[lululxvi]

I want to ask @lululxvi , if there is any other way to implement traction conditions on the arc?

Yes, you can use OperatorBC. You need to compute the normal direction inside OperatorBC.

for the earlier question, for plane surface, the traction given through OperatorBC and DirichletBC would be the same right?

OperatorBC can be used to implement everything.