2D uniaxial compression using linear elastic equation

[MinfeiLiang]

Hi all,

I am practicing deepxde for solving solid mechanics problems. I adapted the linear elastic example (https://deepxde.readthedocs.io/en/latest/demos/pinn_forward/elasticity.plate.html) to a simpler case: a uniaxial compression case of a 1-by-1 rectangular domain, with the bottom boundary fixed in x and y and an imposed deformation of -0.001m in y direction at the top boundary. However, even with such a simple case, I cannot get correct answers. Please see my code blow. Thanks a lot for the suggestions and help!

import torch
import deepxde as dde
import numpy as np
import matplotlib.pyplot as plt

E = 30
v = 0.2

lmbd = E*v/(1+v)/(1-2*v)
mu = 0.5*E/(1+v)

geom = dde.geometry.Rectangle([0, 0], [1, 1])

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)

ux_bottom_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_bottom, component=0)
uy_bottom_bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary_bottom, component=1)
uy_top_bc = dde.icbc.DirichletBC(geom, lambda x: -1e-3, boundary_top, component=1)

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]

data = dde.data.PDE(
    geom,
    pde,
    [
        ux_bottom_bc,
        uy_bottom_bc,
         uy_top_bc,
    ],
    num_domain=10000,
    num_boundary=10000
    )

layers = [2, [100] * 5, [100] * 5, [100] * 5, [100] * 5, 5]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.PFNN(layers, activation, initializer)

model = dde.Model(data, net)
model.compile("adam", lr=1e-4)

losshistory, train_state = model.train(epochs=10000,display_every=500)
dde.saveplot(losshistory, train_state, issave=True, isplot=True)

x = np.linspace(0, 1, 100)
y = np.linspace(0, 1, 100)
X, Y = np.meshgrid(x, y)
points = np.vstack((X.flatten(), Y.flatten())).T

# post-processing
predictions = model.predict(points)

ux = predictions[:, 0].reshape(X.shape)
uy = predictions[:, 1].reshape(X.shape)
Sxx = predictions[:, 2].reshape(X.shape)
Syy = predictions[:, 3].reshape(X.shape)
Sxy = predictions[:, 4].reshape(X.shape)

# displacement
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.contourf(X, Y, ux, levels=100, cmap='jet')
plt.colorbar()
plt.title('ux displacement')

plt.subplot(1, 2, 2)
plt.contourf(X, Y, uy, levels=100, cmap='jet')
plt.colorbar()
plt.title('uy displacement')

plt.show()

# stress
plt.figure(figsize=(18, 5))
plt.subplot(1, 3, 1)
plt.contourf(X, Y, Sxx, levels=100, cmap='jet')
plt.colorbar()
plt.title('Sxx stress')

plt.subplot(1, 3, 2)
plt.contourf(X, Y, Syy, levels=100, cmap='jet')
plt.colorbar()
plt.title('Syy stress')

plt.subplot(1, 3, 3)
plt.contourf(X, Y, Sxy, levels=100, cmap='jet')
plt.colorbar()
plt.title('Sxy stress')

plt.show()

[ww382]

You may also need to make sure appropriate traction conditions are applied.

[ruchitkini]

Hi, I have solved my problem and the issue was with the application of traction conditions. Maybe you can check your problem definition again. I referred to this issue #263.

[MinfeiLiang]

Hi,

I tried to apply the traction condition but it still does not solve the issue... I was wondering if you could share a snippet of your code related to the traction condition application. Seeing your approach would be helpful for me to understand where I might be going wrong.

[ruchitkini]

Hi, you can have a look at my complete code in issue #1570. I have applied traction as DirichletBC on plane boundaries and curved boundaries you need to use OperatorBC and define it via a function.

[MinfeiLiang]

Thanks! I also added the traction condition but it does not lead to the right solution. It seems scale of the mechanical constants (lmbd and mu) also influences a lot, because when I implemented your code with another set of mechanical constants (E=30000MPa and nu = 0.2), it leads to a problematic solution.

[ruchitkini]

Hey, yes for large parameters, scaling the way I did will make the loss very low and it will lose its importance. In this case, do not scale by E, rather keep it as it is and do the output transform.

[whyzjuhit]

Hi,
To my experience, I suggest that:

1. num_domain=10000, num_boundary=10000 is too large, 1000 is enough
2. lr=1e-4 is too small, 0.001 maybe set to 0.001
3. epochs=10000 is too small at least set to 50000-100000

[MinfeiLiang]

Thanks so much for your kind suggestions!