Implement layered Green's function

[sannant]

@pastewka, would it make sense to implement it as a separate Green's function object here ?

https://contactengineering.github.io/ContactMechanics/source/ContactMechanics.GreensFunctions.html?highlight=green+s#module-ContactMechanics.GreensFunctions

Then we would as a first step use it to compute the elastic energy for full contact by integration with the PSD.

On long term, or if needed we will probably want to interface PeriodicElasticHalspace to to this Green's function construct .

[pastewka]

Yes this makes sense. We need to update the half-space classes to accept varying Green's functions

[sannant]

@prs513rosewood implemented it in his code and pointed me again to the reference. But this is only plain-strain.

I think the result is the same for plain stress but we need a reference.

def normal_layered_green(L, h, shape, substrate, layer):
    """
    Plane-strain Green's function from Li & Popov (2020)
    https://doi.org/10.1177/1350650119854250
    Equations (30) and (31)
    """
    q = np.fft.fftfreq(shape[1], d=L[1] / shape[1]) * 2 * np.pi
    q = np.abs(q)

    E2, nu2 = substrate
    E1, nu1 = layer
    E_star = E1 / (1 - nu1**2)

    A = ((E2 * (3 - 4 * nu1) * (1 + nu1) - E1 * (3 - 4 * nu2) * (1 + nu2)) *
         (E1 * (1 + nu2) - E2 * (1 + nu1)))

    B = 4 * ((E2 * (1 + nu1) + E1 * (3 - 4 * nu2) * (1 + nu2)) *
             (E1 * (1 + nu2) - E2 * (1 + nu1)))

    C = (E1**2 * (4 * nu2 - 3) * (1 + nu2)**2 - 2 * E1 * E2 * (1 + nu1) *
         (2 * nu1 - 1) * (nu2 + 1) * (2 * nu2 - 1) + E2**2 *
         (8 * nu1**2 - 12 * nu1 + 5) * (1 + nu1)**2)

    D = ((E2 * (1 + nu1) + E1 * (3 - 4 * nu2) * (1 + nu2)) *
         (E2 * (3 - 4 * nu1) * (1 + nu1) + E1 * (1 + nu2)))

    with np.errstate(divide='ignore'):
        G = 2 / E_star * (
            (A * np.exp(-4 * q * h) + B * q * h * np.exp(-2 * q * h) + D) /
            (q *
             (-A * np.exp(-4 * q * h) - B * q**2 * h**2 * np.exp(-2 * q * h) +
              2 * C * np.exp(-2 * q * h) + D)))

    G[0] = 0

    return G

[sannant]

@prs513rosewood, this also works in 3D correct ?

https://github.com/sitangshugk95 needs that and might implement it in our code

[prs513rosewood]

@prs513rosewood, this also works in 3D correct ?

This works if you decompose the 2D DFT into plane waves and recompose after applying the Green's function, but I've never done that before. I'd try it on a homogeneous half-space first, with the plane strain homogeneous solution.

[sannant]

Ok thanks. I think that we just hace to replace q with sqrt(qx2 + qy2) and it should work, but we have to think that through.
As far as I understand, the 2D FFT is already a superposition of plane waves (with wavevectors of different orientations). I think we already had that discussion.

[prs513rosewood]

Given the isotropy and the normal load, q = |q| is a reasonable starting point. I can give you an independent 3d code if you're interested to have a reference simulation.