Line data Source code
1 : /****************************************************************************/
2 : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3 : // Copyright (C) 2001-2024 German Aerospace Center (DLR) and others.
4 : // This program and the accompanying materials are made available under the
5 : // terms of the Eclipse Public License 2.0 which is available at
6 : // https://www.eclipse.org/legal/epl-2.0/
7 : // This Source Code may also be made available under the following Secondary
8 : // Licenses when the conditions for such availability set forth in the Eclipse
9 : // Public License 2.0 are satisfied: GNU General Public License, version 2
10 : // or later which is available at
11 : // https://www.gnu.org/licenses/old-licenses/gpl-2.0-standalone.html
12 : // SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
13 : /****************************************************************************/
14 : /// @file HelpersEnergy.cpp
15 : /// @author Daniel Krajzewicz
16 : /// @author Jakob Erdmann
17 : /// @author Michael Behrisch
18 : /// @date Mon, 10.05.2004
19 : ///
20 : // Helper methods for HBEFA-based emission computation
21 : /****************************************************************************/
22 : #include <config.h>
23 :
24 : #include <complex> // due to sqrt of complex
25 : #include <utils/common/MsgHandler.h> // due to WRITE_WARNING
26 : #include <utils/common/SUMOTime.h>
27 : #include <utils/common/ToString.h>
28 : #include <utils/common/PolySolver.h>
29 : #include "HelpersEnergy.h"
30 :
31 :
32 : // ===========================================================================
33 : // method definitions
34 : // ===========================================================================
35 56178 : HelpersEnergy::HelpersEnergy():
36 56178 : PollutantsInterface::Helper("Energy", ENERGY_BASE, ENERGY_BASE)
37 56178 : { }
38 :
39 :
40 : double
41 2330524 : HelpersEnergy::compute(const SUMOEmissionClass /* c */, const PollutantsInterface::EmissionType e, const double v, const double a, const double slope, const EnergyParams* param) const {
42 2330524 : if (e != PollutantsInterface::ELEC) {
43 : return 0.;
44 : }
45 655948 : if (param == nullptr) {
46 0 : param = EnergyParams::getDefault();
47 : }
48 : //@ToDo: All formulas below work with the logic of the euler update (refs #860).
49 : // Approximation order could be improved. Refs. #2592.
50 :
51 655948 : const double lastV = v - ACCEL2SPEED(a);
52 655948 : const double mass = param->getTotalMass(myDefaultMass, 0.);
53 :
54 : // calculate power needed for potential energy difference
55 655948 : double power = mass * GRAVITY * sin(DEG2RAD(slope)) * v;
56 :
57 : // power needed for kinetic energy difference of vehicle
58 655948 : power += 0.5 * mass * (v * v - lastV * lastV) / TS;
59 :
60 : // add power needed for rotational energy diff of internal rotating elements
61 655948 : power += 0.5 * param->getDoubleOptional(SUMO_ATTR_ROTATINGMASS, myDefaultRotatingMass) * (v * v - lastV * lastV) / TS;
62 :
63 : // power needed for Energy loss through Air resistance [Ws]
64 : // Calculate energy losses:
65 : // EnergyLoss,Air = 1/2 * rho_air [kg/m^3] * myFrontSurfaceArea [m^2] * myAirDragCoefficient [-] * v_Veh^2 [m/s] * s [m]
66 : // ... with rho_air [kg/m^3] = 1,2041 kg/m^3 (at T = 20C)
67 : // ... with s [m] = v_Veh [m/s] * TS [s]
68 : // divided by TS to get power instead of energy
69 655948 : power += 0.5 * 1.2041 * param->getDoubleOptional(SUMO_ATTR_FRONTSURFACEAREA, myDefaultFrontSurfaceArea) * param->getDoubleOptional(SUMO_ATTR_AIRDRAGCOEFFICIENT, myDefaultAirDragCoefficient) * v * v * v;
70 :
71 : // power needed for Energy loss through Roll resistance [Ws]
72 : // ... (fabs(veh.getSpeed())>=0.01) = 0, if vehicle isn't moving
73 : // EnergyLoss,Tire = c_R [-] * F_N [N] * s [m]
74 : // ... with c_R = ~0.012 (car tire on asphalt)
75 : // ... with F_N [N] = myMass [kg] * g [m/s^2]
76 : // divided by TS to get power instead of energy
77 655948 : power += param->getDoubleOptional(SUMO_ATTR_ROLLDRAGCOEFFICIENT, myDefaultRollDragCoefficient) * GRAVITY * mass * v;
78 :
79 : // Energy loss through friction by radial force [Ws]
80 : // If angle of vehicle was changed
81 655948 : const double angleDiff = param->getAngleDiff();
82 655948 : if (angleDiff != 0.) {
83 : // Compute new radio
84 5500 : double radius = SPEED2DIST(v) / fabs(angleDiff);
85 :
86 : // Check if radius is in the interval [0.0001 - 10000] (To avoid overflow and division by zero)
87 5500 : if (radius < 0.0001) {
88 : radius = 0.0001;
89 5347 : } else if (radius > 10000) {
90 : radius = 10000;
91 : }
92 : // EnergyLoss,internalFrictionRadialForce = c [m] * F_rad [N];
93 : // Energy loss through friction by radial force [Ws]
94 : // divided by TS to get power instead of energy
95 5500 : power += param->getDoubleOptional(SUMO_ATTR_RADIALDRAGCOEFFICIENT, myDefaultRadialDragCoefficient) * mass * v * v / radius / TS;
96 : }
97 :
98 : // EnergyLoss,constantConsumers
99 : // Energy loss through constant loads (e.g. A/C) [W]
100 655948 : power += param->getDoubleOptional(SUMO_ATTR_CONSTANTPOWERINTAKE, myDefaultConstantPowerIntake);
101 :
102 : // E_Bat = E_kin_pot + EnergyLoss;
103 655948 : if (power > 0) {
104 : // Assumption: Efficiency of myPropulsionEfficiency when accelerating
105 601082 : power /= param->getDoubleOptional(SUMO_ATTR_PROPULSIONEFFICIENCY, myDefaultPropulsionEfficiency);
106 : } else {
107 : // Assumption: Efficiency of myRecuperationEfficiency when recuperating
108 54866 : power *= param->getDoubleOptional(SUMO_ATTR_RECUPERATIONEFFICIENCY, myDefaultRecuperationEfficiency);
109 54866 : if (a != 0) {
110 : // Fiori, Chiara & Ahn, Kyoungho & Rakha, Hesham. (2016).
111 : // Power-based electric vehicle energy consumption model: Model
112 : // development and validation. Applied Energy. 168. 257-268.
113 : // 10.1016/j.apenergy.2016.01.097.
114 : //
115 : // Insaf Sagaama, Amine Kchiche, Wassim Trojet, Farouk Kamoun
116 : // Improving The Accuracy of The Energy Consumption Model for
117 : // Electric Vehicle in SUMO Considering The Ambient Temperature
118 : // Effects
119 19542 : power *= (1 / exp(param->getDoubleOptional(SUMO_ATTR_RECUPERATIONEFFICIENCY_BY_DECELERATION, myDefaultRecuperationEfficiencyByDeceleration) / fabs(a)));
120 : }
121 : }
122 :
123 : // convert from [W], to [Wh/s] (3600s / 1h):
124 655948 : return power / 3600.;
125 : }
126 :
127 :
128 : double
129 70 : HelpersEnergy::acceleration(const SUMOEmissionClass /* c */, const PollutantsInterface::EmissionType e, const double v, const double P, const double slope, const EnergyParams* param) const {
130 70 : if (e != PollutantsInterface::ELEC) {
131 : return 0.;
132 : }
133 :
134 70 : if (param == nullptr) {
135 0 : param = EnergyParams::getDefault();
136 : }
137 :
138 : // Inverse formula for the function compute()
139 :
140 : // Computes the achievable acceleration using the given speed and drive power in [Wh/s]
141 : // It does not consider friction losses by radial force and the acceleration-dependent
142 : // recuperation efficiency (eta_recuperation is considered constant)
143 : //
144 : // The power `Prest` used for acceleration is given by a cubic polynomial in acceleration,
145 : // i.e.the equation
146 : //
147 : // Prest = const1*a + const2*a^2 + const3*a^3
148 : //
149 : // and we need to find `a` given `Prest`.
150 : //
151 : // The solutions of `a(P)` is stored in variables `x1`, `x2`, and `x3` returned by
152 : // the method `PolySolver::cubicSolve()`, see below.
153 : //
154 : // Power used for accelerating, `Prest`, is the total used power minus power wasted by running resistances.
155 :
156 70 : const double mass = param->getTotalMass(myDefaultMass, 0.);
157 70 : const double rotMass = param->getDoubleOptional(SUMO_ATTR_ROTATINGMASS, myDefaultRotatingMass);
158 : double const1, const2, const3;
159 : double Prest;
160 : int numX;
161 : double x1, x2, x3;
162 :
163 : // Power delivered by the drive
164 70 : if (P > 0) {
165 : // Assumption: Efficiency of myPropulsionEfficiency when accelerating
166 70 : Prest = P * 3600 * param->getDoubleOptional(SUMO_ATTR_PROPULSIONEFFICIENCY, myDefaultPropulsionEfficiency);
167 : } else {
168 : // Assumption: Efficiency of myRecuperationEfficiency when recuperating
169 0 : Prest = P * 3600 / param->getDoubleOptional(SUMO_ATTR_RECUPERATIONEFFICIENCY, myDefaultRecuperationEfficiency);
170 : }
171 :
172 : // calculate power drop due to a potential energy difference
173 : // in inverse original formula: power = mass * GRAVITY * sin(DEG2RAD(slope)) * v;
174 : // i.e. in terms of 'lastV' and 'a': power = mass * GRAVITY * sin(DEG2RAD(slope)) * (lastV + TS * a);
175 70 : Prest -= mass * GRAVITY * sin(DEG2RAD(slope)) * (v);
176 70 : const1 = mass * GRAVITY * sin(DEG2RAD(slope)) * (TS);
177 :
178 : // update coefficients of a(P) equation considering power loss through Roll resistance
179 : // in inverse original formula: power += param->getDoubleOptional(SUMO_ATTR_ROLLDRAGCOEFFICIENT) * GRAVITY * mass * v;
180 : // i.e. in terms of 'lastV' and 'a': power += param->getDoubleOptional(SUMO_ATTR_ROLLDRAGCOEFFICIENT) * GRAVITY * mass * (lastV + TS * a);
181 70 : Prest -= param->getDoubleOptional(SUMO_ATTR_ROLLDRAGCOEFFICIENT, myDefaultRollDragCoefficient) * GRAVITY * mass * v;
182 70 : const1 += param->getDoubleOptional(SUMO_ATTR_ROLLDRAGCOEFFICIENT, myDefaultRollDragCoefficient) * GRAVITY * mass * (TS);
183 :
184 : //Constant loads are omitted. We assume P as the max limit for the main traction drive. Constant loads are often covered by an auxiliary drive
185 : //Power loss through constant loads (e.g. A/C) [W]
186 : //Prest -= param->getDoubleOptional(SUMO_ATTR_CONSTANTPOWERINTAKE) / TS;
187 :
188 : // update coefficients of a(P) equation considering kinetic energy difference of vehicle
189 : // in inverse original formula: power += 0.5 * mass * (v * v - lastV * lastV) / TS;
190 : // i.e. in terms of 'lastV' and 'a': power += 0.5 * mass * (2 * lastV * a + TS * a * a);
191 70 : const1 += 0.5 * mass * (2 * v);
192 70 : const2 = 0.5 * mass * (TS);
193 :
194 : // update coefficients of a(P) equation considering rotational energy diff of internal rotating elements
195 : // in inverse original formula: power += 0.5 * rotMass * (v * v - lastV * lastV) / TS;
196 : // i.e. in terms of 'lastV' and 'a': power += 0.5 * rotMass * (2 * lastV * a + TS * a * a);
197 70 : const1 += 0.5 * rotMass * (2 * v);
198 70 : const2 += 0.5 * rotMass * (TS);
199 :
200 : // update coefficients of a(P) equation considering energy loss through Air resistance
201 : // in inverse original formula: power += 0.5 * 1.2041 * param->getDoubleOptional(SUMO_ATTR_FRONTSURFACEAREA) * param->getDoubleOptional(SUMO_ATTR_AIRDRAGCOEFFICIENT) * v * v * v;
202 : // i.e. in terms of 'lastV' and 'a': power += 0.5 * rotMass * (lastV^3 + 3* lastV^2 * TS * a + 3* lastV * TS^2 *a^2 + TS^3 * a^3);
203 70 : const double drag = 0.5 * 1.2041 * param->getDoubleOptional(SUMO_ATTR_FRONTSURFACEAREA, myDefaultFrontSurfaceArea) * param->getDoubleOptional(SUMO_ATTR_AIRDRAGCOEFFICIENT, myDefaultAirDragCoefficient);
204 70 : Prest -= drag * (v * v * v);
205 70 : const1 += drag * (3 * v * v * TS);
206 70 : const2 += drag * (3 * v * TS * TS);
207 70 : const3 = drag * (TS * TS * TS);
208 :
209 : // Prest = const1*a + const2*a^2 + const3*a^3
210 : // Solve cubic equation in `a` for real roots, and return the number of roots in `numX`
211 : // and the roots in `x1`, `x2`, and `x3`.
212 70 : std::tie(numX, x1, x2, x3) = PolySolver::cubicSolve(const3, const2, const1, -Prest);
213 :
214 :
215 70 : switch (numX) {
216 : case 1:
217 : return x1;
218 0 : case 2:
219 : return MAX2(x1, x2);
220 70 : case 3:
221 : return MAX3(x1, x2, x3);
222 0 : default:
223 0 : WRITE_ERROR(TL("An acceleration given by the power was not found."));
224 0 : return 0;
225 : }
226 :
227 : }
228 :
229 :
230 : /****************************************************************************/
|
