Carbon fiber H2 tank optimization problem

tank_opt.py

import numpy as np
from matplotlib import pyplot as plt
import openmdao.api as om

def rotations(theta):
    # define rotation matrices per http://web.mit.edu/course/3/3.11/www/modules/laminates.pdf
    s = np.sin(theta)
    c = np.cos(theta)

    # stress in fiber axis = A * stress in xy axis
    A = np.array([[c**2, s**2, 2*s*c],[s**2, c**2, -2*s*c],[-s*c, s*c, c**2-s**2]])
    Ainv = np.array([[c**2, s**2, -2*s*c],[s**2, c**2, 2*s*c],[s*c, -s*c, c**2-s**2]])
    R = np.diag([1, 1, 2])
    Rinv = np.diag([1, 1, 1/2])
    return A, Ainv, R, Rinv

def analyze_tank_structure(theta_design, t_overall, length, radius, pressure_design=70):
    """
    Computes stresses and failure criterion in a composite cylindrical pressure tank
    Inputs
    ------
    theta_design : float
        Composite fiber angle (in rad) relative to hoop direction.
        Angle relative to axial direction is pi/2 - theta_design
    t_overall : float
        Thickness (in m) of the overall composite laminate (all plies)
    length : float
        Length (in m) of the cylindrical portion of the tank (not including the spherical endcaps)
    radius : float
        Radius (in m) of the spherical endcap and cylindrical tank.
    pressure_design : float
        Design pressure rating (in MPa). 2.35 safety factor applied. 
    
    Outputs
    -------
    TW_failure : float
        Vector of Tsai-Wu failure criterion in each ply (vector, 4 plies)
    weight : float
        Tank weight (in kg)
    """

    # Material properties for Toray 3960 material system (with T1100G UD fibers)
    # https://www.toraycma.com/files/library/0b1f847f1abb9142.pdf
    # all MPa
    Xt = 3896 # tensile strength, fiber direction
    Yt = 79.1 # tensile strength, transverse dir
    Xc = 1868 # compressive strength, fiber direction
    Yc = 265 # compressive strength, transverse dir
    ShearYield = 124 # shear strength in plane, ultimate
    E11 = 173000 # tensile modulus, fiber dir
    E22 = 9500 # tensile modulus, transverse dir
    G12 = 6200 # shear modulus in plane
    v12 = 0.3 # poisson's ratio

    # define constants from Tsai-Wu failure criterion
    F11 = 1/Xt/Xc
    F22 = 1/Yt/Yc
    F12 = -(1/2)*(Xt*Xc*Yt*Yc)**(-1/2)
    F66 = ShearYield**-2
    F1 = 1/Xt - 1/Xc
    F2 = 1/Yt - 1/Yc
    F6 = 0.0

    rho = 1573 # material density, kg/m3

    P = pressure_design * 2.35 # tank pressure, MPa, safety factor 2.35
    # define stress state in the cylindrical portion of the tank
    load_hoop = P*radius # MN / m 
    load_ax = P*radius/2 # MN / m

    # define a symmetric and balanced laminate with one design angle
    thetas = [theta_design, -theta_design, -theta_design, theta_design]

    # Composite analysis notes from:
    # http://web.mit.edu/course/3/3.11/www/modules/laminates.pdf

    # theta is defined from the hoop direction (90 degrees from the axial) as illustrated below
    #
    #[cylinder wall]
    #======================= ) 
    #  |theta/                  )
    #  |    /                    )
    #  |   /                      )  [tank end cap]
    #  |  /                      )
    #  | /                      )  
    #======================= )

    t_one_layer = t_overall / 4

    # define ply interface height (0.0 is midline)                
    z = np.array([-t_overall/2, -t_one_layer, 0.0,
                   t_one_layer,  t_overall / 2])  # meters
                    
    # material property for carbon fiber system in the fiber axis
    Q11 = E11 ** 2 / (E11 - v12**2 * E22)
    Q12 = v12 * E11 * E22 / (E11 - v12**2 * E22)
    Q22 = E11 * E22 / (E11 - v12**2 * E22)
    Q66 = G12

    # ply stiffness matrix in fiber axis
    D = np.array([[Q11, Q12, 0.],[Q12, Q22, 0.],[0., 0., Q66]])
    # ply compliance matrix in fiber axis
    S = np.array([[1/E11, -v12/E11, 0.],[-v12/E11, 1/E22, 0.],[0.,0.,1/G12]])
    # these are inverses of each other

    # assemble the ABD matrices of the laminate stack
    Amat = np.zeros((3,3))
    Bmat = np.zeros((3,3))
    Dmat = np.zeros((3,3))
    for i in range(len(thetas)):
        theta_this_layer = thetas[i]
        A, Ainv, R, Rinv = rotations(theta_this_layer)

        # strain in xy axis = Sbar * stress in xy
        # Sbar = R @ Ainv @ Rinv @ S @ A

        # stress in xy axis = Dbar * strain in xy 
        Dbar = Ainv @ D @ R @ A @ Rinv
        Amat += Dbar * (z[i+1]-z[i])
        # these are zero in this balanced symm layups
        Bmat += 0.5 * Dbar * (z[i+1]**2-z[i]**2)
        # don't care about the moments for this problem
        Dmat += (1/3) * Dbar * (z[i+1]**3-z[i]**3)

    if np.sum(np.abs(Bmat)) > 1e-12:
        raise ValueError('B matrix is nonzero - something is wrong in laminate calculations')

    # N = Nx, Ny, Nxy
    # define the load state on the laminate in terms of hoop load (force / meter) and axial load
    N = np.array([load_hoop, load_ax, 0.0])

    # get the whole-laminate strain state
    strain_xy = np.linalg.inv(Amat) @ N

    # get stress and compute Tsai-Wu failure criterion in one layer (they'll all be identical)
    TW_failure = np.zeros(4)
    for i in range(4):
        theta_this_layer = thetas[i]
        A, Ainv, R, Rinv = rotations(theta_this_layer)
        stress_mat_axes = D @ R @ A @ Rinv @ strain_xy
        s_1 = stress_mat_axes[0]
        s_2 = stress_mat_axes[1]
        t_12 = stress_mat_axes[2]
        TW = (F1 * s_1 + 
            F2 * s_2 + 
            F11 * s_1 * s_1 + 
            F22 * s_2 * s_2 +
            2 * F12 * s_1 * s_2 + 
            F66 * t_12**2)
        TW_failure[i] = TW
    # these should all be uniform unless something is wrong
    if np.sum(np.abs(TW_failure-np.mean(TW_failure))) > 1e-10:
        raise ValueError('All the failure criteria in each ply should be identical for this layup')
    # compute tank weight
    surface_area = np.pi * (length * 2 * radius + 4 * radius ** 2)
    weight = surface_area * t_overall * rho
    return TW_failure, weight

class TankComp(om.ExplicitComponent):
    """
    A thin wrapper of the analyze_tank_structure method
    Finite difference derivatives because the model is cheap
    """

    def setup(self):
        self.add_input('theta_design', val=20.0, units='deg')
        self.add_input('thickness', val=0.04, units='m')
        self.add_input('length', val=1.0, units='m')
        self.add_input('radius', val=0.25, units='m')
        self.add_input('design_pressure', val=70.0, units='MPa')
        self.add_input('hydrogen_density', val=42.0, units='kg/m**3')

        self.add_output('tsai_wu_failure', val=np.ones((4,)))
        self.add_output('tank_weight', val=10., units='kg')
        self.add_output('vol_wt_ratio', val=1.0, units='m**3/kg')
        self.add_output('volume', val=1.0, units='m**3')
        self.add_output('grav_pct', val=0.05)
        self.add_output('hydrogen_weight', val=1., units='kg')

        self.declare_partials(['*'],['*'], method='fd',step=1e-6)

    def compute(self, inputs, outputs):
        # do the structural analysis
        tws, tank_weight = analyze_tank_structure(inputs['theta_design'][0]*np.pi/180, 
                                 inputs['thickness'][0], 
                                 inputs['length'][0], 
                                 inputs['radius'][0],
                                 inputs['design_pressure'][0]
                                 )
        outputs['tsai_wu_failure'] = tws
        outputs['tank_weight'] = tank_weight

        # compute some auxiliary numbers
        # tank volume
        volume = np.pi * (inputs['radius']**2 * inputs['length'] + 4/3*inputs['radius']**3)
        outputs['volume'] = volume

        # volume to weight ratio
        outputs['vol_wt_ratio'] = volume / tank_weight
        
        # fuel weight in the tank at design pressure
        H2_wt = volume * inputs['hydrogen_density'][0]
        outputs['hydrogen_weight'] = H2_wt

        # gravimetric percentage (fuel / tank weight)
        outputs['grav_pct'] = H2_wt / tank_weight
        
if __name__ == "__main__":
    
    # The purpose of this model is to find the optimal 
    # design of a composite wrapped Hydrogen tank
    # We are only considering the tank structure for simplicity
    # (not liner or fittings, which are significant)

    prob = om.Problem()

    # the design pressure is an important parameter
    ivc = om.IndepVarComp()
    ivc.add_output('design_pressure', val=70., units='MPa')
    prob.model.add_subsystem('iv', ivc)

    # compute hydrogen density at 298K using a curve fit
    # data from https://h2tools.org/hyarc/hydrogen-data/hydrogen-density-different-temperatures-and-pressures
    h2density = om.MetaModelStructuredComp()
    h2density.add_input('pressure', 
                        training_data=np.array([0.1, 1., 5., 10., 30., 50., 70., 100.]), 
                        units='MPa')
    h2density.add_output('density', 
                         training_data=np.array([0.0813, 0.8085, 3.9490, 7.6711, 20.537, 30.811, 42.0, 49.424]), 
                         units='kg/m**3')
    prob.model.add_subsystem('h2density', h2density)
    
    # do the structural and weight analysis
    prob.model.add_subsystem('tank', TankComp())

    prob.model.connect('iv.design_pressure', ['tank.design_pressure','h2density.pressure'])
    prob.model.connect('h2density.density', 'tank.hydrogen_density')

    # setup the optimization
    prob.driver = om.ScipyOptimizeDriver()
    prob.driver.options['optimizer'] = 'SLSQP'

    # the tank fiber orientation angle in degrees
    prob.model.add_design_var('tank.theta_design', lower=0.0, upper=75., ref=40.)
    # tank wall thickness in m
    # this must be minimum-gauged for crash and projectile resistance hence 1cm lower bound
    prob.model.add_design_var('tank.thickness', lower=0.01, upper=0.1, ref=0.02)
    # tank length (wants to be upper-bounded, so comment out for now)
    # prob.model.add_design_var('tank.length', lower=0.2, upper=2.0, ref=1.0)
    
    # radius of the cylinder in m
    prob.model.add_design_var('tank.radius', lower=0.02, upper=20.0, ref=0.2)
    # tank design pressure in MPa (multiply by 10 to get bar)
    prob.model.add_design_var('iv.design_pressure', lower=2., upper=70., ref=10.)

    # minimize the weight of the tank subject to structural sizing and fuel wt requirement
    prob.model.add_objective('tank.tank_weight', scaler=0.05)
    prob.model.add_constraint('tank.tsai_wu_failure', upper=np.ones((4,)))
    prob.model.add_constraint('tank.hydrogen_weight', lower=100.0)

    prob.setup()

    # save some parameters
    thicknesses = []
    weights = []
    radii = []
    pressures = []
    fiber_angles = []
    h2_weights = []

    # run a series of optimizations at different tank lengths and see what happens 
    # to optimal pressure and tank radius
    lengths = np.linspace(0.5,2.0,10)
    for length in lengths:
        prob['tank.length'] = length
        prob.run_driver()
        thicknesses.append(prob['tank.thickness'][0].copy())
        weights.append(prob['tank.tank_weight'][0].copy())
        radii.append(prob['tank.radius'][0].copy())
        pressures.append(prob['iv.design_pressure'][0].copy())
        fiber_angles.append(prob['tank.theta_design'][0].copy())
        h2_weights.append(prob['tank.hydrogen_weight'][0].copy())

    # plot a bunch of design parameters and outputs
    fig, axs = plt.subplots(3, 2, sharex='col')
    axs[0, 0].plot(100*lengths, 1000*np.array(thicknesses))
    axs[0, 0].set(ylabel='Wall thickness (mm)')
    axs[0, 1].plot(100*lengths, np.array(weights))
    axs[0, 1].set(ylabel='Cylinder weight (kg)')
    axs[1, 0].plot(100*lengths, 100*np.array(radii))
    axs[1, 0].set(ylabel='Cylinder radius (cm)')
    axs[1, 1].plot(100*lengths, 10*np.array(pressures))
    axs[1, 1].set(ylabel='Tank pressure (bar)')
    axs[2, 0].plot(100*lengths, np.array(fiber_angles))
    axs[2, 0].set(ylabel='Fiber angle (deg)')
    axs[2, 0].ticklabel_format(useOffset=False, style='plain')
    axs[2, 0].set_ylim(35, 36)
    axs[2, 1].plot(100*lengths, np.array(h2_weights))
    axs[2, 1].set(ylabel='Hydrogen weight (kg)')
    axs[2, 1].ticklabel_format(useOffset=False, style='plain')
    axs[2, 1].set_ylim(99, 101)
    axs[2, 0].set(xlabel='Tank cylinder length (cm)')
    axs[2, 1].set(xlabel='Tank cylinder length (cm)')    
    plt.tight_layout()
    plt.show()

    print('Pressure ', prob['iv.design_pressure'])
    print('Theta ', prob['tank.theta_design'])
    print('Thickness ', 1000*prob['tank.thickness'], ' mm')
    print('Length ', prob['tank.length'])
    print('Radius ', prob['tank.radius'])
    print('Tank Weight ', prob['tank.tank_weight'])
    print('Grav pct ', prob['tank.grav_pct'])
    print('Failure crit ', prob['tank.tsai_wu_failure'])
    print('Tank Volume', prob['tank.volume'])
    prob.model.list_outputs(units=True)

Great work