ChannelGeometry.java [src/java/hydraulics] Revision: ca784ea44b9b35192fe845b26152967df2964653 Date: Wed May 03 10:00:06 MDT 2017
package hydraulics;
import java.io.IOException;
import java.util.ArrayList;
/**
* Last Updated: 3-May-2017
* @author Tyler Wible
* @since 1-August-2014
*/
public class ChannelGeometry {
String depthType = "user";//"normal";//"critical";// //The depth to be used for the cross-section properties (if normal, Manning's equation is used, if 'user' the provided depth is used)
double depth = 5; //the depth of flow in the cross-section (meters)
boolean irregularGeometry = false; //Whether the cross-section is regular like trapezoidal, triangular, rectangular (false) or irregular (true)
double bottomWidth = 6.096; //the bottom width of the channel (set to zero for triangular cross-sections, meters)
double sideSlope = 2; //the side slope (z), calculated with a right triangle with a vertical unit of 1 and a horizontal unit of "z"
public double discharge = 11.3267; //discharge in the cross-section (in cubic meters per second, m^3/s) (only used for normal and critical depth calculations)
double bedSlope = 0.0016; //channel bed slope (m/m) (only used for normal depth calculations)
public double n = 0.025; //manning's roughness coefficient (only used for normal depth calculations)
String crossSection = "x\ty\n0\t3250\n10\t3240\n41.4\t3242\n52.3\t3243\n61.8\t3248\n72\t3241\n84\t3239\n87\t3242\n95\t3248\n123\t3252";//a list of x-y points for the cross section ;
//Global variable (no set needed)
double[][] xyPoints = {{0,368}, //a list of x-y points for the cross section, is filled when "setCrossSection" is called.
{5,350}, //xyPoints[i][0] = x points (measured from left to right starting at zero, in meters),
{10,350}, //xyPoints[i][1] = y points (total elevation of the ground surface, in meters)
{12,350}};
//Outputs
double area = -1;
double perimeter = -1;
double topWidth = -1;
double hydraulicRadius = -1;
double hydraulicDepth = -1;
double normalDischarge = -1;
String warningMessage = "?";
//Gets
public double getDepth(){
return depth;
}
public double getArea(){
return area;
}
public double getWettedPerimeter(){
return perimeter;
}
public double getTopWidth(){
return topWidth;
}
public double getHydraulicRadius(){
return hydraulicRadius;
}
public double getHydraulicDepth(){
return hydraulicDepth;
}
public double getNormalDischarge(){
return normalDischarge;
}
public String getWarningMessage(){
return warningMessage;
}
//Sets
public void setDepthType(String depthType){
this.depthType = depthType;
}
public void setIrregularGeometryTF(boolean irregularGeometry){
this.irregularGeometry = irregularGeometry;
}
public void setDepth(double depth){
this.depth = depth;
}
public void setBottomWidth(double bottomWidth){
this.bottomWidth = bottomWidth;
}
public void setSideSlope(double sideSlope){
this.sideSlope = sideSlope;
}
public void setDischarge(double discharge){
this.discharge = discharge;
}
public void setBedSlope(double bedSlope){
this.bedSlope = bedSlope;
}
public void setManningRoughness(double n){
this.n = n;
}
public void setCrossSection(String crossSection) throws IOException{
this.crossSection = crossSection;
//Break up formatted cross-section points and create a double array for easier use
String[] xyPointList = crossSection.split("\n");
if(xyPointList.length > 1){
//If the list of points is of sufficient size and format, store it
double[][] xyPoints_new = new double[xyPointList.length - 1][2];
for(int i=1; i<xyPointList.length; i++){
String[] xy = xyPointList[i].split("\t");
xyPoints_new[i-1][0] = Double.parseDouble(xy[0]);
xyPoints_new[i-1][1] = Double.parseDouble(xy[1]);
}
this.xyPoints = xyPoints_new;
}else{
//if the list is empty or formatted wrong, thrown an error
ArrayList<String> errorMessage = new ArrayList<String>();
errorMessage.add("The list of XY points provided is not in the expected tab-delimited format like: 'x\ty\n0\t350\n5\t368\n10\t376\n12\t350\n'");
writeError(errorMessage);
}
}
public void setXYPoints(double[][] xyPoints) throws IOException{
this.xyPoints = xyPoints;
String crossSection = "";
for(int i=0; i<xyPoints.length; i++){
if(i==0){
crossSection = "x\ty\n" + xyPoints[i][0] + "\t" + xyPoints[i][1];
}else{
crossSection = crossSection + "\n" + xyPoints[i][0] + "\t" + xyPoints[i][1];
}
}
this.crossSection = crossSection;
}
/**
* Writes out the error message, if any, and then exits the program
* @param error string array to be written as each line of an error message
* @throws IOException
*/
public void writeError(ArrayList<String> error) throws IOException{
//Output data to text file
String errorContents = error.get(0);
for(int i=1; i<error.size(); i++){
errorContents = errorContents + "\n" + error.get(i);
}
throw new IOException("Error encountered. Please see the following message for details: \n" + errorContents);
}
/**
* Calculates the area of a trapezoidal, triangular with bottomWidth = 0, or
* rectangular with sideSlope = 0 channel cross section.
*
* Note: this assumes the same units for depth and bottom width resulting
* in an area with those units squared
* @return the area of the trapezoidal cross-section
*/
public double calcTrapezoidalArea(){
this.area = (bottomWidth + sideSlope * depth) * depth;
return area;
}
/**
* Calculates the wetted perimeter of a trapezoidal, triangular with
* bottomWidth = 0, or rectangular with sideSlope = 0 channel cross section.
*
* Note: this assumes the same units for depth and bottom width resulting
* in an wetted perimeter in those same units
* @return the wetted perimeter of the trapezoidal cross-section
*/
public double calcTrapezoidalPerimeter(){
this.perimeter = bottomWidth + 2 * depth * Math.sqrt(1 + (sideSlope * sideSlope));
return perimeter;
}
/**
* Calculates the top width of a trapezoidal, triangular with
* bottomWidth = 0, or rectangular with sideSlope = 0 channel cross section.
*
* Note this assumes the same units for depth and bottom width resulting
* in a top width in those same units
* @return the top width of the trapezoidal cross-section
*/
public double calcTrapezoidalTopWidth(){
this.topWidth = bottomWidth + 2 * sideSlope * depth;
return topWidth;
}
/**
* Calculates the hydraulic radius of a trapezoidal, triangular with
* bottomWidth = 0, or rectangular with sideSlope = 0 channel cross section.
*
* Note: this assumes the same units for depth and bottom width resulting
* in a hydraulic radius in those same units
* @return the hydraulic radius of the trapezoidal cross-section
*/
public double calcTrapezoidalHydraulicRadius(){
this.area = calcTrapezoidalArea();
this.perimeter = calcTrapezoidalPerimeter();
this.hydraulicRadius = area / perimeter;
return hydraulicRadius;
}
/**
* Calculates the hydraulic depth of a trapezoidal, triangular with
* bottomWidth = 0, or rectangular with sideSlope = 0 channel cross section.
*
* Note: this assumes the same units for depth and bottom width resulting
* in a hydraulic depth in those same units
* @return the hydraulic depth of the trapezoidal cross-section
*/
public double calcTrapezoidalHydraulicDepth(){
this.area = calcTrapezoidalArea();
this.topWidth = calcTrapezoidalTopWidth();
this.hydraulicDepth = area / topWidth;
return hydraulicDepth;
}
/**
* Calculates a power relationship = a * X ^ b which can be used for things
* ex. a stage-discharge relationship would be setup like this:
* discharge = a * stage ^ b
* @param a the coefficient multiplying the whole term
* @param b the exponent that the know parameter is raised to
* @param x the know parameter to be used in the power relationship (ex. stage)
* @return the result of the power relationship calculation (ex. discharge)
*/
public double calcPowerRelation(double a, double b, double x){
double value = a * Math.pow(x, b);
return value;
}
/**
* Calculates the area of an irregular cross-section via the list of
* provided x (perpendicular to flow), y (elevation) points simulating the
* ground surface
*
* Note: if the water elevation is higher than the edge points provided for
* the cross-section (the starting left point and ending right point) this
* method assumes a vertical edge wall and calculates a trapezoidal area for
* that section
*
* Note: that this assumes the same units for depth and bottom width
* resulting in an area with those units squared
* @return the area of the irregular cross-section
*/
public double calcIrregularArea(){
GeneralFunctions general = new GeneralFunctions();
double thalwegElevation = general.min(general.getColumn(xyPoints, 1));
double waterElevation = thalwegElevation + depth;
//Loop through and determine the area of flow in between each set of points moving from left (i=0) to right (i=n) across the cross-section
double x1 = 0, y1 = 0;
double newArea = 0;
for(int i=0; i<xyPoints.length; i++){
if(i > 0){
//Get current set of points (to the right of the last set)
double x2 = xyPoints[i][0];
double y2 = xyPoints[i][1];
//Check if water elevation is above the ground surface here
if(waterElevation > y1 && waterElevation > y2){
//Calculate a trapezoidal area between the two points
double width = x2 - x1;
double height1 = waterElevation - y1;
double height2 = waterElevation - y2;
newArea = newArea + width * ((height1 + height2) / 2);
}else if(waterElevation > y1){
//Calculate a triangular area on the right side
double[] xArray = {x1, x2};
double[] yArray = {y1, y2};
double xValue = general.linearInterpolation(yArray, xArray, waterElevation);//Note the x and y arrays are switched since I need to interpolate an "x" value on the cross-section which is a "y" value in the interpolation tool
double width = xValue - x1;
double height = waterElevation - y1;
newArea = newArea + 0.5 * width * height;
}else if(waterElevation > y2){
//Calculate a triangular area on the left side
double[] xArray = {x1, x2};
double[] yArray = {y1, y2};
double xValue = general.linearInterpolation(yArray, xArray, waterElevation);//Note the x and y arrays are switched since I need to interpolate an "x" value on the cross-section which is a "y" value in the interpolation tool
double width = x2 - xValue;
double height = waterElevation - y2;
newArea = newArea + 0.5 * width * height;
}
if ( (waterElevation > y1 && i == 1) || (waterElevation > y2 && i == xyPoints.length-1) ){
this.warningMessage = "The depth of flow exceeds one or both edges of this cross-section. "
+ "A vertical wall has been assumed at the edge(s) of the cross-section for this calculation. "
+ "More points should be added to the cross-section, and these hydraulic properties should be re-calculated if this assumption is not met.";
}
}
//Get last set of points
x1 = xyPoints[i][0];
y1 = xyPoints[i][1];
}
this.area = newArea;
return area;
}
/**
* Calculates the wetted perimeter of an irregular cross-section via the list of
* provided x (perpendicular to flow), y (elevation) points simulating the
* ground surface
*
* Note: if the water elevation is higher than the edge points provided for
* the cross-section (the starting left point and ending right point) this
* method assumes a vertical edge wall and calculates wetted perimeter
* accordingly
*
* Note: that this assumes the same units for depth and x-y width-elevations
* resulting in a perimeter with those units
* @return the wetted perimeter of the irregular cross-section
*/
public double calcIrregularPerimeter(){
GeneralFunctions general = new GeneralFunctions();
double thalwegElevation = general.min(general.getColumn(xyPoints, 1));
double waterElevation = thalwegElevation + depth;
double xMax = general.max(general.getColumn(xyPoints, 0));
//Loop through and determine the wetted perimeter of flow in between each set of points moving from left (i=0) to right (i=n) across the cross-section
double x1 = 0, y1 = 0;
double length = 0;
for(int i=0; i<xyPoints.length; i++){
if(i > 0){
//Get current set of points (to the right of the last set)
double x2 = xyPoints[i][0];
double y2 = xyPoints[i][1];
//Check if water elevation is above the ground surface here
if(waterElevation > y1 && waterElevation > y2){
//Calculate the line length between the two points
double height = y2 - y1;
double width = x2 - x1;
length = length + Math.sqrt(height*height + width*width);
//Check if this is an edge point of the overall cross-section
if(x1 == 0){
//Calculate the left vertical edge up to the water surface
length = length + (waterElevation - y1);
}else if(x2 == xMax){
//Calculate the right vertical edge up to the water surface
length = length + (waterElevation - y2);
}
}else if(waterElevation > y1){
//Calculate the line length from the last point up to the water surface
double[] xArray = {x1, x2};
double[] yArray = {y1, y2};
double xValue = general.linearInterpolation(yArray, xArray, waterElevation);//Note the x and y arrays are switched since I need to interpolate an "x" value on the cross-section which is a "y" value in the interpolation tool
double width = xValue - x1;
double height = waterElevation - y1;
length = length + Math.sqrt(height*height + width*width);
}else if(waterElevation > y2){
//Calculate the line length from water surface down to the first point
double[] xArray = {x1, x2};
double[] yArray = {y1, y2};
double xValue = general.linearInterpolation(yArray, xArray, waterElevation);//Note the x and y arrays are switched since I need to interpolate an "x" value on the cross-section which is a "y" value in the interpolation tool
double width = x2 - xValue;
double height = waterElevation - y2;
length = length + Math.sqrt(height*height + width*width);
}
if ( (waterElevation > y1 && i == 1) || (waterElevation > y2 && i == xyPoints.length-1) ){
this.warningMessage = "The depth of flow exceeds one or both edges of this cross-section. "
+ "A vertical wall has been assumed at the edge(s) of the cross-section for this calculation. "
+ "More points should be added to the cross-section, and these hydraulic properties should be re-calculated if this assumption is not met.";
}
}
//Get last set of points
x1 = xyPoints[i][0];
y1 = xyPoints[i][1];
}
this.perimeter = length;
return perimeter;
}
/**
* Calculates the top width of an irregular cross-section via the list of
* provided x (perpendicular to flow), y (elevation) points simulating the
* ground surface
*
* Note: if the water elevation is higher than the edge points provided for
* the cross-section (the starting left point and ending right point) this
* method assumes a vertical edge wall and calculates a top width which
* terminates at this edge.
*
* Note: that this assumes the same units for depth and x-y width-elevations
* resulting in a top width with those units
* @return the top width of the irregular cross-section
*/
public double calcIrregularTopWidth(){
GeneralFunctions general = new GeneralFunctions();
double thalwegElevation = general.min(general.getColumn(xyPoints, 1));
double waterElevation = thalwegElevation + depth;
//Loop through and determine the wetted perimeter of flow in between each set of points moving from left (i=0) to right (i=n) across the cross-section
double x1 = 0, y1 = 0;
double width = 0;
for(int i=0; i<xyPoints.length; i++){
if(i > 0){
//Get current set of points (to the right of the last set)
double x2 = xyPoints[i][0];
double y2 = xyPoints[i][1];
//Check if water elevation is above the ground surface here
if(waterElevation > y1 && waterElevation > y2){
//Calculate the width between the two points
width = width + x2 - x1;
}else if(waterElevation > y1){
//Calculate the width from the last point to the edge of the water on the right side
double[] xArray = {x1, x2};
double[] yArray = {y1, y2};
double xValue = general.linearInterpolation(yArray, xArray, waterElevation);//Note the x and y arrays are switched since I need to interpolate an "x" value on the cross-section which is a "y" value in the interpolation tool
width = width + xValue - x1;
}else if(waterElevation > y2){
//Calculate the width from the last point to the edge of the water on the left side
double[] xArray = {x1, x2};
double[] yArray = {y1, y2};
double xValue = general.linearInterpolation(yArray, xArray, waterElevation);//Note the x and y arrays are switched since I need to interpolate an "x" value on the cross-section which is a "y" value in the interpolation tool
width = width + x2 - xValue;
}
if ( (waterElevation > y1 && i == 1) || (waterElevation > y2 && i == xyPoints.length-1) ){
this.warningMessage = "The depth of flow exceeds one or both edges of this cross-section. "
+ "A vertical wall has been assumed at the edge(s) of the cross-section for this calculation. "
+ "More points should be added to the cross-section, and these hydraulic properties should be re-calculated if this assumption is not met.";
}
}
//Get last set of points
x1 = xyPoints[i][0];
y1 = xyPoints[i][1];
}
this.topWidth = width;
return topWidth;
}
/**
* Calculates the hydraulic radius of an irregular cross-section via the list of
* provided x (perpendicular to flow), y (elevation) points simulating the
* ground surface
*
* Note: if the water elevation is higher than the edge points provided for
* the cross-section (the starting left point and ending right point) this
* method assumes a vertical edge wall and calculates a trapezoidal area for
* that section
*
* Note: that this assumes the same units for depth and x-y width-elevations
* resulting in an area with those units squared
* @return the hydraulic radius of the irregular cross-section
*/
public double calcIrregularHydraulicRadius(){
//Get properties for hydraulic radius calculation
this.area = calcIrregularArea();
this.perimeter = calcIrregularPerimeter();
this.hydraulicRadius = area / perimeter;
return hydraulicRadius;
}
/**
* Calculates the hydraulic depth of an irregular cross-section via the list of
* provided x (perpendicular to flow), y (elevation) points simulating the
* ground surface
*
* Note: if the water elevation is higher than the edge points provided for
* the cross-section (the starting left point and ending right point) this
* method assumes a vertical edge wall and calculates a trapezoidal area for
* this section
*
* Note: that this assumes the same units for depth and x-y width-elevations
* resulting in a hydraulic depth with those units
* @return the hydraulic depth of the irregular cross-section
*/
public double calcIrregularHydraulicDepth(){
this.area = calcIrregularArea();
this.topWidth = calcIrregularTopWidth();
this.hydraulicDepth = area / topWidth;
return hydraulicDepth;
}
/**
* Calculates normal depth based on Manning's equation for a regular or
* irregular cross-section geometry based on the provided discharge and
* channel cross-section
* @return
* @throws IOException
*/
public double calcNormalDepth() throws IOException{
//Use bisection method to find the "normal" depth that results in the specified discharge
double a = 0.02;//minimum depth
double b = 999;//maximum depth
double error = 1;
int ctr = 0;
//Use bisection method to find new guess based on average of current limits
double guessDepth = 0;
while(Math.abs(error) > 0.00001 && ctr < 50){
//Estimate a new guess for depth
guessDepth = (a + b) / 2.;
if(irregularGeometry){
//Calculate discharge using Manning's equation (in metric) for a "irregular" cross-section (list of xy points)
ChannelGeometry geo = new ChannelGeometry();
geo.setIrregularGeometryTF(true);
geo.setDepth(guessDepth);
geo.setCrossSection(crossSection);
//Calculate discharge using Manning's equation (in metric)
double A = geo.calcIrregularArea();
double R = geo.calcIrregularHydraulicRadius();
double Q_manning = (1./n) * A * Math.pow(R, 2./3.) * Math.sqrt(bedSlope);
error = discharge - Q_manning;
}else{
//Calculate discharge using Manning's equation (in metric) for a "regular" cross-section (trapezoid, triangle, rectangle)
ChannelGeometry geo = new ChannelGeometry();
geo.setIrregularGeometryTF(false);
geo.setDepth(guessDepth);
geo.setBottomWidth(bottomWidth);
geo.setSideSlope(sideSlope);
//Calculate discharge using Manning's equation (in metric)
double A = geo.calcTrapezoidalArea();
double R = geo.calcTrapezoidalHydraulicRadius();
double Q_manning = (1./n) * A * Math.pow(R, 2./3.) * Math.sqrt(bedSlope);
error = discharge - Q_manning;
}
//Check if the result (discharge) based on the guessDepth is positive
//(like discharge - manning(a)) or negative (like discharge - manning(b))
//and reset the limits of the bisection method to interpolate between
//the two new points
if(error < 0){
b = guessDepth;
}else{
a = guessDepth;
}
ctr++;
}
return guessDepth;
}
/**
* Calculates normal discharge based on Manning's equation for a regular or
* irregular cross-section geometry based on the provided depth and
* **currently not used but useful for future development**
* channel cross-section
* @return the discharge for the provided depth and cross-section
* @throws IOException
*/
public double calcNormalDischarge() throws IOException{
double Q_manning = -1;
if(irregularGeometry){
//Calculate discharge using Manning's equation (in metric) for a "irregular" cross-section (list of xy points)
ChannelGeometry geo = new ChannelGeometry();
geo.setIrregularGeometryTF(true);
geo.setDepth(depth);
geo.setCrossSection(crossSection);
//Calculate discharge using Manning's equation (in metric)
double A = geo.calcIrregularArea();
double R = geo.calcIrregularHydraulicRadius();
Q_manning = (1./n) * A * Math.pow(R, 2./3.) * Math.sqrt(bedSlope);
}else{
//Calculate discharge using Manning's equation (in metric) for a "regular" cross-section (trapezoid, triangle, rectangle)
ChannelGeometry geo = new ChannelGeometry();
geo.setIrregularGeometryTF(false);
geo.setDepth(depth);
geo.setBottomWidth(bottomWidth);
geo.setSideSlope(sideSlope);
//Calculate discharge using Manning's equation (in metric)
double A = geo.calcTrapezoidalArea();
double R = geo.calcTrapezoidalHydraulicRadius();
Q_manning = (1./n) * A * Math.pow(R, 2./3.) * Math.sqrt(bedSlope);
}
return Q_manning;
}
/**
* Calculates critical depth based Froude Number set equal 1 (critical open
* channel flow) for a regular or irregular cross-section geometry based on
* the provided discharge and channel cross-section
* @return
* @throws IOException
*/
public double calcCriticalDepth() throws IOException{
//Use bisection method to find the "critical" depth that results in the specified discharge
double a = 0.02;//minimum depth
double b = 999;//maximum depth
double error = 1;
int ctr = 0;
//Use bisection method to find new guess based on average of current limits
double guessDepth = 0;
while(Math.abs(error) > 0.00001 && ctr < 50){
//Estimate a new guess for depth
guessDepth = (a + b) / 2.;
if(irregularGeometry){
//Calculate Froude number (in metric) for a "irregular" cross-section (list of xy points)
ChannelGeometry geo = new ChannelGeometry();
geo.setIrregularGeometryTF(true);
geo.setDepth(guessDepth);
geo.setCrossSection(crossSection);
//Calculate Froude number (in metric)
double A = geo.calcIrregularArea();
double hydDepth = geo.calcIrregularHydraulicDepth();
double fr = (discharge / A) / Math.sqrt(9.81 * hydDepth);
error = 1 - fr;
}else{
//Calculate Froude number (in metric) for a "regular" cross-section (trapezoid, triangle, rectangle)
ChannelGeometry geo = new ChannelGeometry();
geo.setIrregularGeometryTF(false);
geo.setDepth(guessDepth);
geo.setBottomWidth(bottomWidth);
geo.setSideSlope(sideSlope);
//Calculate Froude number (in metric)
double A = geo.calcTrapezoidalArea();
double hydDepth = geo.calcTrapezoidalHydraulicDepth();
double fr = (discharge / A) / Math.sqrt(9.81 * hydDepth);
error = 1 - fr;
}
//Check if the result (Froude Number) based on the guessDepth is positive
//(like 'a') or negative (like 'b') and reset the limits of the
//bisection method to interpolate between the two new points
if(error > 0){
b = guessDepth;
}else{
a = guessDepth;
}
ctr++;
}
return guessDepth;
}
public void run() throws IOException{
//Determine if which depth to use for cross-section properties
if(depthType.equalsIgnoreCase("normal")){
//Use inputs to calculate a "normal" depth using Manning's equation
this.depth = calcNormalDepth();
}else if(depthType.equalsIgnoreCase("critical")){
//Use inputs to calculate a "critical" depth
this.depth = calcCriticalDepth();
}else if(depthType.equalsIgnoreCase("user")){
//Do nothing and just use the user-specified depth
}
//Calculate parameters of channel cross-section
if(irregularGeometry){
//Use irregular cross-section geometry calculations
calcIrregularArea();
calcIrregularPerimeter();
calcIrregularTopWidth();
calcIrregularHydraulicRadius();
calcIrregularHydraulicDepth();
}else{
//Use regular trapezoidal (rectangular and triangular included) cross-section geometry calculations
calcTrapezoidalArea();
calcTrapezoidalPerimeter();
calcTrapezoidalTopWidth();
calcTrapezoidalHydraulicRadius();
calcTrapezoidalHydraulicDepth();
}
//Round answers for significant digits
GeneralFunctions genFuncs = new GeneralFunctions();
this.depth = genFuncs.round(this.depth,3);
this.area = genFuncs.round(this.area,3);
this.perimeter = genFuncs.round(this.perimeter,3);
this.topWidth = genFuncs.round(this.topWidth,3);
this.hydraulicRadius = genFuncs.round(this.hydraulicRadius,3);
this.hydraulicDepth = genFuncs.round(this.hydraulicDepth,3);
}
public static void main(String[] args) throws IOException{
ChannelGeometry model = new ChannelGeometry();
model.run();
}
}