Tan lib

tan_lib

Functions:

Calc_S12_eff

Calc_S12_eff(C_strat: ndarray) -> tuple[complex, complex]

Return the effective s1 and s2 for a given stiffness matrix.

Source code in src\vimseo\lib_vimseo\tan_lib.py

def Calc_S12_eff(C_strat: np.ndarray) -> tuple[complex, complex]:
    """Return the effective s1 and s2 for a given stiffness matrix."""

    # Determination of the appropriate flexibility matrix

    S_strat_omega, omega = Calc_S_matrix(C_strat)

    # Solving characteristic equation

    s1_st, s2_st = S_12(S_strat_omega)

    s1 = (s1_st * np.cos(omega) + np.sin(omega)) / (
        np.cos(omega) - s1_st * np.sin(omega)
    )

    s2 = (s2_st * np.cos(omega) + np.sin(omega)) / (
        np.cos(omega) - s2_st * np.sin(omega)
    )

    # Regularise the isotropic double-root singularity (see _separate_roots).
    return _separate_roots(s1, s2)

Calc_S_matrix

Calc_S_matrix(C_strat: ndarray) -> list[ndarray, float]

Return the effective compliance matrix and the angle of the orthotropic axes.

Source code in src\vimseo\lib_vimseo\tan_lib.py

def Calc_S_matrix(C_strat: np.ndarray) -> list[np.ndarray, float]:
    """Return the effective compliance matrix and the angle of the orthotropic axes."""

    # Degenerate (in-plane isotropic) case: when both the extension-shear
    # coupling and the C11/C22 asymmetry are at the rounding-noise level relative
    # to the stiffness magnitude, the orthotropy axes are undefined -- every
    # rotation leaves the stiffness invariant. The angle otherwise computed as
    # arctan(coupling / asymmetry) is then a 0/0 ratio of noise, not reproducible
    # across platforms/BLAS. We make the arbitrary but deterministic choice
    # omega = 0 (keep the current frame), which is physically exact here.
    scale = np.linalg.norm(C_strat)

    coupling = np.abs(C_strat[0, 2]) + np.abs(C_strat[1, 2])

    asymmetry = np.abs(C_strat[0, 0] - C_strat[1, 1])

    if coupling <= _ISOTROPY_RTOL * scale and asymmetry <= _ISOTROPY_RTOL * scale:
        C_final = C_strat.copy()

        C_final[-1, :2] = np.zeros(2)

        C_final[:2, -1] = np.zeros(2)

        return np.linalg.inv(C_final), 0.0

    # Identification of possible reference points

    if np.abs(C_strat[0, 0] - C_strat[1, 1]) > 1e-9:
        frac = 2 * (C_strat[0, 2] + C_strat[1, 2]) / (C_strat[0, 0] - C_strat[1, 1])

        omega1 = np.arctan(frac) / 2

    else:
        omega1 = np.pi / 4

    denom = np.abs(
        C_strat[0, 0] + C_strat[1, 1] - 2 * C_strat[0, 1] - 4 * C_strat[2, 2]
    )

    if denom > 1e-9:
        frac_bis = 4 * (C_strat[0, 2] + C_strat[1, 2]) / denom

        omega2 = np.arctan(frac_bis) / 4

    else:
        omega2 = np.pi / 8

    # Possible rotation matrices

    T_omega_inv1 = Mat_rot(omega1)

    T_omega_inv2 = Mat_rot(omega2)

    # stiffness matrices in rotated coordinate systems

    C_try1 = np.linalg.multi_dot([T_omega_inv1, C_strat, T_omega_inv1.T])

    C_try2 = np.linalg.multi_dot([T_omega_inv2, C_strat, T_omega_inv2.T])

    if np.sum(np.abs(C_try1[-1, :2])) > np.sum(np.abs(C_try2[-1, :2])):
        # on cherche le repère d orthotropie avec les couplages les moins forts

        C_final = C_try2

        omega = omega2

    else:
        C_final = C_try1

        omega = omega1

    # Clean-up numerical errors

    C_final[-1, :2] = np.zeros(2)

    C_final[:2, -1] = np.zeros(2)

    return np.linalg.inv(C_final), omega

Mat_rot

Mat_rot(angle: float) -> ndarray

Return the rotation matrix for a given angle.

Source code in src\vimseo\lib_vimseo\tan_lib.py

def Mat_rot(angle: float) -> np.ndarray:
    """Return the rotation matrix for a given angle."""

    s = np.sin(angle)

    c = np.cos(angle)

    return np.array([
        [c**2, s**2, 2 * c * s],
        [s**2, c**2, -2 * c * s],
        [-c * s, c * s, c**2 - s**2],
    ])

S_12

S_12(S_strat: ndarray) -> tuple[complex, complex]

Return the two roots of the characteristic equation.

Source code in src\vimseo\lib_vimseo\tan_lib.py

def S_12(S_strat: np.ndarray) -> tuple[complex, complex]:
    """Return the two roots of the characteristic equation."""

    gamma0 = np.sqrt(S_strat[1, 1] / S_strat[0, 0])

    beta0 = (2 * S_strat[0, 1] + S_strat[2, 2]) / (2 * S_strat[0, 0])

    if gamma0 > beta0:
        s1 = np.sqrt((gamma0 - beta0) / 2) + 1j * np.sqrt((gamma0 + beta0) / 2)

        s2 = -np.sqrt((gamma0 - beta0) / 2) + 1j * np.sqrt((gamma0 + beta0) / 2)

        return s1, s2

    s1 = 1j * (np.sqrt((beta0 - gamma0) / 2) + np.sqrt((gamma0 + beta0) / 2))

    s2 = 1j * (-np.sqrt((beta0 - gamma0) / 2) + np.sqrt((gamma0 + beta0) / 2))

    return s1, s2

force_fluxes

force_fluxes(
    R: float,
    s1: complex,
    s2: complex,
    p: ndarray,
    x: float,
    y: float,
) -> ndarray

Return the force fluxes at a given point.

Source code in src\vimseo\lib_vimseo\tan_lib.py

def force_fluxes(
    R: float,
    s1: complex,
    s2: complex,
    p: np.ndarray,
    x: float,
    y: float,
) -> np.ndarray:
    """Return the force fluxes at a given point."""

    Phi_12 = phi(x, y, p, s1, s2, R)  # tuple de (phi1,phi2)

    N_I = np.zeros(3)

    N_II = np.zeros(3)

    N_I[0] = p[0] + 2 * np.real(s1**2 * Phi_12[0][0] + s2**2 * Phi_12[1][0])

    N_I[1] = 2 * np.real(Phi_12[0][0] + Phi_12[1][0])

    N_I[2] = -2 * np.real(s1 * Phi_12[0][0] + s2 * Phi_12[1][0])

    N_II[0] = 2 * np.real(s1**2 * Phi_12[0][1] + s2**2 * Phi_12[1][1])

    N_II[1] = p[1] + 2 * np.real(Phi_12[0][1] + Phi_12[1][1])

    N_II[2] = -2 * np.real(s1 * Phi_12[0][1] + s2 * Phi_12[1][1])

    return np.vstack((N_I, N_II))

phi

phi(
    x: float,
    y: float,
    p: ndarray,
    s1: complex,
    s2: complex,
    R: float,
) -> tuple[ndarray, ndarray]

Return the phi1 and phi2 values for given x, y, p, s1, s2 and R.

Source code in src\vimseo\lib_vimseo\tan_lib.py

def phi(
    x: float, y: float, p: np.ndarray, s1: complex, s2: complex, R: float
) -> tuple[np.ndarray, np.ndarray]:
    """Return the phi1 and phi2 values for given x, y, p, s1, s2 and R."""

    phi1 = np.zeros(2, dtype=complex)

    phi2 = np.zeros(2, dtype=complex)

    z1 = x + s1 * y

    z2 = x + s2 * y

    zeta_ = zeta(z1, z2, s1, s2, R)

    phi1[0] = -1j * p[0] * (1 - zeta_[0]) / (2 * (s1 - s2) * (1 + s1 * 1j))

    phi1[1] = p[1] * s2 * (1 - zeta_[0]) / (2 * (s1 - s2) * (1 + s1 * 1j))

    phi2[0] = 1j * p[0] * (1 - zeta_[1]) / (2 * (s1 - s2) * (1 + s2 * 1j))

    phi2[1] = -p[1] * s1 * (1 - zeta_[1]) / (2 * (s1 - s2) * (1 + s2 * 1j))

    return phi1, phi2

principal_stress

principal_stress(N: ndarray) -> float

Return the angle of the principal stress directions.

Source code in src\vimseo\lib_vimseo\tan_lib.py

def principal_stress(N: np.ndarray) -> float:
    """Return the angle of the principal stress directions."""

    if N[0] == N[1]:
        a1 = np.pi / 4 * (N[2] > 0) - np.pi / 4 * (N[2] < 0)

    else:
        a1 = 1 / 2 * np.arctan(2 * N[2] / (N[0] - N[1]))

    return a1

tan_model

tan_model(
    r: float,
    theta: float,
    N: ndarray,
    C_strat: ndarray,
    R: float,
    w: float,
) -> ndarray

Return the stress at a given point for a given set of parameters.

Parameters:

• r (float) –

  The distance from the hole center to the point where the stress is calculated.

• theta (float) –

  The angle between the x-axis and the line connecting the hole center to the point where the stress is calculated.

• N (ndarray) –

  The applied load vector (N_x, N_y, N_xy).

• C_strat (ndarray) –

  The stiffness matrix of the laminate.

• R (float) –

  The radius of the hole.

• w (float) –

  The width of the plate.

Source code in src\vimseo\lib_vimseo\tan_lib.py

def tan_model(
    r: float, theta: float, N: np.ndarray, C_strat: np.ndarray, R: float, w: float
) -> np.ndarray:
    """Return the stress at a given point for a given set of parameters.

    Args:
        r: The distance from the hole center to the point where the stress is calculated.
        theta: The angle between the x-axis and the line connecting the hole center to the point where the stress is calculated.
        N: The applied load vector (N_x, N_y, N_xy).
        C_strat: The stiffness matrix of the laminate.
        R: The radius of the hole.
        w: The width of the plate.
    """
    # angle of the principal axis
    a = principal_stress(N)  

    M_rot = Mat_rot(a)

    # stress in the principal frame
    p = np.dot(M_rot, N)  

    p[2] = 0

    x = r * np.cos(theta - a)

    y = r * np.sin(theta - a)

    # stiffness in the principal frame
    C_pli = np.linalg.multi_dot([
        M_rot,
        C_strat,
        M_rot.T,
    ])  

    s1, s2 = Calc_S12_eff(C_pli)  

    # Infinite force fluxes along direction 1 and 2
    F_ = force_fluxes(R, s1, s2, p, x, y)

    u = 2 * R / w

    # Correction for a finite plate
    coeff = (2 + (1 - u) ** 3) / (3 - 3 * u)

    return coeff * np.linalg.solve(M_rot, np.sum(F_, axis=0))

tan_model_grid

tan_model_grid(
    r: ndarray,
    theta: ndarray,
    N: ndarray,
    C_strat: ndarray,
    R: float,
    w: float,
) -> ndarray

Vectorised :func:tan_model over arrays of (r, theta).

Everything that does not depend on the evaluation point (principal frame, effective roots s1/s2, finite-width correction) is computed once instead of once per point, and the force fluxes are evaluated with array operations. The result is numerically identical to calling tan_model on each point.

Parameters:

• r (ndarray) –

  Distances from the hole centre, shape (K,).

• theta (ndarray) –

  Angles from the x-axis, shape (K,).

• N (ndarray) –

  The applied load vector (N_x, N_y, N_xy).

• C_strat (ndarray) –

  The stiffness matrix of the laminate.

• R (float) –

  The radius of the hole.

• w (float) –

  The width of the plate.

Returns:

• ndarray –

  The stresses (N_xx, N_yy, N_xy) at each point, shape (K, 3).

Source code in src\vimseo\lib_vimseo\tan_lib.py

def tan_model_grid(
    r: np.ndarray,
    theta: np.ndarray,
    N: np.ndarray,
    C_strat: np.ndarray,
    R: float,
    w: float,
) -> np.ndarray:
    """Vectorised :func:`tan_model` over arrays of ``(r, theta)``.

    Everything that does not depend on the evaluation point (principal frame,
    effective roots ``s1``/``s2``, finite-width correction) is computed once
    instead of once per point, and the force fluxes are evaluated with array
    operations. The result is numerically identical to calling ``tan_model`` on
    each point.

    Args:
        r: Distances from the hole centre, shape ``(K,)``.
        theta: Angles from the x-axis, shape ``(K,)``.
        N: The applied load vector ``(N_x, N_y, N_xy)``.
        C_strat: The stiffness matrix of the laminate.
        R: The radius of the hole.
        w: The width of the plate.

    Returns:
        The stresses ``(N_xx, N_yy, N_xy)`` at each point, shape ``(K, 3)``.
    """
    r = np.asarray(r, dtype=float)

    theta = np.asarray(theta, dtype=float)

    # --- quantities constant over the grid (computed once) ---
    a = principal_stress(N)

    M_rot = Mat_rot(a)

    p = np.dot(M_rot, N)

    p[2] = 0

    C_pli = np.linalg.multi_dot([M_rot, C_strat, M_rot.T])

    s1, s2 = Calc_S12_eff(C_pli)

    u = 2 * R / w

    coeff = (2 + (1 - u) ** 3) / (3 - 3 * u)

    # --- per-point force fluxes (vectorised force_fluxes + phi) ---
    x = r * np.cos(theta - a)

    y = r * np.sin(theta - a)

    z1 = x + s1 * y

    z2 = x + s2 * y

    zeta1, zeta2 = zeta(z1, z2, s1, s2, R)

    denom1 = 2 * (s1 - s2) * (1 + s1 * 1j)

    denom2 = 2 * (s1 - s2) * (1 + s2 * 1j)

    # phi1 uses zeta1 / denom1, phi2 uses zeta2 / denom2 (see ``phi``).
    phi1_0 = -1j * p[0] * (1 - zeta1) / denom1

    phi1_1 = p[1] * s2 * (1 - zeta1) / denom1

    phi2_0 = 1j * p[0] * (1 - zeta2) / denom2

    phi2_1 = -p[1] * s1 * (1 - zeta2) / denom2

    # N_I (column 0) + N_II (column 1), matching ``force_fluxes`` then summing.
    flux = np.empty((3, r.size))

    flux[0] = (
        p[0]
        + 2 * np.real(s1**2 * phi1_0 + s2**2 * phi2_0)
        + 2 * np.real(s1**2 * phi1_1 + s2**2 * phi2_1)
    )

    flux[1] = (
        2 * np.real(phi1_0 + phi2_0)
        + p[1]
        + 2 * np.real(phi1_1 + phi2_1)
    )

    flux[2] = (
        -2 * np.real(s1 * phi1_0 + s2 * phi2_0)
        - 2 * np.real(s1 * phi1_1 + s2 * phi2_1)
    )

    # Back to the global frame; solve broadcasts over the columns of ``flux``.
    return (coeff * np.linalg.solve(M_rot, flux)).T

zeta

zeta(
    z1: complex,
    z2: complex,
    s1: complex,
    s2: complex,
    R: float,
) -> tuple[complex, complex]

Return the zeta values for given z1, z2, s1, s2 and R.

Source code in src\vimseo\lib_vimseo\tan_lib.py

def zeta(
    z1: complex, z2: complex, s1: complex, s2: complex, R: float
) -> tuple[complex, complex]:
    """Return the zeta values for given z1, z2, s1, s2 and R."""

    val1 = z1 / np.sqrt(z1**2 - R**2 * (1 + s1**2))

    val2 = z2 / np.sqrt(z2**2 - R**2 * (1 + s2**2))

    # Select the root with non-negative real part, independently for each value.
    # Written with ``np.where`` so it works for both scalars and grid arrays; the
    # original four if/else cases were exactly this per-value sign flip.
    val1 = np.where(np.real(val1) >= 0, val1, -val1)

    val2 = np.where(np.real(val2) >= 0, val2, -val2)

    return val1, val2