Tapered element issue?

[HaukeWittich]

Hello, I have the impression that there is a small bug in the shaft element part, more specifically in the second part of the stiffness matrix Kh2, which reflects the shear deformation. This part is depending on the surface A. For a tapered element the surface is variable described via a quadratic polynomial function with a constant A_l, which is the surface of the left hand side of the element. Further there is a transformation with phi, which is defined with the mean surface A in the middle of the element. The quotient between this two surfaces is not respected in my eyes. Thank you.

[raphaeltimbo]

Hi @HaukeWittich.
Thanks for looking at this. It is always good to have users double checking the code.
From what I understand, you are talking about this area defined here:

self.A = A_l * (1 + a1 * 0.5 + b1 * 0.5**2)

# which we could also calculate with

self.A = np.pi*(rom**2-rim**2)
# where rom is the outer radius mean between left and right nodes, and rim the inner.

Notice that this area is not a mean value between the area at the left node and the area in the right node.
This is the actual area in the middle of the element.

From this, I could not identify an error in the code.
Let me know if I am missing something here.

[HaukeWittich]

Hi
I guess, I become a real pain now, but I just saw the result in the code and unfortunately the correction there is the other way around. The code is:
K = E * Ie_l / (105 * L ** 3 * (1 + phi) ** 2) * (K1 + 105 * phi * K2 * A_l / A)
But that does make the situation even worse and the right correction is
K = E * Ie_l / (105 * L ** 3 * (1 + phi) ** 2) * (K1 + 105 * phi * K2 * A / A_l)
Sorry for being so persistant. Thank you.

[raphaeltimbo]

I am sorry for the error. I have performed the tests described in #880 in a jupyter notebook and messed up things when changing the actual code 🤦.
Should have added the tests to the test suite. Going to fix and do this now.
Thank you for your persistence!

[raphaeltimbo]

@HaukeWittich , I have just merged #883 to fix this.
Thanks again!

[HaukeWittich]

Hello, it is me again. I found another issue with hollow tapered elements also in the shear stiffness and would propose now a new approach:

[raphaeltimbo]

Hi @HaukeWittich.
Thank you once again for your contribution.
I have just changed the code with PR #888.
Will soon make a new release with the updated code.

[HaukeWittich]

Here is a test program with an alternative direct integration method:

import numpy as np

n = 6
x, w = np.polynomial.legendre.leggauss(n)


def HG(xi, Phi, l):
    xi2 = xi ** 2
    xi3 = xi ** 3
    k = 1 / (1 + Phi)

    H = np.zeros(4)
    G = np.zeros(4)

    H[0] = 2 * xi3 - 3 * xi2 + 1 + Phi * (1 - xi)
    H[1] = (xi3 - 2 * xi2 + xi + 0.5 * Phi * (xi - xi2)) * l
    H[2] = 3 * xi2 - 2 * xi3 + Phi * xi
    H[3] = (xi3 - xi2 + 0.5 * Phi * (xi2 - xi)) * l

    G[0] = (xi2 - xi) * 6 / l
    G[1] = 3 * xi2 - 4 * xi + 1 + Phi * (1 - xi)
    G[2] = (-xi2 + xi) * 6 / l
    G[3] = 3 * xi2 - 2 * xi + Phi * xi

    H = H * k
    G = G * k

    return H, G


def HGs(xi, Phi, l):
    xi2 = xi ** 2
    k = 1 / (1 + Phi)
    l2 = l * l

    Hs = np.zeros(4)
    Gs = np.zeros(4)

    Hs[0] = (6 * xi2 - 6.0 * xi - Phi) / l
    Hs[1] = 3 * xi2 - 4.0 * xi + 1.0 + 0.5 * Phi * (1.0 - 2.0 * xi)
    Hs[2] = (6 * xi - 6.0 * xi2 + Phi) / l
    Hs[3] = 3 * xi2 - 2.0 * xi + 0.5 * Phi * (2.0 * xi - 1.0)

    Gs[0] = (2 * xi - 1) * 6 / l2
    Gs[1] = (6 * xi - 4 - Phi) / l
    Gs[2] = (-2 * xi + 1) * 6 / l2
    Gs[3] = (6 * xi - 2 + Phi) / l

    Hs = Hs * k
    Gs = Gs * k

    return Hs, Gs


def TRfct(xi, Phi, l):
    H, G = HG(xi, Phi, l)

    T = np.zeros((2, 8))
    R = np.zeros((2, 8))

    T[0, 0] = H[0]
    T[1, 1] = H[0]
    T[1, 2] = -H[1]
    T[0, 3] = H[1]
    T[0, 4] = H[2]
    T[1, 5] = H[2]
    T[1, 6] = -H[3]
    T[0, 7] = H[3]

    R[0, 0] = G[0]
    R[1, 1] = -G[0]
    R[1, 2] = G[1]
    R[0, 3] = G[1]
    R[0, 4] = G[2]
    R[1, 5] = -G[2]
    R[1, 6] = G[3]
    R[0, 7] = G[3]

    return T, R


def TRsfct(xi, Phi, l):
    Hs, Gs = HGs(xi, Phi, l)

    Ts = np.zeros((2, 8))
    Rs = np.zeros((2, 8))

    Ts[0, 0] = Hs[0]
    Ts[1, 1] = Hs[0]
    Ts[1, 2] = -Hs[1]
    Ts[0, 3] = Hs[1]
    Ts[0, 4] = Hs[2]
    Ts[1, 5] = Hs[2]
    Ts[1, 6] = -Hs[3]
    Ts[0, 7] = Hs[3]

    Rs[0, 0] = Gs[0]
    Rs[1, 1] = -Gs[0]
    Rs[1, 2] = Gs[1]
    Rs[0, 3] = Gs[1]
    Rs[0, 4] = Gs[2]
    Rs[1, 5] = -Gs[2]
    Rs[1, 6] = Gs[3]
    Rs[0, 7] = Gs[3]

    return Ts, Rs


class ShaftElement:
    def __init__(self, L, idl, odl, idr=None, odr=None):
        if idr is None:
            idr = idl
        if odr is None:
            odr = odl

        self.L = float(L)
        self.idl = float(idl)
        self.odl = float(odl)
        self.idr = float(idr)
        self.odr = float(odr)

        self.A_l = np.pi * (odl ** 2 - idl ** 2) / 4
        self.Ie_l = np.pi * (odl ** 4 - idl ** 4) / 64

        roj = odl / 2
        rij = idl / 2
        rok = odr / 2
        rik = idr / 2

        # geometrical coefficients
        delta_ro = rok - roj
        delta_ri = rik - rij
        a1 = 2 * np.pi * (roj * delta_ro - rij * delta_ri) / self.A_l
        a2 = np.pi * (roj ** 3 * delta_ro - rij ** 3 * delta_ri) / self.Ie_l
        b1 = np.pi * (delta_ro ** 2 - delta_ri ** 2) / self.A_l
        b2 = (
            3
            * np.pi
            * (roj ** 2 * delta_ro ** 2 - rij ** 2 * delta_ri ** 2)
            / (2 * self.Ie_l)
        )
        gama = np.pi * (roj * delta_ro ** 3 - rij * delta_ri ** 3) / self.Ie_l
        delta = np.pi * (delta_ro ** 4 - delta_ri ** 4) / (4 * self.Ie_l)

        self.a1 = a1
        self.a2 = a2
        self.b1 = b1
        self.b2 = b2
        self.gama = gama
        self.delta = delta

        self.A = self.Af(0.5)
        self.Ie = self.If(0.5)

        self.rho = 7.85e-6
        Poisson = 0.3
        self.E = 206000
        self.G_s = self.E / 2 / (1 + Poisson)

        r = ((idl + idr) / 2) / ((odl + odr) / 2)
        r2 = r * r
        r12 = (1 + r2) ** 2

        self.kappa = (
            6
            * r12
            * ((1 + Poisson) / (r12 * (7 + 6 * Poisson) + r2 * (20 + 12 * Poisson)))
        )

        self.Phi = 12 * self.E * self.Ie / (self.G_s * self.kappa * self.A * L ** 2)

    def Af(self, xi):
        return self.A_l * (1 + self.a1 * xi + self.b1 * xi ** 2)

    def If(self, xi):
        return self.Ie_l * (
            1
            + self.a2 * xi
            + self.b2 * xi ** 2
            + self.gama * xi ** 3
            + self.delta * xi ** 4
        )

    def integrate(self):
        W = np.zeros((2, 8))
        Met = np.zeros((8, 8))
        Mer = np.zeros((8, 8))
        Ket = np.zeros((8, 8))
        Ker = np.zeros((8, 8))
        Ge = np.zeros((8, 8))

        h1 = 0.5
        h2 = 0.5

        for i in range(n):
            xn = h1 * x[i] + h2  # coordinate transformation
            T, R = TRfct(xn, self.Phi, self.L)
            Ts, Rs = TRsfct(xn, self.Phi, self.L)

            # Met
            h3 = self.rho * self.L * self.Af(xn) * w[i]
            h = np.matmul(T.T, T) * h3
            Met += h1 * h

            # Mer
            h3 = self.rho * self.L * self.If(xn) * w[i]
            h = np.matmul(R.T, R) * h3
            Mer += h1 * h

            # Ker
            h3 = self.E * self.L * self.If(xn) * w[i]
            h = np.matmul(Rs.T, Rs) * h3
            Ker += h1 * h

            # Ket
            h3 = self.G_s * self.L * self.kappa * self.Af(xn) * w[i]
            for j in range(2):
                for k in range(8):
                    if j == 0:
                        W[j, k] = Ts[j, k] - R[j, k]
                    else:
                        W[j, k] = Ts[j, k] + R[j, k]

            h = np.matmul(W.T, W) * h3
            Ket += h1 * h

            # Ges
            h3 = 2 * self.rho * self.L * self.If(xn) * w[i] * h1
            for j in range(8):
                for k in range(8):
                    h[j, k] = R[1, j] * R[0, k] - R[0, j] * R[1, k]
            Ge += h3 * h

        return Met, Mer, Ker, Ket, Ge, Met + Mer, Ket + Ker


se = ShaftElement(25, 8, 16.188212, 8, 25)
Met, Mer, Ker, Ket, Ge, M, K = se.integrate()

#%%
print(M[:4, :4])
print(K[:4, :4])
print(Ge[:4, :4])

[raphaeltimbo]

Thanks for sharing it @HaukeWittich .
I have edited your text just to format/highlight the python code.

[HaukeWittich]

Here a better readable version:
elements.pdf