2. Paraxial element theory#
In this part, we present the principle of paraxial approximation in the case of linear elastodynamics. Two theoretical approaches make it possible to identify the spirit and the practical application of elastic paraxial elements: the first is due to Cohen and Jennings [bib2] _ and the second is due to Modaressi [bib3] _. The application of the theory of paraxial elements to the fluid case will be made in the following part.
Throughout the following, as presented on r4.02.05-domaine_isfs, it is assumed that the border of the ground mesh is located in a domain with elastic behavior. The Modaressi approach implemented in code_aster makes it possible both to build absorbent boundaries and to introduce the incident seismic field.
2.1. Boundary spectral impedance#
To obtain the paraxial equation, we must first determine the shape of the diffracted displacement field in the vicinity of the border. To do this, we start from the equations of 3D elastodynamics:
- frac {{partial} ^ {2}disp} {partial {t} ^ {2}}
{tensTwo {E_ {11}}}} {soundSpeed} ^ {2}frac {partial^2disp} {partialposiTang^ {2}}} - {tensTwo {E_ {12}}}} {soundSpeed} ^ {2}frac {partial^2disp} {posiTangpartialposiNorm} - {tensTwo {E_ {22}}}} {soundSpeed} ^ {2}frac {partial^2disp} {partialposinorm^ {2}}} = VectorZero
With displacement \(\disp\) which breaks down into a normal and tangential part:
- available
= leftlbracebegin {array} {cc}
DispTang\ DispNorm
end {array}rightrbrace
and
The constant \(\soundSpeed\), which is homogeneous at one speed, is introduced to make certain quantities dimensionless. The equations and their solutions are of course independent of this constant.
We call \(\posiTang\) and \(\dispTang\) the directions and components of the movement in the tangential plane and \(\posiNorm\) and \(\dispNorm\) according to \(\basisVector{3}\), the normal direction at the border.
Two Fourier transforms are performed, one in relation to time, the other in relation to the space variables in the border plane. We limit ourselves to the case of a flat border without corners:
The equations are then written. According to the tangential component:
- left (SoundSpeedComp^ {2} -SoundSpeedCISA^ {2}right)
- left [
-vector {xi} _ {tau}cdotwidehat {dispTang} + i,frac {partialwidehat {dispNorm}} {partialposiNorm}
right],vector {xi} _ {tau} + SoundSpeedCISA^ {2} left [
normeUCL {vector {xi} _ {tau}} ^ {2} + frac {partial^2} {partialposiNorm^ {2}}
right],widehat {dispTang} + {omega} ^ {2},widehat {dispTang} = VectorZero
and according to the normal component
- left (SoundSpeedComp^ {2} -SoundSpeedCISA^ {2}right)
- left [
-I , vector {xi} _ {tau} cdot frac {partialwidehat {dispTang}} {partialPosiNorm} + frac {partial^2widehat {dispNorm}} {partialposiNorm^ {2}}
right] + SoundSpeedCISA^ {2} left [
normeUCL {vector {xi} _ {tau}} ^ {2} + frac {partial^2} {partialposiNorm^ {2}}
right],widehat {dispNorm} + {omega} ^ {2},widehat {dispNorm} = 0
where \(\widehat{\disp}\) and \(\widehat{\dispNorm}\) refer to Fourier transforms and \(\vector{\xi}_{\tau}\) the wave vector associated with \(\posiTang\).
It is a differential system in \(\posiNorm\) that we know how to solve by diagonalizing it. From this we deduce:
- leftlbrace
- begin {array} {ll}
frac {left (widehat {dispTang}wedgevector {xi} _ {tau}right),basisVector {3}} {normeUCL {vector {xi} _ {tau}}}} &=& A,expleft (-i, {xi} _ {S}\,posiNormright)\ {dispTangcdotvector {xi}} _ {tau}} &=&normeUCL {vector {xi} _ {tau}} _ {tau}}
- left [
{A} _ {P},exp (-i {xi} _ {P}posiNorm) + {A} _ {S}exp (-i {xi} _ {S}posiNorm)
right]\
normeUCL {vector {xi} _ {tau}}widehat {dispNorm} &=& - {A} _ {P}, {xi} _ {P} _ {P},expleft (-i {xi} _ {P}posiNormright) - {A} _ {S}\, {xi}\, {xi} _ {S}posiNormright)\
end {array}
right.
With \(\vector{\xi}_P=\sqrt{\frac{{\omega}^{2}}{\soundSpeedComp^{2}}-{\normEucl{\vector{\xi}_{\tau}}}^{2}}\) and \(\vector{\xi}_S=\sqrt{\frac{{\omega}^{2}}{\soundSpeedCisa^{2}}-{\normEucl{\vector{\xi}_{\tau}}}^{2}}\).
To determine the constants \(A\), \({A}_{S}\), and \({A}_{P}\), it is assumed that \(\dispTang \left(\vector{\xi}, 0\right)\) is known on the border of the finite element domain. These are expressed in terms of \(\widehat{\dispTang}\left(\vector{\xi}_{\tau}, 0\right)={\widehat{\dispTang}}_{,0}\) and \(\widehat{\dispNorm}\left(\vector{\xi}_{\tau},0\right)={\widehat{\dispNorm}}_{,0}\).
We will now evaluate the stress vector \(\sigmVector\) on a facet of normal \(\basisVector{3}\) in \(\posiNorm=0\), which will give us the impedance of the border. \(\sigmVector(\posiTang, \posiNorm)\) is subjected to the same Fourier transform in space as for the elastodynamic equations, so that:
- sigmVectorTang (vector {xi} _ {tau},posiNorm)
= left [
i,lameLambda,widehat {dispTang}cdotvector {xi} _ {tau} + (lameLambda +2laMemu)frac {partialwidehat {dispNorm}} {partialPosiNorm}} {partialPosiNorm}
right]basisVector {3} + LaMemuleft (
frac {partialwidehat {dispTang}} {partialPosiNorm} + i,widehat {dispNorm},posiTang
right)
\(\lameLambda\) and \(\lameMu\) are Lamé’s elastic coefficients.
In \(\posiNorm=0\), we want to get rid of terms containing derivatives in \(\posiNorm\). The system obtained previously allows us to do this according to \({\widehat{\dispTang}}_{,0}\) and \({\widehat{\dispNorm}}_{,0}\).
- leftlbrace
- begin {array} {ll}
frac {partialwidehat {dispTang} _ {,0}} {partialposiNorm},vector {xi} _ {tau} _ {tau} &=& i {normeUCL {vector {xi}} _ {tau}}}} ^ {2} {widehat {dispNorm} _ {,0}}}\ left (frac {partialwidehat {dispTang} _ {dispTang} _ {dispTang} _ {tau}right),basisVector {3}wedgevector {xi} _ {tau}right),basisVector {3} &=& -ivector {xi} _Sleft (widehat {dispTang} _ {,0}wedgevector {xi} _ {tau}right),basisVector {3}\ frac {partial {widehat {disp}}} _ {text {30}}} {partialPosiNorm} &=& ileft [-vector {xi} _Pvector {xi} _Svector {xi} _Sfrac {widehat {dispTang} _ {xi} _ {tau}} {{tau}}} {{tau}}} {{tau}},vector {xi}},vector {xi} _ {tau}},vector {xi} _P+vector {xi}}} {{tau}} +left (vector {xi}} _P+vector {xi}}} {tau}} +left (vector {xi}} _P+vector {xi}} _P+vector {xi}} xi} _Sright) {widehat {dispNorm} _ {,0}}right]\
end {array}
right.
The spectral impedance of the border is thus obtained:
- widehat {sigmVector} _ {0}
= {a} ^ {0}basisVector {3} + {b} ^ {0}vector {xi} _ {tau} + {c} ^ {0}vector {xi} _ {tau}wedgebasisVector {3}
where \({a}^{0}\), \({b}^{0}\), and \({c}^{0}\) are functions of \(\normEucl{\vector{\xi}_{\tau}}\) and \(\omega\) that depend linearly on \(\widehat{\dispTang}_{,0}\) and \({\widehat{\dispNorm}_{,0}}\).
We can then write:
- widehat {sigmVector} _ {0}:
= Aleft (
normeUCL {vector {xi} _ {tau}},omega
right) , {widehat {disp}} _ {0} left (vector {xi} _ {tau},omegaright)
where \(A\) refers to the global spectral impedance operator. We return to physical space by two inverse Fourier transforms.
2.2. Paraxial approximation of impedance#
The spectral impedance calculated previously is not local in space or time since it involves \(\left( \vector{\xi}_{\tau},\omega\right)\), the Fourier transform of \({\disp}_{0} \left(\posiTang, t\right)\) for all \(\posiTang\) and all \(t\).
The idea is then to develop \({\vector{\xi}}_{P}\) and \({\vector{\xi}}_{S}\) according to the powers of \(\frac{\normEucl{\vector{\xi}_{\tau}}}{\omega}\). This approximation will be good either at high frequencies or for \(\normEucl{\vector{\xi}_{\tau}}\) small.
Let’s look at dependency in \(\posiNorm\), for example from \(\widehat{\dispNorm}\). For \(\dispNorm\left(\posiTang, \posiNorm,t\right)\), we will have terms in the form \(\exp\left[i\left(\vector{\xi}_{\tau}\posiTang+\omega t-{\vector{\xi}}_{P}\posiNorm\right)\right]\).
With the development of \(\vector{\xi}_P\):
- vector {xi} _P
= frac {omega} {SoundSpeedComp} left [
left (frac {SoundSpeedCompnormeUCL {vector {xi} _ {tau}}} {omega}right) ^ {2} + dots
right]
We show that for \(\normEucl{\vector{\xi}_{\tau}}\) small, we will have waves propagating in directions similar to normal \(\basisVector{3}\) at the border, because the exponential is written:
- expleft{
i,omega left [
- left (
t-frac {posiNorm} {SoundSpeedComp}right) + i,mathcal {O}left (
frac {normeUCL {vector {xi} _ {tau}}} {omega}} {omega} right)
right]
right}
Therefore, with an asymptotic expansion of \({\vector{\xi}}_{P}\) and \({\vector{\xi}}_{S}\), by multiplying by a suitable power of \(\omega\) to remove this quantity from the denominator, we obtain:
where \({A}_{0}\) and \({A}_{1}\) are polynomial functions in \(\vector{\xi}_{\tau}\) and \(\omega\).
Or, after the two inverse Fourier transforms:
We thus obtain the final form of the local transient impedance approximated as a function of the last term in \(\frac{\normEucl{\vector{\xi}_{\tau}}}{\omega}\) selected. The detailed calculation of \({A}_{i}\) can be found in [bib4] _.
For example, for order 0:
This corresponds to viscous dampers distributed along the border of the finite element domain, with the corrective coefficient \(C^{a}\) entered behind the keyword COEF_AMOR assigned by [DEFI_MATERIAU] to the material of the absorbent boundary elements.
By analogy, for order 0, we can add a term in motion corresponding to rigidities distributed along the border of the finite element domain:
This additional formulation includes Lamé coefficients \(\lameLambda\) and \(\lameMu\), as well as a characteristic dimension \(L\).
In order 1:
- begin {multiline}
- frac {partial {sigmVector} _ {0}} {partial t} =
rhoSoundSpeedCompfrac {partial^2dispNorm} {partial {t} ^ {2}}basisVector {3} + rhoSoundSpeedCisafrac {partial^2dispTang} {partial {t} ^ {2}}} +\ rhoSoundSpeedCisaleft [
- left (2soundspeedCISA-soundspeedCompright),
frac {partial^2dispNorm} {partialposiTangpartial t}basisVector {3}
- left (SoundSpeedComp-2SoundSpeedCisaright),
frac {partial^2dispTang} {partialposiTangpartial t}right] +\
- rhoSoundSpeedComp^ {2}left (SoundSpeedCISA-frac {SoundSpeedComp} {2}right)
frac {partial^2dispNorm,basisVector {3}} {{partialposiTang} ^ {2}}} +
- rhoSoundSpeedCISA^ {2}left (SoundSpeedComp-frac {SoundSpeedCisa} {2}right)
frac {partial^2dispTang} {{partialposiTang} ^ {2}}
end {multiline}
We see the time derivative of the stress vector \(\sigmVector\) appearing. In digital processing, it will be necessary to use an integration of this term on border elements.
To conclude, it should be noted that the paraxial approximation leads to a transient local impedance essentially involving only derivatives in time and in the plane tangential to the border, as well as possibly an additional term without a time derivative.
Symbolically, we write in order 0:
and in order 1:
- frac {partial {sigmVector} _ {0}} {partial t}
= {A} _ {1}left (
frac {partial^2disp} {partial {t} ^ {2}}, frac {partialdisp} {partial t}, available
right)
2.3. Taking into account the incident seismic field#
It is recalled that the behavior of the ground is assumed to be elastic at least in the vicinity of the border. The diffracted field \({\disp}_r\) is therefore introduced such that:
- available
= {disp} _ {i} + {disp} _ {r}
Infinitely, the total field \(\disp\) must be equal to the incident field \({\disp}_{i}\) (one of the consequences of the Sommerfeld radiation condition), that is to say:
- lim_ {xto +infty}
{{disp} _ {r}} = VectorZero
At the border of the finite element mesh, we write the absorption condition for the diffracted field in order 0:
We deduce the total stress vector on the border of the finite element mesh:
We thus obtain the variational formulation of the problem in the vicinity of the border for order 0:
- rhoint_ {domain} {domain} {frac {partial^2disp} {partial {t} ^ {2}}dispVirtmeasDomain}
int_ {domain} {stress,left,left (dispright):strainleft (dispVirtright)measDomain} - int_ {bound} {{A} _ {0}\,left (frac {partialdisp} {partial t}right)dispVirtmeasBound} = int_ {bound} {
- left [
sigmVectorleft ({disp} _ {i}right) - {A} _ {0}left (frac {partial {disp} _ {i}} {partial t}right)
right],dispVirtmeasBound}
For any \(\dispVirt\) field that is kinematically eligible.
The absorption condition for the field diffracted to order 1 is written as:
- frac {partial {sigmVector} _ {0} ({disp} _ {r})} {partial t}
= {A} _ {1},left (
frac {partial^2 {disp} _ {r}} {partial {t} ^ {2}}}, frac {partial {disp} _ {r}} {partial t}, {disp} _ {r}
right)
For order 1, the conventional formulation is retained:
- rhoint_ {domain} {domain} {frac {partial^2disp} {partial {t} ^ {2}}dispVirtmeasDomain}
int_ {domain} {stressleft (uright):strainleft (dispVirtright)measDomain} - int_ {bound} {sigmVectorleft (dispright)dispVirtmeasBound} = 0
where \(\sigmVector \left(\disp\right)\) follows the following law of evolution:
- frac {partialsigmVectorleft (dispright)} {partial t}
=frac {partialsigmVectorleft ({disp} _ {i}right)} {partial t} + {A} _ {1}left (frac {partial^2disp} {partial {t} ^ {2}}},frac {partialdisp} {partialdisp} {partial^2disp},dispright) - {A} _ {1}left (frac {partial^2 {partial^2 {disp}} _ {i}} {partial} ^ {2}},frac {partial {disp} _ {disp} _ {disp} _ {disp}disp} _ {disp}right)
The solicitation due to the incident field appears explicitly in the case of order 0, but it is contained in the law of evolution of \(\sigmVector \left(\disp\right)\) for order 1.