Mercurial > dive4elements > river

comparison artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/tkhcalculation/TkhCalculator.java @ 8920:29442c03c6e3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
d50 aggregation by median instead of arithmetic mean, negative tkh replaced by 0
author mschaefer
date Thu, 01 Mar 2018 09:01:58 +0100
parents 5d5d0051723f
children 9c02733a1b3c
comparison
equal deleted inserted replaced
8916:5d5d0051723f 8920:29442c03c6e3
56 final Integer soundingYear = bedHeightsProvider.getInfo().getYear(); 56 final Integer soundingYear = bedHeightsProvider.getInfo().getYear();
57 final BedQualityD50KmValueFinder bedMeasurementsFinder = BedQualityD50KmValueFinder.loadBedMeasurements(river, calcRange, soundingYear, 57 final BedQualityD50KmValueFinder bedMeasurementsFinder = BedQualityD50KmValueFinder.loadBedMeasurements(river, calcRange, soundingYear,
58 VALID_BED_MEASUREMENT_YEARS); 58 VALID_BED_MEASUREMENT_YEARS);
59 59
60 if (bedMeasurementsFinder == null) { 60 if (bedMeasurementsFinder == null) {
61 final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingSoilKind", null, label); 61 final String message = Resources.getMsg(context.getMeta(), "sinfo_calc_flow_depth.warning.missingD50", null, label);
62 problems.addProblem(message); 62 problems.addProblem(message);
63 return null; 63 return null;
64 } 64 }
65 65
66 // FIXME: wie wird ggf. interpoliert? prüfung ob werte vorhanden? 66 // FIXME: wie wird ggf. interpoliert? prüfung ob werte vorhanden?
157 if (Double.isNaN(d50)) 157 if (Double.isNaN(d50))
158 return new Tkh(km, wst, meanBedHeight, discharge, kind, Double.NaN, Double.NaN, Double.NaN); 158 return new Tkh(km, wst, meanBedHeight, discharge, kind, Double.NaN, Double.NaN, Double.NaN);
159 159
160 if (!this.flowVelocitiesFinder.findKmQValues(km, discharge)) { 160 if (!this.flowVelocitiesFinder.findKmQValues(km, discharge)) {
161 // TODO: ggf. station in Fehlermeldung? 161 // TODO: ggf. station in Fehlermeldung?
162 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, this.problemLabel); 162 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingVelocity", null, this.problemLabel);
163 this.problems.addProblem(km, message); 163 this.problems.addProblem(km, message);
164 // FIXME: cumulate problems to one message? 164 // FIXME: cumulate problems to one message?
165 return new Tkh(km, wst, meanBedHeight, discharge, kind, Double.NaN, Double.NaN, Double.NaN);
165 } 166 }
166 167
167 final double tkh = calculateTkh(wst - meanBedHeight, this.flowVelocitiesFinder.getFindVmainFound(), d50, this.flowVelocitiesFinder.getFindTauFound()); 168 final double tkh = calculateTkh(wst - meanBedHeight, this.flowVelocitiesFinder.getFindVmainFound(), d50, this.flowVelocitiesFinder.getFindTauFound());
168 // FIXME: noch mal prüfen, im alten code wurde hier immer auf 0 gesetzt
169 if (Double.isNaN(tkh) || (tkh < 0)) {
170 // TODO: ggf. station in Fehlermeldung?
171
172 // FIXME: Fehlermeldung nicht korrekt, passiert mit Wasserspiegel 'MHQ' und 'QP-1993
173 final String message = Resources.getMsg(this.context.getMeta(), "sinfo_calc_flow_depth.warning.missingQ", null, this.problemLabel);
174 this.problems.addProblem(km, message);
175
176 return new Tkh(km, wst, meanBedHeight, discharge, kind, Double.NaN, Double.NaN, Double.NaN);
177 }
178
179 /*
180 * log.debug(String.format("calculateTkh km %.3f q %.0f w %.2f mbh %.2f vm %.1f tau %.1f d50(mm) %.1f tkh(cm) %.1f",
181 * km, discharge, wst, meanBedHeight, flowVelocitiesFinder.getFindVmainFound(), flowVelocitiesFinder.getFindTauFound(),
182 * d50*1000, tkh));
183 */
184
185 double tkhUp; 169 double tkhUp;
186 double tkhDown; 170 double tkhDown;
187 switch (kind) { 171 switch (kind) {
188 case starr: 172 case starr:
189 tkhUp = tkh; 173 tkhUp = tkh;
209 * flow velocity in m 193 * flow velocity in m
210 * @param d50 194 * @param d50
211 * grain diameter D50 in m (!) 195 * grain diameter D50 in m (!)
212 * @param tau 196 * @param tau
213 * shear stress in N/m^2 197 * shear stress in N/m^2
214 * @return transport body height in cm (!) 198 * @return transport body height in cm (!), never negative
215 */ 199 */
216 private double calculateTkh(final double h, final double vm, final double d50, final double tau) { 200 private double calculateTkh(final double h, final double vm, final double d50, final double tau) {
217 final double PHYS_G = 9.81; 201 final double PHYS_G = 9.81;
218 final double PHYS_SPECGRAV_S = 2.6; 202 final double PHYS_SPECGRAV_S = 2.6;
219 final double PHYS_VELOCCOEFF_N = 6; 203 final double PHYS_VELOCCOEFF_N = 6;
224 208
225 final double froude = vm / Math.sqrt(PHYS_G * h); 209 final double froude = vm / Math.sqrt(PHYS_G * h);
226 final double partReynolds = Math.sqrt((PHYS_SPECGRAV_S - 1) * PHYS_G * d50) / PHYS_VISCOSITY_NUE * d50; 210 final double partReynolds = Math.sqrt((PHYS_SPECGRAV_S - 1) * PHYS_G * d50) / PHYS_VISCOSITY_NUE * d50;
227 final double critShields = 0.22 * Math.pow(partReynolds, -0.6) + 0.06 * Math.pow(10, 7.7 * Math.pow(partReynolds, -0.6)); 211 final double critShields = 0.22 * Math.pow(partReynolds, -0.6) + 0.06 * Math.pow(10, 7.7 * Math.pow(partReynolds, -0.6));
228 final double critTau = critShields * (PHYS_GRAIN_DENSITY_RHOS - PHYS_WATER_DENSITY_RHO) * PHYS_G * d50; 212 final double critTau = critShields * (PHYS_GRAIN_DENSITY_RHOS - PHYS_WATER_DENSITY_RHO) * PHYS_G * d50;
229 return 100 * h * (1 - Math.pow(froude, 2)) / (2 * PHYS_VELOCCOEFF_N * PHYS_FORMCOEFF_ALPHA) * (1 - critTau / tau); 213 final double tkh = 100 * h * (1 - Math.pow(froude, 2)) / (2 * PHYS_VELOCCOEFF_N * PHYS_FORMCOEFF_ALPHA) * (1 - critTau / tau);
214 // Some regular input values may give a negative calculation result; that is unwanted
215 if (tkh < 0.0)
216 return 0.0;
217 else
218 return tkh;
230 } 219 }
231 } 220 }
http://dive4elements.wald.intevation.org