axisymmetric.py

Go to the documentation of this file.

import numpy as np
import collections
import itertools



def solver(layerinfo,r,t,b,A,alpha,n,E,nu,gap_ini,intpresslayer,extpresslayer
           , T_eff, p_int, p_ext):

    # calculation section
    # setting up the equation system
    K,R,eps11index,u3index,deltatindex,pindex,fillfactor_in,fillfactor_out = systemsetup(r,t,b,A,n,E,nu,alpha,gap_ini,intpresslayer,extpresslayer,layerinfo,T_eff,p_int,p_ext)
    solution=np.linalg.solve(K,R)

    np.savetxt('k.csv', K, fmt='%.18e', delimiter=',', newline='\n')
    np.savetxt('solution.csv', solution[:,0], fmt='%.18e', delimiter=',', newline='\n')
    np.savetxt('r.csv', R[:,0], fmt='%.18e', delimiter=',', newline='\n')

    Nlayers=len(layerinfo)
    # correcting the solution if there should be gaps
    solution=corrector(r,t,b,A,n,E,nu,alpha,gap_ini,intpresslayer,extpresslayer,layerinfo,T_eff,p_int,p_ext,solution,pindex,Nlayers)

    # finding the stresses
    stress=np.zeros((Nlayers,len(T_eff),3))
    for i in range(Nlayers):
        if layerinfo[i]=='wires':
            stress[i,:,0]=solution[i,:]*E[i]
            if i>0:
                stress[i,:,0]+=nu[i]*solution[pindex+i,:]/(2.0*fillfactor_in[i])
            if i<Nlayers-1:
                stress[i,:,0]+=nu[i]*solution[pindex+1+i,:]/(2.0*fillfactor_out[i])
            if i== intpresslayer:
                stress[i,:,0]+=nu[i]*p_int/(2.0*fillfactor_in[i])
            if i== extpresslayer:
                stress[i,:,0]+=nu[i]*p_ext/(2.0*fillfactor_out[i])
            stress[i,:,1]=0.0
            stress[i,:,2]=solution[deltatindex+i,:]/t[i]*E[i]-nu[i]*stress[i,:,0]

        elif layerinfo[i]=='sheath':
            router=r[i]+t[i]/2.0
            rinner=r[i]-t[i]/2.0
            #Finding coefficients from standard thick shell solution (https://en.wikipedia.org/wiki/Cylinder_stress)
            if i>0:
                stress[i,:,1]+=(solution[pindex+i+1,:]*(r[i]+t[i]/2.0)**2*(r[i]-t[i]/2.0)**2
                                / (((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2)*r[i]**2) +
                                solution[pindex+i,:]*(r[i]-t[i]/2.0)**2
                                /((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2))
                stress[i,:,2]+=(-solution[pindex+i,:]*(r[i]+t[i]/2.0)**2*(r[i]-t[i]/2.0)**2
                                / (((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2)*r[i]**2) +
                                solution[pindex+i,:]*(r[i]-t[i]/2.0)**2
                                /((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2))

            if i< Nlayers-1:
                stress[i,:,1]+=(-solution[pindex+i+1,:]*(r[i]+t[i]/2.0)**2*(r[i]-t[i]/2.0)**2
                                / (((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2)*r[i]**2) -
                                solution[pindex+i,:]*(r[i]+t[i]/2.0)**2
                                /((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2))
                stress[i,:,2]+=(solution[pindex+i+1,:]*(r[i]+t[i]/2.0)**2*(r[i]-t[i]/2.0)**2
                                / (((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2)*r[i]**2) -
                                solution[pindex+i,:]*(r[i]+t[i]/2.0)**2
                                /((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2))
            if i==intpresslayer:
                stress[i,:,1]+=(p_int*(r[i]+t[i]/2.0)**2*(r[i]-t[i]/2.0)**2
                                / (((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2)*r[i]**2) +
                                p_int*(r[i]-t[i]/2.0)**2
                                /((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2))
                stress[i,:,2]+=(-p_int*(r[i]+t[i]/2.0)**2*(r[i]-t[i]/2.0)**2
                                / (((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2)*r[i]**2) +
                                p_int*(r[i]-t[i]/2.0)**2
                                /((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2))
            if i== extpresslayer:
                stress[i,:,1]+=(-p_ext*(r[i]+t[i]/2.0)**2*(r[i]-t[i]/2.0)**2
                                / (((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2)*r[i]**2) -
                                p_ext*(r[i]+t[i]/2.0)**2
                                /((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2))
                stress[i,:,2]+=(p_ext*(r[i]+t[i]/2.0)**2*(r[i]-t[i]/2.0)**2
                                / (((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2)*r[i]**2) -
                                p_ext*(r[i]+t[i]/2.0)**2
                                /((r[i]+t[i]/2.0)**2-(r[i]-t[i]/2.0)**2))

            stress[i,:,0]+=E[i]*solution[-1,:]+nu[i]*(stress[i,:,1]+stress[i,:,2])

    contpress=np.zeros((len(T_eff),Nlayers))
    contpress[:,1:]=solution[pindex+1:pindex+Nlayers].T
    axialstress = stress[:,:,0].T
    hoopstress  = stress[:,:,1].T
    radialstress= stress[:,:,2].T
    # printout for testing
    np.savetxt('stress.csv', stress[:,:,0], fmt='%.18e', delimiter=',', newline='\n')
    np.savetxt('contpress.csv', contpress.T, fmt='%.18e', delimiter=',', newline='\n')
    np.savetxt('finalsolution.csv', solution[:,-1], fmt='%.18e', delimiter=',', newline='\n')
    return axialstress,hoopstress,radialstress,contpress


def combfinder(combinations):
    sol=[]
    for i in range(1,len(combinations)+1):
        sol.append(itertools.combinations(combinations,i))
    return sol

def systemsetup(r,t,b,A,n,E,nu,alpha,gap_ini,intpresslayer,extpresslayer,layerinfo,T_eff,p_int,p_ext,gaplayers=[0]):
    # finding constants and setting up matrices and counters
    fillfactor_in=n*b/(2*np.pi*(r-t)*np.cos(alpha))
    fillfactor_out=n*b/(2*np.pi*(r+t)*np.cos(alpha))
    Nlayers=len(E)

    linecounter=0
    presscounter=1
    temp=collections.Counter(layerinfo)
    Nhelic=temp['wires']
    Nsheath=temp['sheath']
    eps11index=0
    u3index=Nhelic
    deltatindex=Nhelic+Nlayers
    pindex=Nhelic+2*Nlayers-1
    ihelic=-1

    if gaplayers==[0]:
        layersequence=[(0,Nlayers)]
        nequations=Nhelic*4+3*Nsheath
    else:
        tmp=list(gaplayers)
        nequations=Nhelic*4+3*Nsheath-len(gaplayers)
        layersequence=[]
        layersequence.append((0,tmp[0]+1))
        for i in range(len(gaplayers)):
            if i<len(gaplayers)-1:
                layersequence.append((tmp[i]+1,tmp[i+1]+1))
            else:
                layersequence.append((tmp[i]+1,Nlayers))
    K=np.zeros((nequations,nequations))
    R=np.zeros((nequations,len(T_eff)))
    # setting up equations system
    for currentlayers in layersequence:
        presscounter-=1

        for i in range(currentlayers[0],currentlayers[1]):
            if layerinfo[i]=='wires':
                ihelic+=1
                # solving for
                #local strain, radial expansion, change in thickness, contact pressure,global strain
                #eps1,eps2...,u31,u32,u33...,deltat1,deltat2,deltat3...,p1,p2,p3,..epsglobal
                #  kinematics local/global

                K[linecounter,eps11index+ihelic]+=-r[i]
                K[linecounter,u3index+i]+=np.sin(alpha[i])**2
                K[linecounter,-1]+=np.cos(alpha[i])**2*r[i]
                #Last row is total axial equilibrium
                K[-1,eps11index+ihelic]+=n[i]*E[i]*A[i]*np.cos(alpha[i])
                if i>currentlayers[0]:
                    K[-1,pindex+presscounter]+=-nu[i]*n[i]*A[i]*np.cos(alpha[i])/(2.0*fillfactor_in[i])
                if i<currentlayers[1]-1:
                    K[-1,pindex+presscounter+1]  +=-nu[i]*n[i]*A[i]*np.cos(alpha[i])/(2.0*fillfactor_out[i])
                if i==0:
                    R[-1,:]+=T_eff+np.pi*p_int*(r[intpresslayer]-t[intpresslayer]/2.0)**2-np.pi*p_ext*(r[extpresslayer]+t[extpresslayer]/2.0)**2
                if i==intpresslayer:
                    R[-1,:]+=nu[i]*n[i]*A[i]*np.cos(alpha[i])/(2.0*fillfactor_in[i])*p_int
                if i==extpresslayer:
                    R[-1,:]+=nu[i]*n[i]*A[i]*np.cos(alpha[i])/(2.0*fillfactor_out[i])*p_ext
                # Second row is pressure equilibrium
                # check row number
                linecounter+=1
                K[linecounter,eps11index+ihelic]+=np.sin(alpha[i])**2*A[i]*E[i]*n[i]/(np.cos(alpha[i])*r[i])
                if i>currentlayers[0]:
                    K[linecounter,pindex+presscounter]+=-nu[i]*np.sin(alpha[i])**2*A[i]*n[i]/(fillfactor_in[i]*2*r[i]*np.cos(alpha[i]))-\
                                   2*np.pi*(r[i]-t[i]/2.0)
                if i<currentlayers[1]-1:
                    K[linecounter,pindex+presscounter+1]  +=-nu[i]*np.sin(alpha[i])**2*A[i]*n[i]/(fillfactor_out[i]*2*r[i]*np.cos(alpha[i]))+\
                                   2*np.pi*(r[i]+t[i]/2.0)
                if i == intpresslayer:
                    R[linecounter]+=+nu[i]*np.sin(alpha[i])**2*A[i]*n[i]/(fillfactor_in[i]*2*r[i]*\
                                np.cos(alpha[i]))*p_int+\
                                   2*np.pi*(r[i]-t[i]/2.0)*p_int
                if i == extpresslayer:
                    R[linecounter]+=+nu[i]*np.sin(alpha[i])**2*A[i]*n[i]/(fillfactor_out[i]*2*r[i]*\
                                np.cos(alpha[i]))*p_ext\
                                  -2*np.pi*(r[i]+t[i]/2.0)*p_ext
                # Radial squeeze
                linecounter+=1
                K[linecounter,eps11index+ihelic]+=nu[i]
                K[linecounter,deltatindex+i]+=1.0/t[i]
                if i>currentlayers[0]:
                    K[linecounter,pindex+presscounter]+=(1+nu[i]**2)/(2*E[i]*fillfactor_in[i])
                if i<currentlayers[1]-1:
                    K[linecounter,pindex+presscounter+1]  +=(1+nu[i]**2)/(2*E[i]*fillfactor_out[i])
                if i==intpresslayer:
                    R[linecounter]               +=-(1+nu[i]**2)/(2*E[i]*fillfactor_in[i])*p_int
                if i == extpresslayer:
                    R[linecounter]               +=-(1+nu[i]**2)/(2*E[i]*fillfactor_out[i])*p_ext

                #Radial continuity
                #Only for closed interfaces of this one
                linecounter+=1
                if i>currentlayers[0]:
                    K[linecounter,u3index+i-1]+=-1
                    K[linecounter,u3index+i]+=1
                    K[linecounter,deltatindex+i-1]+=-0.5
                    K[linecounter,deltatindex+i]+=-0.5
                    R[linecounter]+=-gap_ini[i]
                    linecounter+=1
            # add system for plastic layers
            if layerinfo[i]=='sheath':
                router=r[i]+t[i]/2.0
                rinner=r[i]-t[i]/2.0
                # adding to axial equilibrium
                K[-1,-1]+=A[i]*E[i]
                # using eqs 4 and 5 from Custodio
                if i>currentlayers[0]:
                    K[-1,pindex+presscounter]+=   nu[i]*2.0*rinner**2*A[i]/(router**2-rinner**2)
                if i<currentlayers[1]-1:
                    K[-1,pindex+presscounter+1]+=-nu[i]*2.0*router**2*A[i]/(router**2-rinner**2)
                if i==0:
                    R[-1,:]+=T_eff+np.pi*p_int*(r[intpresslayer]-t[intpresslayer]/2.0)**2-np.pi*p_ext*(r[extpresslayer]+t[extpresslayer]/2.0)**2
                if i ==intpresslayer:
                    R[-1,:]+=                    -nu[i]*2.0*rinner**2*A[i]/(router**2-rinner**2)*p_int
                if i==extpresslayer:
                    R[-1,:]+=                     nu[i]*2.0*router**2*A[i]/(router**2-rinner**2)*p_ext



                # radial equilibrium
                K[linecounter,-1]+=r[i]*nu[i]
                K[linecounter,u3index+i]+=1.0
                if i>currentlayers[0]:
                    K[linecounter,pindex+presscounter]  +=-(1+nu[i])*rinner**2*((2*nu[i]-1)*r[i]**2-router**2)/(E[i]*(rinner**2-router**2)*r[i])
                if i<currentlayers[1]-1:
                    K[linecounter,pindex+presscounter+1]+= (1+nu[i])*router**2*((2*nu[i]-1)*r[i]**2-rinner**2)/(E[i]*(rinner**2-router**2)*r[i])
                if i==intpresslayer:
                    R[linecounter]           +=            (1+nu[i])*rinner**2*((2*nu[i]-1)*r[i]**2-router**2)/(E[i]*(rinner**2-router**2)*r[i])*p_int
                if i==extpresslayer:
                    R[linecounter]           +=           -(1+nu[i])*router**2*((2*nu[i]-1)*r[i]**2-rinner**2)/(E[i]*(rinner**2-router**2)*r[i])*p_ext
                linecounter+=1

                # Thickness equilibrium
                K[linecounter,-1]+=t[i]*nu[i]
                K[linecounter,deltatindex+i]+=1.0
                if i>currentlayers[0]:
                    K[linecounter,presscounter+pindex]+=    (1+nu[i])*rinner*((2*nu[i]-1)*rinner+router)/(E[i]*(rinner+router))
                if i < currentlayers[1]-1:
                    K[linecounter,presscounter+pindex+1]+= -(1+nu[i])*router*((2*nu[i]-1)*router+rinner)/(E[i]*(rinner+router))
                if i==intpresslayer:
                    R[linecounter]             +=          -(1+nu[i])*rinner*((2*nu[i]-1)*rinner+router)/(E[i]*(rinner+router))*p_int
                if i==extpresslayer:
                    R[linecounter]             +=           (1+nu[i])*router*((2*nu[i]-1)*router+rinner)/(E[i]*(rinner+router))*p_ext
                linecounter+=1


                # Radial continuity
                # Only for closed interfaces (could be moved out of loop
                if i>currentlayers[0]:
                    K[linecounter,u3index+i-1]+=-1
                    K[linecounter,u3index+i]+=1
                    K[linecounter,deltatindex+i-1]+=-0.5
                    K[linecounter,deltatindex+i]+=-0.5
                    R[linecounter]+=-gap_ini[i]
                    linecounter+=1
            presscounter+=1
    return K,R,eps11index,u3index,deltatindex,pindex,fillfactor_in,fillfactor_out

def corrector(r,t,b,A,n,E,nu,alpha,gap_ini,intpresslayer,extpresslayer,layerinfo,T_eff,p_int,p_ext,solution,pindex,Nlayers):
    # find solutions with negative gap
    negpresindices=np.where(solution[pindex+1:pindex+Nlayers,:] < 0)

    if len(negpresindices[0])>0:
        # found negative pressure
        negtime=np.unique(negpresindices[1])
        #neglayers=np.unique(negpresindices[0])
        neglayers=list(range(0,Nlayers-1))
        # correct solutions is gap in one or more of the neglayers
        solution_options=combfinder(neglayers)
        # brute force through solving all possible solutions
        for sol in reversed(solution_options):
            for sol2 in sol:
                # fill in changed stuff and solve
                K,R,eps11index,u3index,deltatindex,pindex,fillfactor_in,fillfactor_out =\
                systemsetup(r,t,b,A,n,E,nu,alpha,gap_ini,intpresslayer,extpresslayer,layerinfo,T_eff[negtime],
                            p_int[negtime],p_ext[negtime],gaplayers=sol2)

                solutionalternative=np.linalg.solve(K,R)
                np.savetxt('k1.csv', K, fmt='%.18e', delimiter=',', newline='\n')
                np.savetxt('solution1.csv', solutionalternative[:,0], fmt='%.18e', delimiter=',', newline='\n')
                np.savetxt('r1.csv', R[:,0], fmt='%.18e', delimiter=',', newline='\n')
                # checking if solution is correct at some time, if so, fill into original solution
                negpresindicesalternative=np.where(solutionalternative[pindex+1:pindex+Nlayers-len(sol2),:] < 0)
                negtimealternative=np.unique(negpresindicesalternative[1])
                # check if we have times with no negative pressures
                if len(negtimealternative) < len(negtime):
                    nonegative=sorted(set(range(0,negtime[-1]+1)).difference(negtimealternative))
                    # check if no unallowed layer contact
                    interface=list(sol2)
                    u3indices=np.array([x+u3index for x in interface])
                    u3plusoneindices=np.array([x+u3index+1 for x in interface])
                    deltatindices=np.array([x+deltatindex for x in interface])
                    deltatplusoneindices=np.array([x+deltatindex+1 for x in interface])
                    gapindices=np.array([x+1 for x in interface])
                    # Checking for positive gap
                    #now approved indices is the indices of the new solution, not the original solution
                    approvedindices=np.where(-solutionalternative[u3indices,:]
                                                 -0.5*solutionalternative[deltatindices,:]
                                                 -0.5*solutionalternative[deltatplusoneindices,:]
                                                 +solutionalternative[u3plusoneindices,:]
                                                 +gap_ini[gapindices].reshape(-1,1) > np.min((solutionalternative[deltatindices,:],solutionalternative[deltatplusoneindices,:]))*1.0e-3
                                                 )
                    # check the approved indices for each time against the number that should be approved. Move that into the solution from before, and check if we are done.
                    # check if approved for all interfaces
                    approvedtimeindices=[k for k,v in collections.Counter(approvedindices[1]).items() if v==len(interface)]
                    # Putting the approved solutions into the original solution matrix
                    solution[eps11index:pindex+1,negtime[approvedtimeindices]]=solutionalternative[eps11index:pindex+1,approvedtimeindices]
                    count=0
                    for i in range(Nlayers-1):
                        if i in interface:
                            solution[pindex+1+i,negtime[approvedtimeindices]]=0.0
                        else:
                            solution[pindex+1+i,negtime[approvedtimeindices]]=solutionalternative[pindex+1+count,approvedtimeindices]
                            count+=1
                    solution[-1,negtime[approvedtimeindices]]=solutionalternative[-1,approvedtimeindices]
                    # update negative indices
                    negpresindices=np.where(solution[pindex+1:pindex+Nlayers,:] < 0)
                    if negpresindices[0].size == 0:
                        #solution is great eject!!
                        return solution

    else:
        return solution