Coupled Piezoelectric Partial differential equations

[Rajat735]

Hello Lu,
Can this deep learning library solve this equation?

[lululxvi]

Ideally it can. But to achieve a good performance, it may require some hyperparameter tuning, etc.

[Rajat735]

I really appreciate it. Can you help me in defining this coupled PDE into the deep learning library.?
My schematic diagram is

I am considering the case when only force is applied on the cantilever beam.

[lululxvi]

Read the paper/video, also check the examples and the questions other people asked, and you will see that it is easy to code with DeepXDE. Let me know when you have any specific questions or difficulties.

[Rajat735]

Hello again, I am stuck in defining the p.d.e where it is double contracted with the displacement gradient. I am confused in matrix formation.

[Rajat735]

Is there any video as you mentioned? May you give me the link?

[lululxvi]

You can find the video in the document page. Could you provide more details about your problem of "double contracted"?

[Rajat735]

Thank you, Actually [c] is a fourth order tensor and gradient of u is a second order tensor so double contraction is to sum of two tensor orders and cancel it with twice the minimum of the two tensor order. It takes sum of the remaining two orders.

[lululxvi]

Is it c_{ijkl} u_k? I realized that the problem is 3D. I am not familiar with this PDE. If you have problems for the 3D PDE, I think a better way is to start from 1D or 2D cases.

[Rajat735]

Yes you get it right. I have to solve it for 2-d case only. I have to show strains in x and y direction plus the voltage variation when force is applied on the cantilever beam.

[Rajat735]

Hello lu, I am stuck in simplifying this p.d.e as c is a 3x3 matrix and u is a 2x2 matrix. May you help me with that?

[lululxvi]

Sorry, I am not familiar with this PDE, so I don't know the 2D form exactly. You may do some literature/textbook search. But my understanding is that in 2D, c should be a 2x2 matrix, and u is a scalar. You can re-write this original PDE in the scale form, and remove those terms with z.

[Rajat735]

Hello lu, I have simplified the PDE. It is still a little complex one with boundary conditions and partial derivatives. May you help me with that?

[Rajat735]

You were right about the c matrix matrix . Sorry, but I think u is a displacement vector .

[lululxvi]

@Rajat735 Could you post the simplified PDE/BCs, your code, and questions here?
@smao-astro @sanjusoni

[Rajat735]

Eq1

Eq2

Eq3

Refer to the above schematic diagram
BC
1: Left boundary is fixed
2: Potential is zero at the bottom
3: Normal derivative of stress is equal to traction force(some constant I will assume) at the right boundary

I have taken mechanical displacement as u=ux in x-direction + uy in y direction and potential as a function of x and y as f(x,y)

here i am stuck in implementing third boundary condition as stress is not the output components in PDE

I have written the code . Please correct it and scale also the PDE .

from future import absolute_import
from future import division
from future import print_function

import numpy as np
import tensorflow as tf

import deepxde as dde

def main():

def pde(x,y):
    c11 = 1.57 * 10**11
    c12 = 1.47 * 10**11
    c22 = c11
    e31 = -0.51
    ep_1 = 8.52
    ep_2 = ep_1
    Fx = 10
    
    ux,uy,f = y[:,0:1],y[:,1:2],y[:,2:3]
    dux_x = tf.gradients(ux, x)[0]
    duy_x = tf.gradients(uy, x)[0]
    df_x = tf.gradients(f , x)[0]
    
    dux_x, dux_y = dux_x[:, 0:1], dux_x[:, 1:2]
    duy_x, duy_y = duy_x[:, 0:1], duy_x[:, 1:2]
    df_x, df_y = df_x[:, 0:1], df_x[:, 1:2]
    
    dux_xy = tf.gradients(dux_x, x)[0][:, 1:2]
    duy_xy = tf.gradients(duy_x, x)[0][:, 1:2]
    
    
    dux_xx = tf.gradients(dux_x, x)[0][:, 0:1]
    dux_yy = tf.gradients(dux_y, x)[0][:, 1:2]
    
    duy_xx = tf.gradients(duy_x, x)[0][:, 0:1]
    duy_yy = tf.gradients(duy_y, x)[0][:, 1:2]
    
    df_xx = tf.gradients(df_x, x)[0][:, 0:1]
    df_yy = tf.gradients(df_y, x)[0][:, 1:2]

    return [  0.5*( 2*c11*dux_xx + c12*(dux_xy + dux_xx) + c11*(dux_yy + duy_xy) + 2*c12*duy_yy) + e31*df_xx + e31*df_yy - Fx,
          0.5*( 2*c12*dux_xx + c22*(dux_xy + dux_xx) + c12*(dux_yy + duy_xy) + 2*c22*duy_yy) + e31*df_xx + e31*df_yy ,
          ep_1*df_xx + ep_2*df_yy - e31*(2*dux_xx + dux_xy + duy_xx + duy_xy + dux_yy +2*duy_yy) ]

     
def func(x,y):
    c11 = 1.57 * 10**11
    c12 = 1.47 * 10**11
    c22 = c11
    e31 = -0.51
    ep_1 = 8.52
    ep_2 = ep_1

    
    ux,uy,f = y[:,0:1],y[:,1:2],y[:,2:3]
    dux_x = tf.gradients(ux, x)[0]
    duy_x = tf.gradients(uy, x)[0]
    df_x = tf.gradients(f , x)[0]
    
    dux_x, dux_y = dux_x[:, 0:1], dux_x[:, 1:2]
    duy_x, duy_y = duy_x[:, 0:1], duy_x[:, 1:2]
    df_x, df_y = df_x[:, 0:1], df_x[:, 1:2]
    return [ 0.5*(2*c11*dux_x + c12*(dux_y+duy_x))+ e31*df_x , 0.5*(2*c12*dux_x+c22*(dux_y+duy_x)) + e31*df_x ]







def boundary_l(x, on_boundary):

    return on_boundary and np.isclose(x[0], 0 ) 

def boundary_b(x, on_boundary):
    return on_boundary and np.isclose(x[1], 0 )

def boundary_nr(x, on_boundary):
    return on_boundary and np.isclose(x[0], 0.1)


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

bc_l = dde.DirichletBC(geom,lambda x: np.zeros((len(x), 1)), boundary_l,component=0)
bc_b = dde.DirichletBC(geom,lambda x: np.zeros((len(x), 1)), boundary_b,component=2)
bc_nr = dde.OperatorBC(geom ,  func, boundary_nr)


data = dde.data.PDE(geom,
                       3,
                     pde,
                     [bc_l,bc_b,bc_nr,],
                     num_domain=784,
                     num_boundary=216,
                     num_test=1500
                   )

layer_size = [2] + [50] * 2 + [3]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)

model = dde.Model(data, net)

model.compile("adam", lr=0.001)
losshistory, train_state = model.train(epochs=10000)

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

if name == "main":
main()

Third Boundary condition please help with it

[smao-astro]

I guess that you want to split your func (return list) to multiple functions which give the condition for only one output (return only one object), since after inspecting the code for method normal_derivative I found that it only works on one component every time. See

outputs = outputs[:, self.component : self.component + 1]

As a result, func which return a list won't be manipulated properly by NeumannBC/RobinBC. @lululxvi do you agree with that?

[Rajat735]

Thanks @smao-astro . Actually I want three outputs ux,uy which is mechanical displacement and potential function which is f(x,y) .

[lululxvi]

I agree with @smao-astro. @Rajat735 In your func, you have two outputs, 0.5*(2*c11*dux_x + c12*(dux_y+duy_x))+ e31*df_x and 0.5*(2*c12*dux_x+c22*(dux_y+duy_x)) + e31*df_x. Do you want these two to be a constant at the right boundary? Could you write down the third boundary condition in terms of ux, uy and f?

[Rajat735]

Yes @lululxvi the first is equal to some constant ,say 10, and second one is equal to zero. This is the third boundary condition which are expressed in terms of ux , uy and f.
Thank you @lululxvi @smao-astro for your effort.

[lululxvi]

You need to define the two BCs separately. For example, to define the first BC:

def func1(x, y, X):
    ...  # Your code in func related to this BC
    return 0.5 * (2 * c11 * dux_x + c12 * (dux_y + duy_x)) + e31*df_x - 10

bc_r1 = dde.OperatorBC(geom, func1, boundary_nr)

[Rajat735]

Thank you @lululxvi . Code is working fine now but i assume scaling is the problem. Please take a look at this image

Loss is not converging.

[lululxvi]

Yes, in the step 0, the loss is huge.

[Rajat735]

How to deal with this situation?

[lululxvi]

Check the FAQ here https://deepxde.readthedocs.io/en/latest/user/questions_answers.html

[Rajat735]

How to initialize loss weights and how to scale the pde such that each type of losses are in same range?
In #33 you commented about the solution in O(1). What it means ?
In output Graph (in Figure 2), What is the order in which the graph is plotted?

[lululxvi]

• Set loss weights, see Defining BCs #673

model.compile(..., loss_weights=[..., ..., ...])

• For scaling, see more are Intermediate Steps Included for Computation #35
• O(1) means the value of the solution is roughly around 1
• Which Figure 2 do you refer to?

[Rajat735]

3 graphs are generated while running the code. I am referring to the second one?

[lululxvi]

It only plots the first output, i.e., ux in your case.

To get more plots, you can use the output file test.dat.

[Rajat735]

The figure 2 does not shows the train.dat file

[lululxvi]

No, the train.dat is not plotted in the figure, because sometimes the plot would become messy. You can easily plot the figure using train.dat.

[Rajat735]

What if I have to impose a boundary condition in the direction of y axis? I have to impose a condition uy=y along the left boundary
Also i am confused with how many boundary conditions will be required to solve this P.d.e.?

[lululxvi]

• What do you mean boundary condition in the direction of y axis? Could you provide a detailed example?
• The number of BCs depneds on your PDE. Check some PDE book.