# [Getdp] Implementation of PML for linear elastic waves

Julien de Rosny julien.derosny at espci.fr
Sat Sep 27 16:42:19 CEST 2014

```Hello,

I am trying to implement PML for elastic waves (but without success up
to now). I am starting with the classical coordinate transformation
introduced by Berenger

Because only D1 and D2 operator are implemented in getdp, I tried (but I
am not a specialist in finite elements) to use also curl operator and
perform a  transformation to express the PML with the D1 D2 curl operators

To that end I changed the Formulation proposed in the getdep
Elasticity_3D.pro example. See bellow.
A compressional wave is absorbed at the PML, but not a transverse one
that only involves c66.

To workout the C_m_xx matrix, I used Maxima to express  the stiffness
tensor in the D1 D2 and curl basis.

Any idea why it is is not properly working?

Thanks.

Regards.

Julien de Rosny
Insitut Langevin
1 rue Jussieu
75005 PARIS

Formulation {
{ Name Mec3D_u ; Type FemEquation ;
Quantity {
{ Name u ; Type Local ; NameOfSpace H_u_Mec3D ; }
}
Equation {
Galerkin { [L[]* C_m_11[] * Dof{D1 u} , {D1 u} ] ;
Galerkin { [L[]*  C_m_22[] * Dof{D2 u} , {D2 u} ] ;

Galerkin { [ L[]*C_m_32[] * Dof{D2 u} , {Curl u} ] ;
Galerkin { [ L[]*C_m_23[] * Dof{Curl u} , {D2 u} ] ;
Galerkin { [ L[]*C_m_33[] * Dof{Curl u} , {Curl u} ] ;

Galerkin { DtDt [ rho[] * Dof{u} , {u} ];
Galerkin { [ -F[] , {u} ] ;

}
}
}

Where

lambda=nu0*E0/((1+nu0)*(1-2*nu0)); //coef lamé
mu=E0/(2*(1+nu0)); //coef lamé - shear modulus
K=E0/3/(1-2*nu0);

c11 = lambda+2*mu; c12 = lambda; c66 = mu;

DampingProfileX[] = ( (X[]>=Lx) || (X[]<=(-Lx)) )?
1/(Lx+PmlThickness-Fabs[X[]] + epsilon)-1/(PmlThickness + epsilon) : 0;
DampingProfileY[] = ( (Y[]>=Ly) || (Y[]<=(-Ly)) )?
1/(Ly+PmlThickness-Fabs[Y[]] + epsilon)-1/(PmlThickness + epsilon) : 0;
DampingProfileZ[] = ( (Z[]>=Lz) || (Z[]<=(-Lz)) )?
1/(Lz+PmlThickness-Fabs[Z[]] + epsilon)-1/(PmlThickness + epsilon) : 0;

d1[] =Complex[1,  DampingProfileX[]];
d2[] =Complex[1,  DampingProfileY[]];
d3[] =Complex[1,  DampingProfileZ[]];

L[] = d1[]*d2[]*d3[];

C_m_11[   ] = Tensor[ c11/d1[]/d1[], c12/d1[]/d2[], c12/d1[]/d3[],
c12/d1[]/d2[], c11/d2[]/d2[], c12/d2[]/d3[],
c12/d1[]/d3[], c12/d2[]/d3[], c11/d3[]/d3[] ];

C_m_22[  ] = Tensor[ c66*(1/d1[]/d2[]+1/d1[]^2/2+1/d2[]^2/2), 0, 0,
0, c66*(1/d2[]/d3[]+1/d3[]^2/2+1/d2[]^2/2), 0,
0, 0, c66*(1/d1[]/d3[]+1/d3[]^2/2+1/d1[]^2/2) ];

C_m_32[  ] = Tensor[ 0, c66*(1/d2[]^2-1/d3[]^2)/2, 0,
0,                    0,
c66*(1/d3[]^2-1/d1[]^2)/2,
c66*(1/d1[]^2-1/d2[]^2)/2, 0,       0];

C_m_23[  ] = Tensor[ 0, 0,
c66*(1/d1[]^2-1/d2[]^2)/2,
c66*(1/d2[]^2-1/d3[]^2)/2, 0,       0,
0              ,
c66*(1/d3[]^2-1/d1[]^2)/2,                  0];

C_m_33[  ] = Tensor[ c66*(-1/d2[]/d3[]+1/d3[]^2/2+1/d2[]^2/2), 0, 0,
0, c66*(-1/d1[]/d3[]+1/d3[]^2/2+1/d1[]^2/2), 0,
0, 0,
c66*(-1/d2[]/d1[]+1/d1[]^2/2+1/d2[]^2/2) ];

If the
Galerkin { [L[]* C_m_11[] * Dof{D1 u} , {D1 u} ] ;
Galerkin { [L[]*  C_m_22[] * Dof{D2 u} , {D2 u} ] ;

*** MAXIMA CODE

PML:zeromatrix(9,9);PML[1,1]:1/l1;PML[2,2]:1/l2;PML[3,3]:1/l3;PML[4,4]:1/l2;PML[5,5]:1/l3;PML[6,6]:1/l1;PML[7,7]:1/l3;PML[8,8]:1/l1;PML[9,9]:1/l2;

C:matrix([c11,c12,c12,0,0,0,0,0,0],[C12,c11,c12,0,0,0,0,0,0],[C12,c12,c11,0,0,0,0,0,0
],[0,0,0,c66,0,0,c66,0,0],[0,0,0,0,c66,0,0,c66,0],[0,0,0,0,0,c66,0,0,c66],[0,0,0,c66,0,0,c66,0,0],[0,0,0,0,c66,0,0,c66,0],[0,0,0,0,0,c66,0,0,c66]);
P:zeromatrix(9,9);P[1,1]:1;P[2,2]:1;P[3,3]:1;P[4,6]:1;P[4,9]:1;P[5,4]:1;P[5,7]:1;P[6,5]:1;P[6,8]:1;P[7,4]:1;P[7,7]:-1;P[8,5]:1;P[8,8]:-1;P[9,6]:1;P[9,9]:-1;P;
Pi:invert(P);
Cget:expand(P.PML.C.PML.Pi);
C22:submatrix(1,2,3,7,8,9,Cget,1,2,3,7,8,9);C32:submatrix(1,2,3,4,5,6,Cget,1,2,3,7,8,9);C23:submatrix(1,2,3,7,8,9,Cget,1,2,3,4,5,6);C33:submatrix(1,2,3,4,5,6,Cget,1,2,3,4,5,6);

```