coefficients.py

Go to the documentation of this file.

"""Copyright (C) 2016 Joakim A. Taby

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>."""
import numpy as np


def analytical_contact_pressures(sigma,alpha,n,A,R,t,p_int,p_ext):
    """ Finding contact pressure from pressure differencial between layers. Using B1.23 from handbook, assuming
    sliding occurs along midpoint between layers.
    """

    out_in=True

    if len(sigma.shape) ==1:
        h_length=1
    else:
        h_length=len(sigma[:,0])

    R_inner=np.zeros(len(A))
    R_outer=np.zeros(len(A))
    for i in range(len(A)):
        if i ==0:
            R_inner[i]=R[i]-t[i]/2.0
            R_outer[i]=(R[i+1]-t[i+1]/2.0+R[i]+t[i]/2.0)/2.0
        elif i== len(A)-1:
            R_inner[i]=(R[i]-t[i]/2.0+R[i-1]+t[i-1]/2.0)/2.0
            R_outer[i]=R[i]+t[i]/2.0
        else:
            R_inner[i]=(R[i]-t[i]/2.0+R[i-1]+t[i-1]/2.0)/2.0
            R_outer[i]=(R[i+1]-t[i+1]/2.0+R[i]+t[i]/2.0)/2.0

    pressure_diff_out_in = n*sigma*A*np.sin(alpha)**2/(2.0*np.pi*np.cos(alpha)*R*R_inner)

    if out_in:

        pressure_contact_out_in=np.zeros((len(A),h_length))



        for i in range(len(A)-1,-1,-1):
            if i == len(A)-1:
                    pressure_contact_out_in[i]=p_ext*(R[i]+t[i]/2.0)/R_inner[i]+pressure_diff_out_in.T[i]
            else:
                pressure_contact_out_in[i]=pressure_diff_out_in.T[i]+pressure_contact_out_in[i+1]*R_outer[i]/R_inner[i]

        return pressure_contact_out_in.T

    else:

        pressure_diff_in_out=n*sigma*A*np.sin(alpha)**2/(2.0*np.pi*np.cos(alpha)*R*R_outer)
        pressure_contact_in_out=np.zeros((len(A),h_length))
        for i in range(len(A)):
            if i == 0:
                    pressure_contact_in_out[i]=p_int
            else:
                pressure_contact_in_out[i]=(pressure_contact_in_out[i-1]*R_inner[i-1]/R_outer[i-1]
                                           -pressure_diff_in_out.T[i-1])

        return pressure_contact_in_out.T


def critical_curvature(slenderobject, contact_pressure, frictionfactors=np.array([5, 5, 5])):
    """ Finding the critical curvature. Need to find where this is taken from.
    :param frictionfactors: friction factors between layer and the layer inside
    :param contact_pressure: contact line load. A matrix [2,nlayers] for outside and inside of layer
                          respectively
    :param slenderobject: Object of slender type from the sledner class of slenders
    """
    # checking data input
    if isinstance(frictionfactors, np.int) or isinstance(frictionfactors, np.float):
        fricfac = []
        for i in range(len(slenderobject.thickness)):
            fricfac.append(frictionfactors)
    elif max(abs(frictionfactors)) == 5:
        fricfac = slenderobject.fricfac
    elif len(frictionfactors) == 1:
        fricfac = []
        for i in range(len(slenderobject.thickness)):
            fricfac.append(frictionfactors)
    elif len(frictionfactors) == len(slenderobject.fricfac):
        fricfac = frictionfactors
    else:
        print("Cannot parse the friction factors given. Needs to be one single friction "
              "factor, or the same number of friction factors as defined layers in "
              "the slender object")

    tmp = contact_pressure.shape
    dim1len=tmp[0]
    dim2len = len(slenderobject.thickness)
    critcurv = np.zeros((dim2len,dim1len))
    for i in range(len(slenderobject.layer_radii)):
        outerfillfac=slenderobject.comp_number[i]*slenderobject.width[i]/\
                     (2*np.pi*(slenderobject.layer_radii[i]-slenderobject.thickness[i])*
                      np.cos(slenderobject.lay_angle[i]))
        innerfillfac=slenderobject.comp_number[i]*slenderobject.width[i]/\
                     (2*np.pi*(slenderobject.layer_radii[i]+slenderobject.thickness[i])*
                      np.cos(slenderobject.lay_angle[i]))
        if not i == len(slenderobject.fricfac)-1:
            if slenderobject.lay_angle[i] == 0:
                critcurv[i] = 0.0
            else:
                critcurv[i] = ((fricfac[i] * contact_pressure.T[i]/outerfillfac + fricfac[i + 1]/innerfillfac * contact_pressure.T[i + 1]) /
                               (slenderobject.youngs_modulus[i] * slenderobject.thickness[i] *
                                   np.cos(slenderobject.lay_angle[i]) ** 2.0 *
                                   abs(np.sin(slenderobject.lay_angle[i]))))
        else:
            if slenderobject.lay_angle[i] == 0:
                critcurv[i] = 0.0
            else:
                critcurv[i] = (fricfac[i] * contact_pressure.T[i] /
                               (slenderobject.youngs_modulus[i] * slenderobject.thickness[i] *
                                   np.cos(slenderobject.lay_angle[i]) ** 2.0 *
                                   abs(np.sin(slenderobject.lay_angle[i]))))

    return critcurv.T


# @profile
def local_bending(thickness,width, alpha, theta, phi, e,r, type,bend_model="Savik"):
    """ Calculating the local bending coefficients
    Equations from 'On stresses and fatigue in flexible pipes' by Savik 1992 and
    'Validity and limitation of analytical models for the bending stress of a helical
    wire in unbonded flexible pipes' Tang et.al. 2015

    :param bend_model: Which bend model to apply. Savik or Tan implemented
    :param e: Youngs modulus. One for each layer
    :param geometry of the tendon. one item means circular, two means rectangular
    :param alpha: lay_angle of slender structure, list
    :param theta: position around the crossection for calculation, numpy vector
    :param phi: local position for calculation of coefficient, numpy vector
    :return:list of numpy arrays defining the coefficients for each layer
    """
    # reshaping the arrays to given input
    phi = np.reshape(phi, (1, -1))
    bendcoeff = []
    # Choosing bending model
    for i in range(len(thickness)):
        if type[i]=="wires":
            if bend_model.upper() == "SAVIK":
                weak_curvoncurv = -(np.cos(alpha[i]) ** 4) * np.cos(theta)

                strong_curvoncurv = (1 + (np.sin(alpha[i])) ** 2) * np.cos(alpha[i]) * np.sin(theta)
            elif bend_model.upper() == "TAN":
                strong_curvoncurv = np.sin(theta)
                weak_curvoncurv = -np.cos(theta) * np.cos(alpha[i])

            tmp = np.zeros((len(theta[:, 0]), len(theta[0, :]), 4))
            tmp[:, :, 0] = np.squeeze(-weak_curvoncurv * thickness[i] * 0.5 - strong_curvoncurv * width[i] * 0.5)
            tmp[:, :, 1] = np.squeeze(
                weak_curvoncurv * thickness[i] * 0.5 - strong_curvoncurv * width[i] * 0.5)
            tmp[:, :, 2] = np.squeeze(
                +weak_curvoncurv * thickness[i] * 0.5 + strong_curvoncurv * width[i] * 0.5)
            tmp[:, :, 3] = np.squeeze(-weak_curvoncurv * thickness[i] * 0.5 + strong_curvoncurv * width[i] * 0.5)
            bendcoeff.append(tmp * e[i])
        elif type[i]=="sheath":
            bendcoeff.append(e[i]*(r[i]+thickness[i]/2.0)*np.cos(theta))

    return bendcoeff


def tension_stress_history(t_eff, p_int,p_ext,alpha,R,thick,A,E,n):
    """
    Solving for axial stress accorind got Handbook B1.10, B1.23 and B1.13

    :param t_eff: history of effective tension
    :param p_int: history of internal pressure
    :param p_ext: history of external pressure
    :param alpha: lay angles of the layers
    :param R: radiuses of the layers
    :param thick: thickness of the layers
    :param A: Area of the tendons
    :param E: young's modulus of the tendons
    :param n: number of tendons in each layer
    :return: numpy array with tension stress histories
    """
    h_length=len(t_eff)
    r_int=R[0]-thick[0]/2.0
    r_ext=R[-1]+thick[-1]/2.0
    Ka = t_eff+np.pi*p_int*r_int**2-np.pi*p_ext*r_ext**2
    Kp = 2*np.pi*(p_int*r_int-p_ext*r_ext)
    K1 = 0.0
    K2 = 0.0
    K3 = 0.0
    for i in range(len(A)):
        K1+=n[i]*A[i]*E[i]*np.sin(alpha[i])**2*np.cos(alpha[i])/R[i]
        K2+=n[i]*A[i]*E[i]*np.sin(alpha[i])**4/(R[i]**2*np.cos(alpha[i]))
        K3+=n[i]*A[i]*E[i]*np.cos(alpha[i])**3

    u3 = (Ka*K1-K3*Kp)/(K1**2-K3*K2)
    eps_p = (Ka-K1*u3)/K3

    sigma=np.zeros((i+1,h_length))
    for i in range(len(A)):
        sigma[i]=eps_p*E[i]*np.cos(alpha[i])**2+u3*E[i]*np.sin(alpha[i])**2/R[i]

    return sigma.T


def bending_stress_history(slenderobject, curvhist, ntheta=16, nphi=16, bend_model="Savik"):
    """
    Calculates the stress history from time histories of curvatures and axial stress coefficient
    :param bend_model: Which bend model to apply. Savik or Tan implemented
    :param nphi: Number of points around tendon to perform calculation. 4 is used of rectangular
                 tendon given.
    :param ntheta: number of points around the cross section to perform calculation
    :param slenderobject, instande of the Slender class
    :param curvhist: numpy array containg each curvature time series as a column vector
    :return: stress history for each point on defined in bendcoeff
    """
    phi = np.arange(0, 2 * np.pi, 2 * np.pi / nphi)
    # decompose curvature in order to find the relative theta
    theta_load = np.arctan2(curvhist[:, 1], curvhist[:, 0])
    theta_load = theta_load.reshape(-1, 1)
    abscurv = np.sqrt(curvhist[:, 0] ** 2 + curvhist[:, 1] ** 2)
    # augment thetas for call and for call to bendcoeff
    theta0 = np.arange(0, 2 * np.pi, 2 * np.pi / ntheta)
    bendstress = []
    tmp = np.repeat(theta0, len(theta_load))
    theta0_expanded = tmp.reshape((len(theta0), len(theta_load)))
    theta = np.subtract(theta0_expanded.T, theta_load)
    del theta0
    del theta0_expanded
    del tmp
    del theta_load
    del curvhist
    bendcoeff = local_bending(slenderobject.thickness,slenderobject.width, slenderobject.lay_angle,
                              theta, phi, slenderobject.youngs_modulus,
                              slenderobject.layer_radii,slenderobject.layer_type,
                               bend_model=bend_model)
    for j in range(len(slenderobject.thickness)):
        bendstress.append(np.multiply(bendcoeff[j].T, abscurv).T)
    return bendstress


def stick_stress(model, curv, critcurv, theta=None):
    """Calculating the stick stress at position theta on the crossection, limited by the critical curvature.
    Equations from 'Bending Behavior of Helically Wrapped Cables' Hong et.al. 2005
    :param theta: Position around cross section range [0 2*pi]
    :param critcurv: Critical curvature. Curvature at which slip initiates for each layer
    :param curv: Curvature history Dim: 2xN
    :param model: Slender object as defined in slenders.Slender"""
    curv_norm = np.sqrt(curv[:, 0] ** 2 + curv[:, 1] ** 2)
    scaler = []
    scaled_curv = []
    count = 0
    if not theta:
        theta = []
        # find angles from standarddev if we have more than one length
        # if not, we calculate directly
        for i in range(len(model.layer_radii)):
            if len(curv[:, 0]) > 1:
                y = np.std(curv[:, 1])  # +np.pi/2.0
                x = np.std(curv[:, 0])  # +np.pi/2.0
            else:
                y = curv[:, 1]  # +np.pi/2.0
                x = curv[:, 0]  # +np.pi/2.0
            tmp = np.arctan2(y, x)
            if y == 0:
                tmp = np.pi / 2.0
            if x == 0:
                tmp = 0
            if tmp < 0:
                theta.append(tmp + np.pi / 2)
            else:
                theta.append(tmp - np.pi / 2)

    for i in np.transpose(critcurv):
        tmp = np.clip(i / np.clip(curv_norm,0.0000000001,float("inf")), 0, 1)
        scaled_curv.append(curv[:, 0] * tmp * np.sin(theta[count]) - curv[:, 1] * tmp * np.cos(theta[count]))
        count += 1
    stickstress = []
    for i in range(len(model.layer_radii)):
        outer_radius = model.layer_radii[i] + model.width[i]/2.0
        stickstress.append(model.youngs_modulus[i] * np.cos(model.lay_angle[i]) ** 2 * outer_radius * scaled_curv[i])

    return stickstress


def fric_stress(curv, eps_mat, sig_mat, midpoint):
    """
    2D elastoplastic materialmodel with kinematic hardening
    :param midpoint: the midpoint of the yield surface at simulation start
    :param sig_mat: the stress at yield
    :param eps_mat: the strain at yield
    :param curv: curvature history 2xN
    :return: Frictional stress
    """
    # it should take in 2 curvatures, the midpoint and the material model
    # and should return the moment or stress or whatever is the quantity of interest
    # Possibly also the corresponding stiffness
    dimension = len(curv[0, :])
    # taking advantage of that the stick stiffness doesnt vary with time
    # calculating stiffnesses for each elastic-perfectly plastic region
    stiffness = sig_mat[0] / eps_mat[0]
    stresshistory = []
    # Starting timeloop
    for i in range(len(curv[:, 0])):
        stress = 0.0
        # calculate radius from mid
        rad = 0.0
        for j in range(dimension):
            rad += (curv[i, j] - midpoint[j]) ** 2
        rad = np.sqrt(rad)
        diffcurv = curv[i, :] - midpoint
        # finding the distance that the midpoint must be moved
        if rad > 0.0:
            movefactor = ((rad - eps_mat[i]) / rad)
        else:
            movefactor = 0.0
        if movefactor > 0.0:
            # Moving midpoint in direction current curvature relative to mid point
            midpoint += movefactor * diffcurv
            diffcurv = curv[i, :] - midpoint
        # Calculating stress contributions
        tmpstress = stiffness * diffcurv.T
        stress += tmpstress
        stresshistory.append(stress)
    stresshistory = np.array(stresshistory)
    return stresshistory, midpoint


def max_stickstress(model, critcurv):
    max_stick_stress = []
    for i in range(len(model.layer_radii)):
        outer_radius = model.layer_radii[i] + model.thickness[i] / 2.0
        max_stick_stress.append(model.youngs_modulus[i] * np.cos(model.lay_angle[i]) ** 2 * outer_radius * critcurv[i])

    return max_stick_stress