{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Aluminum and Composite Fuel Tank Analysis\n",
"\n",
"Ondrej Fercak, ofercak@pdx.edu\n",
"\n",
"Erin Schmidt, esch2@pdx.edu\n",
"\n",
"Ian Zabel, ibz@pdx.edu\n",
"\n",
"@psas\n",
"\n",
"### Introduction\n",
"The purpose of this analysis is to determine the failure criterea of the aluminum or carbon fiber composite fuel tank walls. An Abaqus 6.14 model is developed assuming that the fuel tank is axisymmetric with some simplifications regarding the internals of the fuel tank. An initial model is developed and run, assuming there is no baffling or plumbing within the fuel tank, to ensure the composite definitions are sensible. A more advanced model is then developed using Solidworks as well as Abaqus that includes the internals expected to be present.\n",
"\n",
"**Internals:**\n",
"- Fuel line\n",
"- Baffling fins\n",
"- Baffling net\n",
"- Transition region ribbing\n",
"\n",
"**Propellant:**\n",
"- **Fuel:** Isopropyl, 70%\n",
"- **Oxidizer:** Liquid Oxygen\n",
"\n",
"### Aluminum Theory\n",
"For the propellant, the oxidizer to fuel ratio $OF$ is set as $1.3$, leading to the following expressions:\n",
"\n",
"$$\\frac{m_{ox}}{m_{f}}=OF=1.3$$\n",
"\n",
"With the mass of the propellant defined by the [LV4 Design Optimization](https://github.com/psas/liquid-engine-analysis/blob/master/optimization/LV4_Optimization.ipynb), the required volume for each of the sections is defined by the following set of expressions, assuming both the fuel and oxidizer tanks have the same hemishperical endcaps:\n",
"\n",
"$$m_{p}=m_{ox}+m_{f}$$\n",
"$$V_{ox} = \\rho_{ox} m_{ox} = 2 \\left(\\frac{4}{3} \\pi\\right)\\left(R_{o}-t\\right)^3 + \\pi L_{ox} \\left(R_{o}-t\\right)^2$$\n",
"$$V_{f} = \\rho_{f} m_{f} = 2 \\left(\\frac{4}{3} \\pi\\right)\\left(R_{o}-t\\right)^3 + \\pi L_{f} \\left(R_{o}-t\\right)^2$$\n",
"\n",
"The goal with using these definitions is to drive the cylindrical tank portion lengths $L_{ox}$ and $L_{f}$ using the ratio $OF=1.3$ and known propellant densities $\\rho_{ox}$ and $\\rho_{f}$. Rearranging to solve for the tank lengths leads to the following expressions:\n",
"\n",
"$$L_{ox}(t)=\\frac{\\rho_{ox} \\left(m_{p}-\\frac{m_{p}}{1+OF}\\right)}{\\pi (R_{o}-t)^2}-2 \\frac{4}{3}(R_{o}-t)$$\n",
"$$L_{f}(t)=\\frac{\\rho_{ox} \\left(\\frac{m_{p}}{1+OF}\\right)}{\\pi (R_{o}-t)^2}-2 \\frac{4}{3}(R_{o}-t)$$\n",
"\n",
"With the lengths of the cylindrical portions defined as functions of the wall thickness, the next step is to determine the wall minimum thickness with a factor of safety. To do this, the properties of the material are considered.\n",
"\n",
"| Material | Failure Criteria |\n",
"| -------- | ---------------- |\n",
"| Ductile | Von Mises |\n",
"| Brittle | Mohr's Theory |\n",
"\n",
"$$\\sqrt{\\frac{(\\sigma_{1}-\\sigma_{2})^2 + (\\sigma_{2}-\\sigma_{3})^2 + (\\sigma_{3}-\\sigma_{1})^2}{2}} \\le \\sigma_{y}$$\n",
"\n",
"### Aluminum Assumptions\n",
"To perform this analysis using composite or aluminum, the pressure vessel is assumed to be thin-walled. To apply this assumption the wall thickness $t$ must be much smaller than a characteristic dimension $R$, in this case the tank radius $r$.\n",
"\n",
"$$\\frac{R}{t}>10$$\n",
"\n",
"The simplifed pressure vessel is assumed to be axisymmetric, and the advanced model vessel is assumed to by symmetric about both the $x,z$ and $y,z$ planes. To simplify manufacturing and design, both the oxidizer and fuel tanks are the same material, wall thickness, radius, and have the same hemispherical endcaps. One of the tanks will have a fuel line running through the center, while the other tank will only have an opening for that specific fuel line at the bottom.\n",
"\n",
"### Composite Overwrap Pressure Vessel\n",
"\n",
"**General Structure:**\n",
"- Liner\n",
"- Hoop wrap\n",
"- Helical wrap\n",
"\n",
"### Composite Theory\n",
"\n",
"([K]+lambdai[S]){psi}i={0}\n",
"K=stiffness\n",
"S=stress stiffness\n",
"lambdai=ith eigenvalue\n",
"psii=ith eigenvector of displacements\n",
"'Block Lanczos'\n",
"\n",
"Regarding the stress and strain of the composite, the specific failure criterea of the fibers and the matrix are considered. The elastic portion of the composite stress-strain curve is simply a combination of the elastic modulus for the fibers and matrix. The plastic portion will be linear, failing when the fiber fails.\n",
"\n",
"$$\\sigma_{composite}^\\star=\\sigma_{fiber}^\\star$$\n",
"\n",
"$$E_{comp}=E_{m} V_{m} + K E_{f} V_{f}, K=1$$\n",
"\n",
"$$\\frac{1}{E_{comp,trans}}=\\frac{V_{m}}{V_{comp}} \\frac{1}{E_{m}} + \\frac{V_{f}}{V_{comp}} \\frac{1}{E_{f}}$$\n",
"\n",
"The failure stresses for the composite in longitudinal and transverse directions are represented by the following equations:\n",
"\n",
"$$\\sigma_{comp,long}^\\star=\\sigma_{m}^p \\frac{1-V_{f}}{V_{c}} + \\sigma_{f}^p \\frac{V_{f}}{V_{c}}$$\n",
"\n",
"$$\\sigma_{comp,trans}^\\star=E_{comp} \\frac{V_{m}}{V_{comp}} \\epsilon_{m} + E_{comp} \\frac{V_{f}}{V_{comp}} \\epsilon_{f}$$\n",
"\n",
"### FEA Procedure\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\n",
"FUEL TANK PARAMETERS\n",
"--------------------\n",
"Tank Radius = 0.13 m\n",
"Tank Wall Thickness = 2.00 mm\n",
"Volume, lox = 0.035 m^3\n",
"Volume, fuel = 0.545 m^3\n",
"\n",
"FUEL TANK MASS\n",
"--------------\n",
"Mass, Lox = 39.57 kg\n",
"Mass, Eth. = 30.43 kg\n",
"Mass, Prop. = 70.00 kg\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import math as m\n",
"import numpy as np\n",
"import matplotlib\n",
"import matplotlib.pyplot as plot\n",
"import pylab\n",
"from matplotlib import rc\n",
"%matplotlib inline\n",
"%config InlineBackend.figure_formats=['svg']\n",
"\n",
"# Ondrej Fercak\n",
"# Ian Zabel\n",
"# Note: - Dimensions from LV4 MDO\n",
"# - Atmo pressure from LV4 MDO\n",
"# - All calculations assume hemisphere endcaps\n",
"\n",
"# Dimensions and Material Properties\n",
"\n",
"Do = 0.254 # Outer tank diameter (m)\n",
"Ro = Do/2 # Outer tank radius (m)\n",
"L = 1.2 # Overall tank length (m)\n",
"TW = 15 # Thin-walled assumption limit\n",
"t = 0.002 # Wall thickness (m)\n",
"FS = 1.1 # Factor of safety\n",
"Sy = 240e6 # Yield strength, Al (Pa)\n",
"Su = 290e6 # Ultimate tensile strength, Al (Pa)\n",
"TL = Ro*0.3 # Transition region length (m)\n",
"\n",
"# Propellant Parameters\n",
"# Note: Volume values TEMP\n",
"\n",
"OF = 1.3 # Oxidizer to fuel ratio\n",
"rhof = 852.3 # Density, ethanol fuel (kg/m^3)\n",
"rhoo = 1141 # Density, lox (kg/m^3)\n",
"\n",
"# Pressure\n",
"# Note: Use either np.all() or np.any() for ineq. def\n",
"\n",
"def atmo(x):\n",
" if np.all(x < 11000):\n",
" T = 15.04-0.0065*x\n",
" P = 101.3*((T+273.1)/288.1)**5.26\n",
" elif np.all(11000 <= x) & np.all(x < 25000):\n",
" T = -56.46\n",
" P = 22.65*exp(1.73-0.00016*x) \n",
" else:\n",
" T = -131.2+0.003*x\n",
" P = 2.488*((T+273.1)/216.6)**(-11.4)\n",
"\n",
" rho = P/(0.29*(T+273.1)) # Density, ambient air (kg/m^3)\n",
" Pa = P*1000 # Pressure, ambient air (Pa)\n",
" Ta = T+273.1 # Temperature, ambient air (K)\n",
" return Pa, rho, Ta\n",
"\n",
"x = np.arange(0, 1e5, 1) # Altitude (m)\n",
"Pa = atmo(x)[0] # Atmospheric pressure (Pa)\n",
"Po = 50-Pa # Pressure, lox, abs. (Pa)\n",
"Pf = 50-Pa # Pressure, fuel, abs. (Pa)\n",
"\n",
"# Required Mass\n",
"\n",
"mp = 70 # Mass, propellant (kg)\n",
"mf = mp/(1+OF) # Mass, fuel (kg)\n",
"mo = mp-mp/(1+OF) # Mass, lox (kg)\n",
"\n",
"# Dome Volume\n",
"\n",
"def dome(R):\n",
" V = (4/3)*m.pi*R**3\n",
" return V\n",
"\n",
"# Cylinder volume\n",
"\n",
"def cyl(R, L):\n",
" V = m.pi*L*R**2\n",
" return V\n",
"\n",
"Endo = 2*dome(Ro) # End cap volume, lox (m^3)\n",
"Endf = 2*dome(Ro) # End cap volume, fuel (m^3)\n",
"\n",
"CLo = rhoo*mo/cyl(Ro-t, L)-(8/3)*(Ro-t) # Cyl. length, lox (m)\n",
"CLf = rhof*mf/cyl(Ro-t, L)-(8/3)*(Ro-t) # Cyl. length, fuel (m)\n",
"Vo = mo/rhoo # Volume, lox (m^3)\n",
"Vf = mf/rhof # Volume, fuel (m^3)\n",
"\n",
"mo = rhoo*Vo # Mass, lox (kg)\n",
"mf = rhof*Vf # mass, fuel (kg)\n",
"mp = mo+mf # Mass, propellant (kg) \n",
"\n",
"# Stress, cylindrical portion\n",
"\n",
"Sho = Po*Ro/t # Hoop stress, cyl, lox (Pa)\n",
"Shf = Pf*Ro/t # Hoop stress, cyl, fuel (Pa)\n",
"Sao = 0.5*Sho # Axial stress, cyl, lox (Pa)\n",
"Saf = 0.5*Shf # Axial stress, cyl, fuel (Pa)\n",
"\n",
"# Stress, endcaps\n",
"# Note: Static equilibrium, stress net resultant\n",
"\n",
"Sno = Sao # Stress, endcaps, lox (Pa)\n",
"Snf = Saf # Stress, endcaps, Fuel (Pa)\n",
"\n",
"# Composite overwrap\n",
"\n",
"# Parameters\n",
"# Note: typical El = 130e9 Pa\n",
"# Et = 11e9 Pa\n",
"# v = 0.27\n",
"\n",
"# Expected values, carbon composite (TEMP): \n",
"# Vf = 0.55\n",
"# E22 = 130 GPa\n",
"# E11 = 11 GPa\n",
"# G12 = 5.5 GPa\n",
"# v12 = 0.27\n",
"\n",
"rhofib = 2500 # Density, fiber (kg/m^3) (TEMP)\n",
"rhom = 1611 # Density, matrix (kg/m^3) (TEMP)\n",
"rhopla = 1000 # Density, internal plastic (kg/m^3) (TEMP)\n",
"wf = 0.65 # Weight, fiber (kg) (TEMP)\n",
"wm = 0.35 # Weight, matrix (kg) (TEMP)\n",
"Ef = 106.3e9 # Modulus, fiber (Pa)\n",
"Em = 3.4e9 # Modulus, matrix (Pa)\n",
"\n",
"Vf = (wf/rhofib)/(wf/rhofib + wm/rhom) # Volume fraction, fiber\n",
"Vm = 1-Vf # Volume fraction, matrix\n",
"\n",
"rhocom = rhofib*Vf + rhom*Vm # Density, composite (kg/m^3)\n",
"\n",
"vf = 0.26 # Poisson's ratio, fiber (TEMP)\n",
"vm = 0.33 # Poisson's ratio, matrix (TEMP)\n",
"\n",
"Gf = Ef/(2*(1+vf)) # Shear modulus, fiber (Pa)\n",
"Gm = Em/(2*(1+vm)) # Shear modulus, matrix (Pa)\n",
"\n",
"# Weave Angle: Orientation of the specific referenced ply\n",
"# Wraps: - Cylindrical, one at 0 deg.\n",
"# - Helical, one at 45 deg., one at -45 deg.\n",
"\n",
"Th = 45 # Weave angle (deg.)\n",
"\n",
"# Major Poisson's: Ratio of trans. contraction strain to long. extension strain in the direction of the stretch force\n",
"# Tensile --> positive\n",
"\n",
"E22 = Ef*Em/(Em*Vf+Ef*Vm) # Modulus, longitudinal (Pa)\n",
"E11 = Ef*Vf+Em*Vm # Modulus, transverse (Pa)\n",
"v12 = vf*Vf+vm*Vm # Poisson's ratio, major\n",
"v21 = E22*v12/E11 # Poisson's ratio, minor\n",
"G12 = Gf*Gm/(Gf*Vm+Gm*Vf) # Shear modulus\n",
"\n",
"# Windenburg and Trilling Equation\n",
"# Carbon Composite Epoxy Helix at 45 deg.\n",
"\n",
"Ex = E22/(m.cos(Th)**4 + (E22/E11)*m.sin(Th)**4 + 0.25*((E22/E11)-2*v12)*m.sin(2*Th)**2) # Modulus, x (Pa)\n",
"Ey = E22/(m.sin(Th)**4 + (E22/E11)*m.cos(Th)**4 + 0.25*((E22/E11)-2*v12)*m.sin(2*Th)**2) # Modulus, y (Pa) \n",
"Pcr = (2.42*Ex*(t/Do)**2.5)/(((1-v12**2)**0.75)*(L/Do - 0.45*m.sqrt(t/Do))) # Critical pressure (Pa)\n",
"\n",
"# Results\n",
"\n",
"print('\\n')\n",
"print('FUEL TANK PARAMETERS')\n",
"print('--------------------')\n",
"print('Tank Radius = {0:.2f} m'.format(Ro))\n",
"print('Tank Wall Thickness = {0:.2f} mm'.format(t*1000))\n",
"print('Volume, lox = {0:.3f} m^3'.format(Vo))\n",
"print('Volume, fuel = {0:.3f} m^3'.format(Vf))\n",
"print('\\nFUEL TANK MASS')\n",
"print('--------------')\n",
"print('Mass, Lox = {0:.2f} kg'.format(mo))\n",
"print('Mass, Eth. = {0:.2f} kg'.format(mf))\n",
"print('Mass, Prop. = {0:.2f} kg'.format(mp))\n",
"\n",
"# Plotting\n",
"rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})\n",
"rc('text', usetex=True)\n",
"pylab.rcParams['figure.figsize'] = (10.0, 10.0)\n",
"f, (ax1, ax2) = plot.subplots(2, sharex=True)\n",
"ax1.plot(x, Po)\n",
"ax1.yaxis.major.locator.set_params(nbins=6)\n",
"ax1.set_title('LV4, Lox Pressure (kPa) vs. Altitude (km)')\n",
"ax2.plot(x, Pf)\n",
"ax2.yaxis.major.locator.set_params(nbins=6)\n",
"ax2.set_title('LV4, Fuel Pressure (kPa) vs. Altitude (km)')\n",
"plot.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 1
}