c car2.f (version 1.0) May 1998; Christian Fellner
c Genral program for calculating hybrid vehicle mileage.
c You may make a copy of this program which you may personally
c and freely use in its unaltered form. You may distribute this
c program subject to the conditions that it be made freely
c available and that any duplication of this program must be
c essentially unaltered and must include this notice.
c We make no warranties, express or implied, that this program
c is free of errors or that it will meet the requirements of your
c application. The author and publisher disclaim all liability for
c direct or consequential damages resulting from use of this program.
c This program calculates the average mileage of a hybrid vehicle.
PROGRAM CAR2
c Simplified from original car program.
c One overpotential fit: ai^2 + bi + c
c Electrode open-circuit potential fit: ax^2 + bx + c
c State of charge overpotential abjustment: ax^2 + bx + c
c No change in open circuit potential correction for overpotential fit.
c - Better for general battery and easier to measure.
c Type 'cat spec# | car2' to run program and read data file.
c Use specg for general battery.
c Type 'cat spec# | car2 >file' to run program, read data file,
c and write output to specified file.
c Combines a simple vehicle model with a simplified battery model.
c Works only for design driving cycle. (6 1 minute cycles)
c Define variable types for battery & vehicle.
REAL Athick, Cthick, MwA, MwC, Aact, Cact, Adensity, Cdensity
REAL ULimit, LLimit, OcA1, OcA2, OcA3, OcC1, OcC2, OcC3
REAL C1, C2, Area, Pbrake, SepStep, SocStep
REAL Apower1a, Apower1b, Apower2a, Apower2b
REAL Apower3a, Apower3b, Apower4a, Apower4b
REAL Vmass, Meng, Emass, Pmass, Bmass
REAL Peff, Eeff, Geff, Beff, Veng
c Define variable types for battery linear approximations.
c State of charge overpotential adjustments.
REAL ScA1, ScA2, ScA3
REAL ScC1, ScC2, ScC3
c Overpotential fit.
REAL Op1, Op2, Op3
c Define variables used internal to the program.
REAL Cpower1, Cpower2, Cpower3, Cpower4
REAL Apower1, Apower2, Apower3, Apower4
REAL Bpower1, Bpower2, Bpower3, Bpower4
REAL Bap1, Bap2, Bap3, Bap4, Bcp1, Bcp2, Bcp3, Bcp4
REAL Bbp1, Bbp2, Bbp3, Bbp4, Bsp
REAL Mass, SocA, SocC, Engine, MA, MC, BatOut, BatIn, BatEff
REAL Delta, Delta1, Delta2, Delta3, Delta4, Xi, Yi, OCP
REAL Separator, Smin, Smax, Vel1, Vel2, Vel3, Vel4
REAL I1, I2, I3, I4, I5, I6, I7, I8, I9, I10, I11, I12, I13
REAL I14, I15, I16, I17, I18, I19, I20, I21, I22, I23, I24
REAL V1, V2, V3, V4, V5, V6, V7, V8, V9, V10, V11, V12, V13
REAL V14, V15, V16, V17, V18, V19, V20, V21, V22, V23, V24
REAL GRAV, AIR, SoccStop, Dump
c Careful with type declarations on different systems (real).
INTEGER FARADY, EGAS, Number, t1, t2, t3, t4, t5, t6, tb, ts
c Enter in constant parameter values.
c FARADY = Faraday's constant; coulombs/equivalent
c GRAV = gravitational constant; m/s2
c EGAS = energy content of gasoline; kJ/L
c AIR = density of air; kg/m3
PARAMETER (FARADY = 96487, GRAV = 9.81)
PARAMETER (EGAS = 31810, AIR = 1.202)
c Use common statement so can get state of charge factors into subroutine
c TUNE, FUEL, OVERPOTENTIAL, OPEN easier, must name since more than 1 common.
COMMON /ad/ ScA1, ScA2, ScA3, ScC1, ScC2, ScC3
COMMON /gas/ C1, C2, Area, Peff, Eeff, Geff, Veng
COMMON /op/ Op1, Op2, Op3
COMMON /oc/ OcA1, OcA2, OcA3, OcC1, OcC2, OcC3
c Read in variables from data file.
READ *, Athick
READ *, Cthick
READ *, Aact
READ *, Cact
READ *, Bmass
READ *, Pbrake
READ *, C1
READ *, C2
READ *, Area
READ *, Adensity
READ *, Cdensity
READ *, Vmass
READ *, Apower1a
READ *, Apower1b
READ *, Apower2a
READ *, Apower2b
READ *, Apower3a
READ *, Apower3b
READ *, Apower4a
READ *, Apower4b
READ *, Meng
READ *, Emass
READ *, Pmass
READ *, Peff
READ *, Eeff
READ *, Geff
READ *, Beff
READ *, MwA
READ *, MwC
READ *, ScA1
READ *, ScA2
READ *, ScA3
READ *, ScC1
READ *, ScC2
READ *, ScC3
READ *, Op1
READ *, Op2
READ *, Op3
READ *, OcA1
READ *, OcA2
READ *, OcA3
READ *, OcC1
READ *, OcC2
READ *, OcC3
c Print out program input parameters.
PRINT *, '*** Vehicle model results. ***'
PRINT *, 'Anode electrode thickness (um)=,', Athick
PRINT *, 'Cathode electrode thickness (um) =,', Cthick
PRINT *, 'Anode % volume active material =,', Aact
PRINT *, 'Cathode % volume active material =,', Cact
PRINT *, 'Battery mass (kg/m2) =,', Bmass
PRINT *, 'Vehicle rolling drag coefficient =,', C1
PRINT *, 'Vehicle air drag coefficient =,', C2
PRINT *, 'Vehicle front surface area (m2) =,', Area
PRINT *, 'Base vehicle mass (kg) =,', Vmass
PRINT *, 'Maximum engine design (kw) =,', Meng
PRINT *, 'Engine mass (kw/kg) =,', Emass
PRINT *, 'Passenger mass (kg) =,', Pmass
PRINT *, 'Powertrain efficiency =,', Peff
PRINT *, 'Engine thermal efficiency =,', Eeff
PRINT *, 'Generator efficiency =,', Geff
PRINT *, 'Engine into battery efficiency =,', Beff
PRINT *, 'Regenerative braking % =,', Pbrake
c ********** Begin main body of the program. **********
c Initialize variables used in program, set by user for ranges.
c Smallest & largest separator area and step size used in program.
Smin = 85
Smax = 95
SepStep = 5
c Set maximum/minimum state of charge for anode and step size.
SocA = 0.45
SoccStop = 0.44
SocStep = 0.05
c Enter % of maximum regenerative braking during hard stop.
Dump = 1.0
PRINT *, '% braking during hard stop =,', Pbrake*Dump
PRINT *, ' '
c Initial engine guess for urban cycle & set variation.
c Engine maximum during acceleration, at level during cruising,
c and minimum during braking and rest.
c Values for Veng range from 0 to 1.
c 0 constant urban cruising engine.
c 1 no engine during braking & stop & double during acceleration.
Engine = 4500
Veng = 0.0
c Set voltage limits on battery.
LLimit = 1.6
ULimit = 2.3
c Cruising velocities in m/s.
Vel1 = 11.111
Vel2 = 15.278
Vel3 = 19.444
Vel4 = 25.000
c Driving cycle time segment lengths in seconds.
t1 = 10
t2 = 30
t3 = 12
t4 = 28
t5 = 15
t6 = 25
tb = 5
ts = 15
c Calculate constants used in program. (mol active in anode & cathode)
c Units are mol/m2 of separator area.
MA = Athick*0.000001*Aact*Adensity*1000/MwA
MC = Cthick*0.000001*Cact*Cdensity*1000/MwC
c First big loop, state of charge.
DO WHILE (SocA .GE. SoccStop)
c Fit cathode to same capacity as anode.
SocC = 1 - SocA*MA/MC
c Calculate open circuit potential
OCP = OPEN (SocA, SocC)
c Print state of charge of electrodes.
PRINT *, '****** Begin new state of charge loop. ******'
PRINT *, ' '
PRINT *, 'Initial x in anode =,', SocA
PRINT *, 'Initial y in cathode =,', SocC
PRINT *, 'Open-circuit cell potential =,', OCP
PRINT *, ' '
c Use Xi & Yi to hold initial values for state of charge
Xi = SocA
Yi = SocC
c Need to re-initialize separator area for next loop.
Separator = Smax
c Second loop, separator area.
DO WHILE (Separator .GE. Smin)
c Add in 1.5 factor for battery and engine support.
Mass = (Meng/Emass)*1.5 + Pmass + Vmass + (Bmass*Separator)*1.5
PRINT *, 'Separator area (m2) =,', Separator
PRINT *, 'Battery mass (no support; kg) =,', Bmass*Separator
PRINT *, 'Total vehicle mass (kg) =,', Mass
CALL POWER (SocA, SocC, Overpot, Separator, LLimit)
c Calculate vehicle wheel power requirements for each driving segment.
Apower1 = Apower1a*Mass + Apower1b
Apower2 = Apower2a*Mass + Apower2b
Apower3 = Apower3a*Mass + Apower3b
Apower4 = Apower4a*Mass + Apower4b
Cpower1 = C1*GRAV*Vel1*Mass + 0.5*AIR*C2*Area*Vel1**3
Cpower2 = C1*GRAV*Vel2*Mass + 0.5*AIR*C2*Area*Vel2**3
Cpower3 = C1*GRAV*Vel3*Mass + 0.5*AIR*C2*Area*Vel3**3
Cpower4 = C1*GRAV*Vel4*Mass + 0.5*AIR*C2*Area*Vel4**3
c Need to divide 1/2 by 5 for time. (5 second stop)
Bpower1 = -0.1*Mass*Vel1**2
Bpower2 = -0.1*Mass*Vel2**2
Bpower3 = -0.1*Mass*Vel3**2
Bpower4 = -0.1*Mass*Vel4**2
c Third loop, solve system to find required engine.
c Need to give Delta an initial non-zero value to enter loop.
Delta = 1
c Initialize engine counting loop.
Number = 1
DO WHILE (ABS(Delta) .GE. 0.01)
c Counting cycle for engine loop.
IF (Number .LE. 100) THEN
Number = Number + 1
ELSE
PRINT *, 'I have been in the engine loop too long!'
STOP
ENDIF
c Calculate battery power requirements at leads for each driving segment.
c For acceleration battery power always assume draining.
c Battery power positive to supply power.
c Battery power negative to accept power.
Bap1 = (Apower1 - Engine*(1+Veng)*Peff)/Geff
Bap2 = (Apower2 - Engine*(1+Veng)*Peff)/Geff
Bap3 = (Apower3 - Engine*(1+Veng)*Peff)/Geff
Bap4 = (Apower4 - Engine*(1+Veng)*Peff)/Geff
c For cruising power must check if draining/charging.
IF (Cpower1 .GT. Engine*Peff) THEN
Bcp1 = (Cpower1 - Engine*Peff)/Geff
ELSE
Bcp1 = (Cpower1/Peff - Engine)*Beff
ENDIF
IF (Cpower2 .GT. Engine*Peff) THEN
Bcp2 = (Cpower2 - Engine*Peff)/Geff
ELSE
Bcp2 = (Cpower2/Peff - Engine)*Beff
ENDIF
IF (Cpower3 .GT. Engine*Peff) THEN
Bcp3 = (Cpower3 - Engine*Peff)/Geff
ELSE
Bcp3 = (Cpower3/Peff - Engine)*Beff
ENDIF
IF (Cpower4 .GT. Engine*Peff) THEN
Bcp4 = (Cpower4 - Engine*Peff)/Geff
ELSE
Bcp4 = (Cpower4/Peff - Engine)*Beff
ENDIF
c Regenerative braking based on energy driving at battery leads.
Bbp1 = -Engine*(1-Veng)*Beff + Bpower1*Pbrake
Bbp2 = -Engine*(1-Veng)*Beff + Bpower2*Pbrake
Bbp3 = -Engine*(1-Veng)*Beff + Bpower3*Pbrake
Bbp4 = -Engine*(1-Veng)*Beff + Bpower4*Pbrake*Dump
Bsp = -Engine*(1-Veng)*Beff
c Divide all battery powers by separator area to get per unit area.
Bap1 = Bap1/Separator
Bap2 = Bap2/Separator
Bap3 = Bap3/Separator
Bap4 = Bap4/Separator
Bcp1 = Bcp1/Separator
Bcp2 = Bcp2/Separator
Bcp3 = Bcp3/Separator
Bcp4 = Bcp4/Separator
Bbp1 = Bbp1/Separator
Bbp2 = Bbp2/Separator
Bbp3 = Bbp3/Separator
Bbp4 = Bbp4/Separator
Bsp = Bsp/Separator
c Determine current requirements during each driving segment.
CALL CURRENT (Bap1,I1,V1,SocA,SocC,Overpot,t1,MA,MC)
CALL CURRENT (Bcp1,I2,V2,SocA,SocC,Overpot,t2,MA,MC)
CALL CURRENT (Bbp1,I3,V3,SocA,SocC,Overpot,tb,MA,MC)
CALL CURRENT (Bsp,I4,V4,SocA,SocC,Overpot,ts,MA,MC)
CALL CURRENT (Bap3,I5,V5,SocA,SocC,Overpot,t3,MA,MC)
CALL CURRENT (Bcp3,I6,V6,SocA,SocC,Overpot,t4,MA,MC)
CALL CURRENT (Bbp3,I7,V7,SocA,SocC,Overpot,tb,MA,MC)
CALL CURRENT (Bsp,I8,V8,SocA,SocC,Overpot,ts,MA,MC)
CALL CURRENT (Bap1,I9,V9,SocA,SocC,Overpot,t1,MA,MC)
CALL CURRENT (Bcp1,I10,V10,SocA,SocC,Overpot,t2,MA,MC)
CALL CURRENT (Bbp1,I11,V11,SocA,SocC,Overpot,tb,MA,MC)
CALL CURRENT (Bsp,I12,V12,SocA,SocC,Overpot,ts,MA,MC)
CALL CURRENT (Bap4,I13,V13,SocA,SocC,Overpot,t5,MA,MC)
CALL CURRENT (Bcp4,I14,V14,SocA,SocC,Overpot,t6,MA,MC)
CALL CURRENT (Bbp4,I15,V15,SocA,SocC,Overpot,tb,MA,MC)
CALL CURRENT (Bsp,I16,V16,SocA,SocC,Overpot,ts,MA,MC)
CALL CURRENT (Bap2,I17,V17,SocA,SocC,Overpot,t3,MA,MC)
CALL CURRENT (Bcp2,I18,V18,SocA,SocC,Overpot,t4,MA,MC)
CALL CURRENT (Bbp2,I19,V19,SocA,SocC,Overpot,tb,MA,MC)
CALL CURRENT (Bsp,I20,V20,SocA,SocC,Overpot,ts,MA,MC)
CALL CURRENT (Bap1,I21,V21,SocA,SocC,Overpot,t1,MA,MC)
CALL CURRENT (Bcp1,I22,V22,SocA,SocC,Overpot,t2,MA,MC)
CALL CURRENT (Bbp1,I23,V23,SocA,SocC,Overpot,tb,MA,MC)
CALL CURRENT (Bsp,I24,V24,SocA,SocC,Overpot,ts,MA,MC)
c Determine if input current = output current.
Delta1 = (I1 + I9 + I21)*t1 + (I2 + I10 + I22)*t2
Delta2 = (I5 + I17)*t3 + (I6 + I18)*t4 +I13*t5 + I14*t6
Delta3 = (I3 + I7 + I11 + I15 + I19 + I23)*tb
Delta4 = (I4 + I8 + I12 + I16 + I20 + I24)*ts
Delta = Delta1 + Delta2 + Delta3 + Delta4
IF (ABS(Delta) .GE. 0.01) THEN
c Recalculate engine and try again. Do better than this.
Engine = Engine + (Delta*OCP*Separator/360)
ENDIF
c Set states of charge back to initial loop values before re-enter.
Socc = Xi
Socmn = Yi
ENDDO
c Calculate battery cycle efficiency, done at battery leads.
c Efficiency = (Total battery energy output) / (Total battery energy input)
c Know braking and stopped charge battery & acceleration drains.
BatIn = ABS((3*Bbp1 + Bbp2 + Bbp3 + Bbp4)*tb + 6*Bsp*ts)
BatOut = 3*Bap1*t1 + (Bap2 + Bap3)*t3 + Bap4*t5
c Need to determine with logical statements if cruising charging/draining.
c Logical statement required for each cruising segment.
IF (Bcp1 .GE. 0) THEN
BatOut = BatOut + 3*Bcp1*t2
ELSE
BatIn = BatIn + ABS(3*Bcp1*t2)
ENDIF
IF (Bcp2 .GE. 0) THEN
BatOut = BatOut + Bcp2*t4
ELSE
BatIn = BatIn + ABS(Bcp2*t4)
ENDIF
IF (Bcp3 .GE. 0) THEN
BatOut = BatOut + Bcp3*t4
ELSE
BatIn = BatIn + ABS(Bcp3*t4)
ENDIF
IF (Bcp4 .GE. 0) THEN
BatOut = BatOut + Bcp4*t6
ELSE
BatIn = BatIn + ABS(Bcp4*t6)
ENDIF
BatEff = BatOut/BatIn
c Print out important results of program, after get proper engine & currents.
PRINT *, 'Battery cycle efficiency =,', BatEff
PRINT *, 'Base Urban engine size (watts) =,', Engine
CALL FUEL (Engine, Mass)
PRINT *, ' '
c Print headings for output from subroutine CURRENT.
PRINT *, 'Average,Average Cell, Battery'
PRINT *, 'Current,Potential, Power'
PRINT *, '(A/m2),(volts), (W/m2)'
PRINT *, I1,',',V1,',',I1*V1
PRINT *, I2,',',V2,',',I2*V2
PRINT *, I3,',',V3,',',I3*V3
PRINT *, I4,',',V4,',',I4*V4
PRINT *, I5,',',V5,',',I5*V5
PRINT *, I6,',',V6,',',I6*V6
PRINT *, I7,',',V7,',',I7*V7
PRINT *, I8,',',V8,',',I8*V8
PRINT *, I9,',',V9,',',I9*V9
PRINT *, I10,',',V10,',',I10*V10
PRINT *, I11,',',V11,',',I11*V11
PRINT *, I12,',',V12,',',I12*V12
PRINT *, I13,',',V13,',',I13*V13
PRINT *, I14,',',V14,',',I14*V14
PRINT *, I15,',',V15,',',I15*V15
PRINT *, I16,',',V16,',',I16*V16
PRINT *, I17,',',V17,',',I17*V17
PRINT *, I18,',',V18,',',I18*V18
PRINT *, I19,',',V19,',',I19*V19
PRINT *, I20,',',V20,',',I20*V20
PRINT *, I21,',',V21,',',I21*V21
PRINT *, I22,',',V22,',',I22*V22
PRINT *, I23,',',V23,',',I23*V23
PRINT *, I24,',',V24,',',I24*V24
c Check to see if outside voltage limits. (Has to be better way to do this)
IF ((V1 .LE. LLimit) .OR. (V1 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V2 .LE. LLimit) .OR. (V2 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V3 .LE. LLimit) .OR. (V3 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V4 .LE. LLimit) .OR. (V4 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V5 .LE. LLimit) .OR. (V5 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V6 .LE. LLimit) .OR. (V6 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V7 .LE. LLimit) .OR. (V7 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V8 .LE. LLimit) .OR. (V8 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V9 .LE. LLimit) .OR. (V9 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V10 .LE. LLimit) .OR. (V10 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V11 .LE. LLimit) .OR. (V11 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V12 .LE. LLimit) .OR. (V12 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V13 .LE. LLimit) .OR. (V13 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V14 .LE. LLimit) .OR. (V14 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V15 .LE. LLimit) .OR. (V15 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V16 .LE. LLimit) .OR. (V16 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V17 .LE. LLimit) .OR. (V17 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V18 .LE. LLimit) .OR. (V18 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V19 .LE. LLimit) .OR. (V19 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V20 .LE. LLimit) .OR. (V20 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V21 .LE. LLimit) .OR. (V21 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V22 .LE. LLimit) .OR. (V22 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V23 .LE. LLimit) .OR. (V23 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSEIF ((V24 .LE. LLimit) .OR. (V24 .GE. ULimit)) THEN
PRINT *, 'Voltage out of bounds!, *****'
ELSE
PRINT *, 'Within voltage limits.'
ENDIF
PRINT *, ' '
PRINT *, '***** End of cycle *****'
PRINT *, ' '
c Adjust separator area & go back into loop.
Separator = Separator - SepStep
ENDDO
c End of separator size loop, second loop.
c Adjust state of charge and go back into loop.
SocA = SocA - SocStep
ENDDO
c End of final loop, state of charge.
c Exit program.
STOP
END
c ************ Define functions and subroutines. ************
c Deterimes overpotential for given current, linear best fit.
REAL FUNCTION OVERPOTENTIAL (i)
REAL i
COMMON /op/ Op1, Op2, Op3
OVERPOTENTIAL = Op1*i**2 + Op2*i + Op3
RETURN
END
c **************************************************
c Determines open circuit potential of battery, Bellcore battery.
REAL FUNCTION OPEN (x, y)
REAL x, y
REAL PA, PC
COMMON /oc/ OcA1, OcA2, OcA3, OcC1, OcC2, OcC3
PA = OcA1*x**2 + OcA2*x + OcA3
PC = OcC1*y**2 + OcC2*y + OcC3
OPEN = PC - PA
RETURN
END
c **************************************************
c Determines state of charge adjustment for overpotential.
c Different overpotentials at different states of charge, slight correction.
REAL FUNCTION TUNE (x, y)
REAL x, y
REAL A1, A2
COMMON /ad/ ScA1, ScA2, ScA3, ScC1, ScC2, ScC3
A1 = ScA1*x**2 + ScA2*x + ScA3
A2 = ScC1*y**2 + ScC2*y + ScC3
TUNE = (A1)*(A2)
RETURN
END
c **************************************************
c Determines new state of charge after current segment, Faraday's law.
SUBROUTINE CHARGE (c, x, y, t, ma, mc)
REAL c, x, y, ma, mc
INTEGER t
REAL MolPA, MolPC
INTEGER FARADY
PARAMETER (FARADY = 96487)
MolPA = ma*x
MolPC = mc*y
MolPA = MolPA - c*t/FARADY
MolPC = MolPC + c*t/FARADY
c New states of charge of electrodes.
x = MolPA/ma
y = MolPC/mc
RETURN
END
c **************************************************
c Determines current for each driving segment.
SUBROUTINE CURRENT (p, c, v, x, y, a, t, ma, mc)
REAL p, c, v, x, y, a, ma, mc
INTEGER t
REAL CHANGE, OVER, FIX1, FIX2, FIX
REAL BatPower1, Batpower2, i1, i2, Error1, Error2
REAL OCP, xi, yi
INTEGER Count
c Guess initial current (P=I*V), Open circuit for initial value.
OCP = OPEN (x,y)
c = p/OCP
Count = 0
Change = 1
xi = x
yi = y
c FIX1 doesn't need to recalculated since based on initial state of charge.
FIX1 = TUNE (x, y)
c All of this stuff is just for an initial point to begin the iteration.
OVER = OVERPOTENTIAL (c)
CALL CHARGE (c, x, y, t, ma, mc)
FIX2 = TUNE (x, y)
FIX = (FIX1 + FIX2)/2
OVER = OVER*FIX
v = OCP - OVER
BatPower1 = c*v
Error1 = BatPower1 - p
c Set i1 = c for initial iteration.
i1 = c
c Reset x & y to initial values to begin loop.
x = xi
y = yi
c New guess for second current value in initial iteration.
c = c - Error1/v
DO WHILE (Change .NE. 0)
c Add counting section to get out if not converging.
Count = Count + 1
IF (Count .GE. 100) THEN
PRINT *, 'I have been in this current loop too long!'
STOP
ENDIF
OVER = OVERPOTENTIAL (c)
CALL CHARGE (c, x, y, t, ma, mc)
FIX2 = TUNE (x, y)
FIX = (FIX1 + FIX2)/2
OVER = OVER*FIX
v = OCP - OVER
BatPower2 = c*v
Error2 = BatPower2 - p
IF (ABS(Error2) .GE. 0.001) THEN
c Need i2 to hold value of c for i1.
i2 = c
c Linear interpolation to guess new current guess.
c = ((i2 - i1)/(Error2 - Error1))*(0 - Error1) + i1
x = xi
y = yi
c Change second variables to first variables for next iteration.
Error1 = Error2
i1 = i2
ELSE
Change = 0
ENDIF
ENDDO
c Print out details after loop has converged. Will give too many details.
c print *, p
c print *, 'Cell Pot (volts), i (A/m2), P (w/m2), x, y'
c print *, v,',',c,',',v*c,',',x,',',y
RETURN
END
c **************************************************
c Calculates fuel efficiency
SUBROUTINE FUEL (e, m)
REAL e, m
REAL Vh, Hpower, Heng, Hkml, Udist, Ukml, Ekml
REAL Eave1a, Eave1b, Eave1, Eave2, Eave
REAL AIR, GRAV
INTEGER EGAS, t1, t2, t3, t4, t5, t6, tb, ts
PARAMETER (EGAS = 31810, AIR = 1.202, GRAV = 9.81)
COMMON /gas/ C1, C2, Area, Peff, Eeff, Geff, Veng
c Driving cycle time segment lengths in seconds.
t1 = 10
t2 = 30
t3 = 12
t4 = 28
t5 = 15
t6 = 25
tb = 5
ts = 15
c Highway cruising speed = 27.778 m/s
Vh = 27.778
c No conversion factor required to get in km/L
Hpower = C1*GRAV*Vh*m + 0.5*AIR*C2*Area*Vh**3
Heng = Hpower/Peff
Hkml = (Vh/Heng)*EGAS*Eeff
Hbat = Hpower/Geff
c Urban distance travelled ~ 3535 meters in 6 minutes
Udist = 3535
c Calculate average urban engine since variable design possible.
Eave1a = (1+Veng)*(3*t1+2*t3+t5)+(3*t2+2*t4+t6)
Eave1b = (1-Veng)*6*(tb+ts)
Eave1 = (Eave1a + Eave1b)*e
Eave2 = 3*(t1+t2)+2*(t3+t4)+t5+t6+6*(tb+ts)
Eave = Eave1/Eave2
Ukml = (Udist/(360*Eave))*EGAS*Eeff
c Equivalent miles per gallon: 55% urban driving & 45% highway driving
Ekml = 0.55*Ukml + 0.45*Hkml
c Print out mileage estimates
PRINT *, 'Engine variability =,', Veng
PRINT *, 'Average urban engine (watts) =,', Eave
PRINT *, 'Highway cruising engine (watts) =,', Heng
PRINT *, 'ZEV battery power (watts) =,', Hbat
PRINT *, 'Highway mileage (km/L) =,', Hkml
PRINT *, 'Urban mileage (km/L) =,', Ukml
PRINT *, 'Average mileage (km/L) =,', Ekml
c Multiple by 0.6214/0.26417 = 2.35227 to convert from km/L to mpg
PRINT *, 'Highway mileage (mpg) =,', Hkml*2.35227
PRINT *, 'Urban mileage (mpg) =,', Ukml*2.35227
PRINT *, 'Average mileage (mpg) =,', Ekml*2.35227
RETURN
END
c **************************************************
c Determines approximate maximum battery power available. (10 seconds)
SUBROUTINE POWER (x, y, a, s, l)
REAL x, y, a, s, l
REAL C, Pnew, Pold, OCP, LESS, VOLT
c Initialize variables and get into loop.
C = 0
OCP = OPEN (x, y)
Pnew = 5
Pold = 0
VOLT = 4
c Loop to determine maximum power of battery with voltage restriction.
DO WHILE ((Pnew .GT. Pold) .AND. (VOLT .GE. l))
Pold = Pnew
C = C + 1
LESS = OVERPOTENTIAL (C)
VOLT = OCP - LESS
Pnew = C*VOLT*s
ENDDO
PRINT *, ' '
PRINT *, 'Voltage at maximum power (volts) =,', VOLT
PRINT *, 'Current at maximum power (A/m2) =,', C
PRINT *, 'Maximum 10 second power (watts) =,', Pnew
PRINT *, ' '
c Loop to determine maximum power of battery w/o voltage restriction.
DO WHILE (Pnew .GT. Pold)
Pold = Pnew
C = C + 1
LESS = OVERPOTENTIAL (C)
VOLT = OCP - LESS
Pnew = C*VOLT*s
ENDDO
PRINT *, 'Voltage at maximum power (volts) =,', VOLT
PRINT *, 'Current at maximum power (A/m2) =,', C
PRINT *, 'Maximum 10 second power (watts) =,', Pnew
PRINT *, ' '
RETURN
END