Reference problem ===================== The objective of this test case is to compare the solution obtained with Code_Aster, using finite element modeling, with an analytical solution. This test case and its analytical solution were developed and detailed by Abousleiman et al. [1]. Here we present the main lines as well as the admissible analogies with our model. Geometry of the problem ---------------------- We consider a volume of rock with a height of :math:`h\mathrm{=}\mathrm{1m}` and a square base of :math:`l=\mathrm{15m}` on a side. In the center, a cylindrical cavity with radius :math:`R\mathrm{=}\mathrm{0.1m}` and infinite extension along its :math:`z` axis is hollowed out. We note that we are modeling here a domain of height :math:`\mathrm{1m}`, which is a restrictive hypothesis (ideally, we would need a much larger domain but which would lead to long calculation times). This cavity is used to model the well dug in the rock mass. This mass is also characterized by isotropy planes perpendicular to the axis of the cylindrical cavity. On the other hand, the well is characterized by two angles in relation to the reference frame of main constraints: * azimuth sound :math:`{\varphi }_{y}\mathrm{=}30°` * its :math:`{\varphi }_{z}\mathrm{=}60°` inclination Which corresponds to nautical angles: * :math:`\alpha =30°` * :math:`\beta =-60°` The Cartesian coordinate system associated with the mass and the stress coordinate system is noted :math:`(Ox\text{'}y\text{'}z\text{'})` and the one associated with the cavity is noted :math:`(\mathit{Oxyz})`. A cylindrical coordinate system :math:`(r,\theta ,z)` is also defined to characterize the cylindrical cavity. In Figure one summarizes the information described above: .. image:: images/Shape1.gif .. _RefSchema_Shape1.gif: **Figure** 1.1-1 **: Representation of the geometry of the inclined well problem** It is important to note that the two geometries described above are strictly equivalent. However, it is chosen to adopt the geometry for which the well is vertical, in order to facilitate the modeling and the design of the mesh. The coordinate system chosen is therefore the transverse isotropy coordinate system and not that of the main constraints. In the end, we will represent only a part of the model such as: .. image:: images/1000020000000400000001C1AA5C87B1F8ED0794.png :width: 4.9165in :height: 1.7165in .. _RefImage_1000020000000400000001C1AA5C87B1F8ED0794.png: **Figure** 1.1-2 **: Modeled Geometry** Symmetry conditions are therefore applied. Material properties -------------------- We refer to the parameters defined in Abu Sleiman's article :ref:`1 `. We start from the anisotropy ratio :math:`\frac{{E}_{L}}{{E}_{N}}\mathrm{=}2` and :math:`\frac{{\nu }_{\mathrm{TL}}}{{\nu }_{\mathrm{NL}}}=2`. We recall that the hypotheses of transverse isotropy imply that :math:`{\nu }_{T}={\nu }_{L}`; :math:`{\nu }_{N}={\nu }_{N}`; :math:`{\nu }_{L}={\nu }_{T}` and :math:`\frac{{\nu }_{N}}{{\nu }_{L}}=\frac{{E}_{N}}{{E}_{L}}`. In addition, and by writing the terms of the Hook matrix, such as: :math:`{M}_{\text{}11}=\frac{{E}_{\mathrm{L.}}({E}_{N}-{E}_{L}{\nu }_{N}^{2})}{(1+{\nu }_{T})({E}_{N}-{E}_{\mathrm{N.}}{\nu }_{T}-{2.E}_{\mathrm{L.}}{\nu }_{N}^{2})}` :math:`{M}_{\text{}12}=\frac{{E}_{L}({E}_{N}{\nu }_{T}+{E}_{L}{\nu }_{N}^{2})}{(1+{\nu }_{T})({E}_{N}-{E}_{\mathrm{N.}}{\nu }_{T}-{2.E}_{\mathrm{L.}}{\nu }_{N}^{2})}` :math:`{M}_{\text{}13}=\frac{{E}_{L}{E}_{N}{\nu }_{N}}{{E}_{N}-{E}_{N}{\nu }_{T}-{2.E}_{L}{\nu }_{N}^{2}}` :math:`{M}_{\text{}33}=\frac{{E}_{N}^{2}(1-{\nu }_{T})}{{E}_{N}-{E}_{N}{\nu }_{T}-{2.E}_{L}{\nu }_{N}^{2}}` The Biot coefficients are then deduced by the following relationship :ref:`1 `: :math:`{b}_{L}=1-\frac{{M}_{\text{}11}+{M}_{\text{}12}+{M}_{\text{}13}}{{\mathrm{3K}}_{S}}` and :math:`{b}_{N}=1-\frac{{\mathrm{2M}}_{\text{}13}+{M}_{\text{}33}}{{\mathrm{3K}}_{S}}` Here, :math:`{K}_{S}=\mathrm{27,5}\mathrm{GPa}`. From this, we indicate here the material properties entered in Code_Aster: .. csv-table:: "", "", "" "Liquid", "Viscosity :math:`{\mu }_{w}` (in :math:`\mathit{Pa.s}`) Compressibility module :math:`\frac{1}{{K}_{w}}` (in :math:`{\mathit{Pa}}^{\mathrm{-}1}`) Liquid density :math:`{\rho }_{w}` (in :math:`\mathit{kg}\mathrm{/}{m}^{3}`) Specific heat at constant pressure", ":math:`{10}^{\text{-3}}` :math:`{5.10}^{\text{-10}}` :math:`1` :math:`4180`" "Elastic parameters", "Young's modulus :math:`{E}_{L}` (in :math:`\mathit{MPa}`) Young's module :math:`{E}_{N}` (in :math:`\mathit{MPa}`) Poisson's Ratio :math:`{\nu }_{\text{LT}}` Poisson's Ratio :math:`{\nu }_{\text{LN}}` Shear Module :math:`{G}_{\mathrm{ln}}` (in :math:`\mathit{MPa}`)", ":math:`9474` :math:`0.5\mathrm{\times }9474` :math:`\mathrm{0,24}` :math:`\mathrm{0,24}` :math:`8880`" "Coupling parameters", "Biot coefficient :math:`{b}_{L}` Biot coefficient :math:`{b}_{N}` Initial homogenized density :math:`{r}_{0}` (in :math:`\mathit{kg}\mathrm{/}{m}^{3}`) Intrinsic permeability :math:`{K}_{L}^{\text{int}}` (in :math:`{m}^{2}\mathrm{/}s`) Intrinsic permeability :math:`{K}_{N}^{\text{int}}` (in :math:`{m}^{2}\mathrm{/}s`)", ":math:`\mathrm{0,81689}` :math:`\mathrm{0,89864}` :math:`2410` :math:`{5.10}^{\text{-20}}` :math:`{5.10}^{\text{-20}}`" **Table** 1.2-1 **: Material Properties** Boundary and initial conditions of the model --------------------------------------------- In the following, we consider the initial and boundary conditions in the coordinate system linked to the well. Travel cases ~~~~~~~~~~~~~~~~~~~~~~ Given the symmetry of the model, we will lock on 3 faces corresponding respectively to the 3 planes :math:`x=0;y=0` and :math:`z=0`. Pressure case ~~~~~~~~~~~~~~~~~~~~~ The initial pore pressure in the massif and for :math:`r\to +\mathrm{\infty }` is fixed at :math:`{p}_{0}\mathrm{=}\mathrm{9,8}\mathit{MPa}`. We will therefore apply this pressure to the 3 faces corresponding respectively to the 3 planes :math:`x=l;y=l` and :math:`z=h`. After the well is filled in, a uniform pressure develops on the wall, which is fixed at :math:`{p}_{w}\mathrm{=}12\mathit{MPa}` on the edge of the cavity. Case of constraints ~~~~~~~~~~~~~~~~~~~~~ In the massif, the stress state *in situ* associated with the :math:`(Ox\text{'}y\text{'}z\text{'})` coordinate system before the well was dug is characterized by the following initial state: :math:`\mathrm{\{}\begin{array}{c}S{\text{'}}_{x}\mathrm{=}\mathrm{-}25\mathit{MPa}\\ S{\text{'}}_{y}\mathrm{=}\mathrm{-}22\mathit{MPa}\\ S{\text{'}}_{z}\mathrm{=}\mathrm{-}29\mathit{MPa}\end{array}` In addition, digging the well involves a disturbance of the initial *in situ* stress state in the massif. In order to take into account this new state of stress in the vicinity of the well, it is necessary to consider two successive rotations to pass from the frame :math:`(Ox\text{'}y\text{'}z\text{'})` associated with the massif, to the frame :math:`(O\mathit{xyz})` associated with the well (cf. Figure). Let :math:`{P}_{1}` and :math:`{P}_{2}` be the passage matrices describing this change of frame of reference: :math:`{P}_{1}=(\begin{array}{ccc}\mathrm{cos}(\alpha )& -\mathrm{sin}(\alpha )& 0\\ \mathrm{sin}(\alpha )& \mathrm{cos}(\alpha )& 0\\ 0& 0& 1\end{array})` :math:`{P}_{2}=(\begin{array}{ccc}\mathrm{cos}(\beta )& 0& -\mathrm{sin}(\beta )\\ 0& 1& 0\\ \mathrm{sin}(\beta )& 0& \mathrm{cos}(\beta )\end{array})` After changing the coordinate system, the initial stress state in frame :math:`(\mathit{Oxyz})` and when :math:`r\to +\mathrm{\infty }` to be applied to the model is as follows: :math:`\{\begin{array}{c}{S}_{\mathrm{xx}}=-\mathrm{27,81}\mathrm{MPa}\\ {S}_{\mathrm{yy}}=-\mathrm{22,75}\mathrm{MPa}\\ {S}_{\mathrm{zz}}=-\mathrm{25,44}\mathrm{MPa}\\ {S}_{\mathrm{xy}}=\mathrm{0,65}\mathrm{MPa}\\ {S}_{\mathrm{xz}}=\mathrm{2,06}\mathrm{MPa}\\ {S}_{\mathrm{yz}}=\mathrm{1,125}\mathrm{MPa}\end{array}` For the boundary conditions, the shear is neglected and the pressure conditions are applied on the 3 faces :math:`x=l;y=l` and :math:`z=h`. At the wall of well :math:`(r\mathrm{=}R)` and the following, the radial stress is applied: :math:`{\sigma }_{\mathrm{rr}}(r=R,t)=-{p}_{w}\mathrm{.}H(t)` where :math:`H(t)` is the Heaviside step function. Benchmark solution --------------------- In this part we recall the theoretical equations useful for calculating radial stress and pore pressure. The boundary conditions of the general problem are summarized: * for :math:`r\to +\mathrm{\infty }` :math:`\{\begin{array}{c}{\sigma }_{\mathrm{xx}}=-{S}_{\mathrm{xx}}\\ {\sigma }_{\mathrm{yy}}=-{S}_{\mathrm{yy}}\\ {\sigma }_{\mathrm{zz}}=-{S}_{\mathrm{zz}}\\ {\sigma }_{\mathrm{xy}}=-{S}_{\mathrm{xy}}\\ {\sigma }_{\mathrm{xz}}=-{S}_{\mathrm{xz}}\\ {\sigma }_{\mathrm{yz}}=-{S}_{\mathrm{yz}}\end{array}` and :math:`p={p}_{0}` * for :math:`r\mathrm{=}R` :math:`{\sigma }_{\mathit{rr}}\mathrm{=}\mathrm{-}{p}_{w}\mathrm{.}H(t)` :math:`{\sigma }_{r\theta }\mathrm{=}{\sigma }_{\mathit{rz}}\mathrm{=}0` :math:`p\mathrm{=}{p}_{w}\mathrm{.}H(t)` Given the linearity of the problem, it is in fact the result of the superposition of 3 sub-problems that allow a three-dimensional study :ref:`1 `: * problem 1: we consider a poro-elastic problem in plane deformations (we only take into account the main stresses along the axes :math:`x` and :math:`y` as well as the shear in the plane :math:`(\mathit{xy})` and the disturbances related to the pore pressure). The boundary conditions for this problem are then written: * for :math:`r\to +\mathrm{\infty }` :math:`{\sigma }_{\mathrm{xx}}=-{S}_{\mathrm{xx}}` :math:`{\sigma }_{\mathrm{yy}}=-{S}_{\mathrm{yy}}` :math:`{\sigma }_{\mathit{xy}}\mathrm{=}\mathrm{-}{S}_{\mathit{xy}}` :math:`{\sigma }_{\mathit{zz}}\mathrm{=}\mathrm{-}{\nu }_{\text{LN}}({S}_{\mathit{xx}}+{S}_{\mathit{yy}})\mathrm{-}({\alpha }_{\text{LN}}\mathrm{-}2{\nu }_{\text{LN}}{\alpha }_{\text{LT}}){p}_{0}\mathrm{-}({\beta }_{3}^{s}\mathrm{-}2{\nu }_{\text{LN}}{\beta }_{1}^{s}){T}_{0}` :math:`{\sigma }_{\mathit{yz}}\mathrm{=}{\sigma }_{\mathit{xz}}\mathrm{=}0` :math:`p\mathrm{=}{p}_{0}` * for :math:`r\mathrm{=}R` :math:`{\sigma }_{\mathit{rr}}\mathrm{=}\mathrm{-}{p}_{w}\mathrm{.}H(t)` :math:`p\mathrm{=}{p}_{w}\mathrm{.}H(t)` * problem 2: we only consider the vertical constraint oriented according to :math:`z`. The boundary conditions for this problem are then written: * for :math:`r\to +\mathrm{\infty }` :math:`{\sigma }_{\mathrm{xx}}={\sigma }_{\mathrm{yy}}={\sigma }_{\mathrm{xy}}={\sigma }_{\mathrm{xz}}={\sigma }_{\mathrm{yz}}=p=0` :math:`{\sigma }_{\mathrm{zz}}=-{S}_{\mathrm{zz}}+({\nu }_{\text{LN}}({S}_{\mathrm{xx}}+{S}_{\mathrm{yy}})+({\alpha }_{\text{LN}}-2{\nu }_{\text{LN}}{\alpha }_{\text{LT}}){p}_{0})` :math:`{\beta }_{3}^{s}` and :math:`{\beta }_{1}^{s}` are defined later in this document. * for :math:`r\mathrm{=}R` :math:`{\sigma }_{\mathrm{rr}}={\sigma }_{r\theta }={\sigma }_{\mathrm{rz}}=p=0` * problem 3: we only consider the shears contained in planes :math:`(\mathit{xz})` and :math:`(\mathit{yz})`. The boundary conditions for this problem are then written: * for :math:`r\to +\mathrm{\infty }` :math:`{\sigma }_{\mathrm{xx}}={\sigma }_{\mathrm{yy}}={\sigma }_{\mathrm{zz}}={\sigma }_{\mathrm{xy}}=p=0` :math:`{\sigma }_{\mathit{xz}}\mathrm{=}\mathrm{-}{S}_{\mathit{xz}}` :math:`{\sigma }_{\mathit{yz}}\mathrm{=}\mathrm{-}{S}_{\mathit{yz}}` * for :math:`r\mathrm{=}R` :math:`{\sigma }_{\mathrm{rr}}={\sigma }_{r\theta }={\sigma }_{\mathrm{rz}}=p=0` In the rest of this document, the solution of problem 1 will be developed. Presentation of the analytical problem ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The expression of radial stress, pore pressure and temperature, for the problem in plane deformations (problem 1 above) is also the result of the superposition of 3 loading modes whose boundary conditions are given below: * Mode 1: taking into account the hydrostatic part of the boundary constraints (in :math:`r\mathrm{=}R`) :math:`{\sigma }_{\mathit{rr}}^{(1)}\mathrm{=}{\sigma }_{0}\mathrm{-}{p}_{w}` :math:`{p}^{(1)}\mathrm{=}0` * Mode 2: taking into account disturbances related to pore pressure and temperature (in :math:`r\mathrm{=}R`) :math:`{\sigma }_{\mathit{rr}}^{(2)}\mathrm{=}0` :math:`{p}^{(2)}\mathrm{=}{p}_{w}\mathrm{-}{p}_{0}` * Mode 3: taking into account the deviatory part of the boundary constraints (in :math:`r\mathrm{=}R`) :math:`{\sigma }_{\mathit{rr}}^{(3)}\mathrm{=}\mathrm{-}{S}_{0}\mathrm{cos}(2(\theta \mathrm{-}{\theta }_{r}))` :math:`{p}^{(3)}\mathrm{=}0` with: :math:`{\sigma }_{0}\mathrm{=}\frac{({S}_{\mathit{xx}}+{S}_{\mathit{yy}})}{2}` the mean stress, :math:`{S}_{0}\mathrm{=}0.5\sqrt{({({S}_{\mathit{xx}}\mathrm{-}{S}_{\mathit{yy}})}^{2}+4{S}_{\mathit{xy}}^{2})}` the deviatoric stress and :math:`\text{}{\theta }_{r}\mathrm{=}0.5\mathrm{arctan}(\frac{2{S}_{\mathit{xy}}}{{S}_{\mathit{xx}}\mathrm{-}{S}_{\mathit{yy}}})` * The general solution for calculating radial stress is therefore given by the following expression: :math:`{\sigma }_{\mathit{rr}}\mathrm{=}\mathrm{-}{\sigma }_{0}+{S}_{0}\mathrm{cos}(2(\theta \mathrm{-}{\theta }_{r}))+{\sigma }_{\mathit{rr}}^{(1)}+{\sigma }_{\mathit{rr}}^{(2)}+{\sigma }_{\mathit{rr}}^{(3)}` with: :math:`{\sigma }_{\mathit{rr}}^{(1)}\mathrm{=}({\sigma }_{0}\mathrm{-}{p}_{w})\frac{{R}^{2}}{{r}^{2}}` :math:`{\sigma }_{\mathit{rr}}^{(2)}\mathrm{=}{L}^{\mathrm{-}1}\left[\frac{1}{s}({b}_{L}(1\mathrm{-}\frac{{M}_{12}}{{M}_{11}})\left[{F}_{1}\psi (\xi )+{F}_{2}\psi (\omega )\right]+{\beta }_{1}^{s}(1\mathrm{-}\frac{{M}_{12}}{{M}_{11}})\left[({T}_{w}\mathrm{-}{T}_{0})\right]\psi (\omega ))\right]` :math:`{\sigma }_{\mathit{rr}}^{(3)}\mathrm{=}{L}^{\mathrm{-}1}\left[\frac{{S}_{0}}{s}\mathrm{cos}(2(\theta \mathrm{-}{\theta }_{r}))({A}_{1}{C}_{1}(\frac{1}{\xi r}{K}_{1}(\xi r)+\frac{1}{{(\xi r)}^{2}}{K}_{2}(\xi r))\mathrm{-}{A}_{2}{C}_{2}\frac{{R}^{2}}{{r}^{2}}\mathrm{-}3{A}_{3}{C}_{3}\frac{{R}^{4}}{{r}^{4}})\right]` The quantities used are defined at the end of this section. * The general solution for calculating the pore pressure is therefore given by the following expression: :math:`p\mathrm{=}{p}_{0}+{p}^{(2)}+{p}^{(3)}` with: :math:`{p}^{(2)}\mathrm{=}{L}^{\mathrm{-}1}\left[\frac{1}{s}\left[{F}_{1}\phi (\xi )+{F}_{2}\phi (\omega )\right]\right]` :math:`{p}^{(3)}\mathrm{=}{L}^{\mathrm{-}1}\left[\frac{{S}_{0}}{s}\mathrm{cos}(2(\theta \mathrm{-}{\theta }_{r}))({A}_{1}{C}_{2}\frac{{R}^{2}}{{r}^{2}}+\frac{{c}_{f}}{2{G}_{\text{LT}}\kappa }{C}_{1}{K}_{2}(\xi r))\right]` The quantities used are defined at the end of this section. On the other hand, :math:`{K}_{n}` is defined as being the modified Bessel functions of the second kind of order :math:`n` (where :math:`n` is an integer) such as: :math:`{K}_{n}(x)\mathrm{=}\begin{array}{c}\frac{1}{2}{(\frac{1}{2}x)}^{\mathrm{-}n}{\mathrm{\sum }}_{k\mathrm{=}0}^{n\mathrm{-}1}\left[\frac{(n\mathrm{-}k\mathrm{-}1)!}{k!}{(\frac{\mathrm{-}1}{4}{x}^{2})}^{k}\right]+{(\mathrm{-}1)}^{n+1}\mathrm{ln}\left[\frac{1}{2}x\right]{I}_{n}(x)\\ +{(\mathrm{-}1)}^{n}\frac{1}{2}{(\frac{1}{2}x)}^{n}{\mathrm{\sum }}_{k\mathrm{=}0}^{\mathrm{\infty }}\mathrm{[}(\chi (k+1)+\chi (n+k+1))\frac{{\left[\frac{1}{4}{x}^{2}\right]}^{k}}{k!(n+k)!}\mathrm{]}\end{array}` where :math:`{I}_{n}(x)\mathrm{=}\frac{1}{\pi }{\mathrm{\int }}_{0}^{\pi }{e}^{\mathit{xcos}\theta }\mathrm{cos}(n\theta )d\theta` are the modified Bessel functions of the first kind of order :math:`n` (where :math:`n` is an integer) and :math:`\chi` the digamma function. The following quantities can therefore be defined: :math:`\phi (x)\mathrm{=}\frac{{K}_{0}(\mathit{xr})}{{K}_{0}(\mathit{xR})}` :math:`\psi (x)\mathrm{=}\frac{{K}_{1}(\mathit{xr})}{{\mathit{xrK}}_{0}(\mathit{xR})}\mathrm{-}\frac{{\mathit{RK}}_{1}(\mathit{xr})}{{\mathit{xr}}^{2}{K}_{0}(\mathit{xR})}` :math:`{F}_{1}\mathrm{=}({p}_{w}\mathrm{-}{p}_{0})\mathrm{-}\frac{{c}_{\mathit{hf}}}{1\mathrm{-}\frac{{c}_{f}}{{c}_{h}}}({T}_{w}\mathrm{-}{T}_{0})` :math:`{F}_{2}\mathrm{=}\frac{{c}_{\mathit{hf}}}{1\mathrm{-}\frac{{c}_{f}}{{c}_{h}}}({T}_{w}\mathrm{-}{T}_{0})` :math:`{c}_{f}=\frac{\kappa M{M}_{11}}{{M}_{11}+{b}_{L}^{2}M}` :math:`\kappa =\frac{{K}_{L}^{\text{int}}(\phi )}{{\mu }_{w}}` :math:`{c}_{\mathit{hf}}=\frac{{c}_{f}}{\kappa }\left[{\beta }^{\mathit{sf}}-\frac{{b}_{L}{\beta }_{1}^{S}}{{M}_{11}}\right]` :math:`{C}_{1}=\frac{4}{{\mathrm{2A}}_{1}({B}_{3}-{B}_{2})-{A}_{2}{B}_{1}}` :math:`{C}_{1}=\frac{4}{{\mathrm{2A}}_{1}({B}_{3}-{B}_{2})-{A}_{2}{B}_{1}}` :math:`{C}_{3}=\frac{{\mathrm{2A}}_{1}({B}_{3}+{B}_{2})+{\mathrm{3A}}_{2}{B}_{1}}{3({\mathrm{2A}}_{1}({B}_{3}-{B}_{2})-{A}_{2}{B}_{1})}` :math:`{A}_{1}\mathrm{=}\frac{{b}_{L}M}{{M}_{11}+{b}_{L}^{2}M}` :math:`{A}_{2}\mathrm{=}\frac{{M}_{11}+{M}_{12}+2{b}_{L}^{2}M}{{M}_{11}+{b}_{L}^{2}M}` :math:`{B}_{1}\mathrm{=}\frac{{M}_{11}}{2{G}_{\text{LT}}{b}_{L}}{K}_{2}(\xi R)` :math:`{B}_{2}\mathrm{=}\frac{1}{\xi R}{K}_{1}(\xi R)+\frac{6}{{(\xi R)}^{2}}{K}_{2}(\xi R)` :math:`{B}_{3}\mathrm{=}2(\frac{1}{\xi R}{K}_{1}(\xi R)+\frac{3}{{(\xi R)}^{2}}{K}_{2}(\xi R))` :math:`\xi \mathrm{=}\sqrt{\frac{s}{{c}_{f}}}` and :math:`\omega \mathrm{=}\sqrt{\frac{s}{{c}_{h}}}` where :math:`s` represents the Laplace variable. Analytical solution resolution method ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ As we have seen previously, the solutions for modes 2 and 3 are expressed in Laplace space. Given the complexity of inverting the solution into time space, we use an inversion algorithm called *Stehfest-Graver.* This method is recognized for its numerical stability. The principle of the method consists in approximating the function to be inverted :math:`\tilde{f}(s)` using a delta-series using the following formula: :math:`f(t)\mathrm{\simeq }\frac{\mathrm{ln}(2)}{t}{\mathrm{\sum }}_{i\mathrm{=}1}^{N}{V}_{i}\tilde{f}(i\frac{\mathrm{ln}(2)}{t})` Where :math:`{V}_{i}\mathrm{=}{(\mathrm{-}1)}^{i+(N\mathrm{/}2)}\frac{{\mathrm{\sum }}_{k\mathrm{=}E((i+1)\mathrm{/}2)}^{\mathit{min}(i,N\mathrm{/}2)}{k}^{N\mathrm{/}2}(\mathrm{2k})!}{(\frac{N}{2}\mathrm{-}k)!k!(k\mathrm{-}1)!(i\mathrm{-}k)!(\mathrm{2k}\mathrm{-}i)!}` with :math:`E(x)` the integer part function and :math:`N` the number of terms in the series. It is generally between 10 and 20. In our case we chose :math:`N\mathrm{=}18`