using System; using System.Collections.Generic; // ======================================================================================================================= // Class name: TSpencerPlane // // Description: // // Copyright (c) 1999-2008 Deltares // // Date ID Modification // 2008-08-12 Best Created // ======================================================================================================================= namespace Deltares.Stability.Calculation.Inner { public class TSpencerPlane : TSlipPlane { // in case of uplift calculation // internal data private const double CDeltaShift = 0.01; // private double FRestForce = 0; // private int FNIterations = 0; private double FDelta = 0; private string FErrorNote = String.Empty; // vaste hoek interlamel krachten private bool FReverseGeometry = false; private double FSafety = 0; // private double FRestMoment = 0; private TPoints[] FSlipPlanePoints; // in case of uplift calculation private double FXActiveCentre = 0; private double FXPassiveCentre = 0; private double FZActiveCentre = 0; private double FZPassiveCentre = 0; private double FZTangent = 0; // ======================================================================================================================= // Date ID Modification // 2008-08-12 Best Created // ======================================================================================================================= //Constructor Create() public TSpencerPlane() : base() { // TODO: Add any constructor code here } public double XActiveCentre { get { return FXActiveCentre; } set { FXActiveCentre = value; } } public double ZActiveCentre { get { return FZActiveCentre; } set { FZActiveCentre = value; } } public double XPassiveCentre { get { return FXPassiveCentre; } set { FXPassiveCentre = value; } } public double ZPassiveCentre { get { return FZPassiveCentre; } set { FZPassiveCentre = value; } } public double ZTangent { get { return FZTangent; } set { FZTangent = value; } } // in case of a default spencer calculation public TPoints[] SlipPlanePoints { get { return FSlipPlanePoints; } set { FSlipPlanePoints = value; } } public double SafetyFactor { get { return GetStabilityFactor(); } } public double Safety { get { return FSafety; } set { FSafety = value; } } public void DetermineResultantMoment() //============================================== { for (int sliceIndex = 0; sliceIndex < FNSlices; sliceIndex++) { FSlices[sliceIndex].DetermineResultantMoment(); } } // ======================================================================================================================= // Date ID Modification // 2008-08-13 Best Created // ======================================================================================================================= public override void CreateSlices() { TSliceDivision LSliceDivision; int i; // Find intersection points with surface LSliceDivision = new TSliceDivision(); LSliceDivision.CalculationType = FInterfaceData.ModelData.CalcType; LSliceDivision.LeftSideIsActive = FInterfaceData.ModelData.LeftSideIsActive; LSliceDivision.StriveSlice = Math.Max(1, FInterfaceData.CalculationOptions.StriveSlice); LSliceDivision.XLeftSurface = FxLeftSurface; LSliceDivision.XRightSurface = FXRightSurface; LSliceDivision.SurfaceLine = FInterfaceData.SurfaceLine; LSliceDivision.GeometryPoints = FInterfaceData.GeometryPoints; LSliceDivision.Layers = FInterfaceData.Layers; LSliceDivision.WaterLines = FInterfaceData.CurrentWaternet.Waterlines; LSliceDivision.PLLines = FInterfaceData.CurrentWaternet.PLLines; if (FInterfaceData.CurrentWaternet.PhreaticLineNumber >= 0) { LSliceDivision.FreaticLine = FInterfaceData.CurrentWaternet.PLLines[FInterfaceData.CurrentWaternet.PhreaticLineNumber]; } if ((FInterfaceData.ModelData.CalcType == TCalculationTypeSet.csLiftSpencer)) { LSliceDivision.XMidActive = FXActiveCentre; LSliceDivision.ZMidActive = FZActiveCentre; LSliceDivision.XMidpassive = FXPassiveCentre; LSliceDivision.ZMidpassive = FZPassiveCentre; LSliceDivision.ZTangent = FZTangent; } else { LSliceDivision.SlipPoints = (TPoints[]) FSlipPlanePoints.Clone(); } LSliceDivision.DivideFailureSurfaceInSlices(); FNSlices = LSliceDivision.NSlices; FSlices = new Tslice[FNSlices]; for (i = 0; i < FNSlices; i++) { FSlices[i] = new Tslice(LSliceDivision.Slices[i]); } } // ======================================================================================================================= // Date ID Modification // 2008-08-12 Best Created // ======================================================================================================================= protected override void RemoveLocalPointers() { FNSlices = 0; FSlipPlanePoints = null; } protected void FindUpliftSurfaceIntersectionPoints(double AXMid, double AZMid, double ARad, List AXIntersect, List AZIntersect) { // find intersection points for (int i = FInterfaceData.SurfaceLine.GetLowerBound(0) + 1; i <= FInterfaceData.SurfaceLine.GetUpperBound(0); i++) { double LXStart = FInterfaceData.SurfaceLine[i - 1].XCoor; double LZStart = FInterfaceData.SurfaceLine[i - 1].ZCoor; double LXEnd = FInterfaceData.SurfaceLine[i].XCoor; double LZEnd = FInterfaceData.SurfaceLine[i - 1].ZCoor; MStabDatafunctions.SnijCirkelLijnstuk(AXMid, AZMid, ARad, LXStart, LZStart, LXEnd, LZEnd, AXIntersect, AZIntersect); } MStabDatafunctions.SortIntersectionPoints(AZIntersect, AXIntersect); MStabDatafunctions.RemoveDoubleIntersectionPoints(AXIntersect, AZIntersect); } // ======================================================================================================================= // Date ID Modification // 2008-06-23 Best Created // ======================================================================================================================= protected void FindUpliftEntreepoints() { var LXIntersect = new List(); var LZIntersect = new List(); double LXMid; double LZMid; double LRad; if (FXActiveCentre > FXPassiveCentre) { LXMid = FXPassiveCentre; LZMid = FZPassiveCentre; LRad = FZPassiveCentre - FZTangent; } else { LXMid = FXActiveCentre; LZMid = FZActiveCentre; LRad = FZActiveCentre - FZTangent; } FindUpliftSurfaceIntersectionPoints(LXMid, LZMid, LRad, LXIntersect, LZIntersect); bool LPointTooHigh = false; if ((LXIntersect.Count >= 2)) { int i = 0; bool LFound = false; // find entree point (left side) closest to center point while ((i <= LXIntersect.Count - 1) && !LFound) { LFound = (LXIntersect[i] <= LXMid); if (LFound) { FxLeftSurface = LXIntersect[i]; FZLeftSurface = LZIntersect[i]; LPointTooHigh = (FZLeftSurface > LZMid); } i++; } } else { FErrorNote = FErrorNote + Constants.cNoExitPointFound; } if (LPointTooHigh) { FErrorNote = FErrorNote + Constants.cExitPointTooHigh; } } protected void FindSpencerEntreePoints() { int LNSlipPoints = FSlipPlanePoints.Length; if (SlipPlanePoints[0].XCoor > SlipPlanePoints[LNSlipPoints - 1].XCoor) { FxLeftSurface = SlipPlanePoints[LNSlipPoints - 1].XCoor; FZLeftSurface = SlipPlanePoints[LNSlipPoints - 1].ZCoor; FXRightSurface = SlipPlanePoints[0].XCoor; FZRightSurface = SlipPlanePoints[0].ZCoor; } else { FxLeftSurface = SlipPlanePoints[0].XCoor; FZLeftSurface = SlipPlanePoints[0].ZCoor; FXRightSurface = SlipPlanePoints[LNSlipPoints - 1].XCoor; FZRightSurface = SlipPlanePoints[LNSlipPoints - 1].ZCoor; } } protected void CheckSpencerForbiddenLine() { TForbiddenLine LForbiddenLine; bool LIntersects = false; int j; int LNDim = FInterfaceData.ForbiddenLines.Length; if ((LNDim > 0)) { int i = 0; while ((i < LNDim) && (!LIntersects)) { LForbiddenLine = FInterfaceData.ForbiddenLines[i]; int LNSlipPoints = SlipPlanePoints.Length; for (j = 0; j < LNSlipPoints - 1; j++) { if (LForbiddenLine.LineIntersectsForbiddenLine(SlipPlanePoints[j].XCoor, SlipPlanePoints[j].ZCoor, SlipPlanePoints[j + 1].XCoor, SlipPlanePoints[j + 1].ZCoor)) { LIntersects = true; } } } if (LIntersects) { FErrorNote = FErrorNote + Constants.cIntersectsForbiddenLine; } } } // ======================================================================================================================= // Date ID Modification // 2008-06-23 Best Created // ======================================================================================================================= protected void FindUpliftExitpoints() { var LXIntersect = new List(); var LZIntersect = new List(); bool LPointTooHigh = false; double LXMid; double LZMid; double LRad; if (FXActiveCentre > FXPassiveCentre) { LXMid = FXActiveCentre; LZMid = FZActiveCentre; LRad = FZActiveCentre - FZTangent; } else { LXMid = FXPassiveCentre; LZMid = FZPassiveCentre; LRad = FZPassiveCentre - FZTangent; } FindUpliftSurfaceIntersectionPoints(LXMid, LZMid, LRad, LXIntersect, LZIntersect); if ((LXIntersect.Count >= 2)) { int i = 1; bool LFound = false; // find exit point (Right side) closest to center point while ((i <= LXIntersect.Count - 1) && !LFound) { LFound = (LXIntersect[i] >= LXMid); if (LFound) { FXRightSurface = LXIntersect[i]; FZRightSurface = LZIntersect[i]; LPointTooHigh = (LZIntersect[i] > LZMid); } i++; } } else { FErrorNote = FErrorNote + Constants.cNoExitPointFound; } if (LPointTooHigh) { FErrorNote = FErrorNote + Constants.cExitPointTooHigh; } } // ======================================================================================================================= // Date ID Modification // 2008-06-25 Best Created // ======================================================================================================================= protected void CheckOnForbiddenLines() { int i; double LRadius; int LNDim; bool LIntersects; LIntersects = false; double LXLeft; double LZLeft; double LXRight; double LZRight; if (FXPassiveCentre > FXActiveCentre) { LXLeft = FXActiveCentre; LZLeft = FZActiveCentre; LXRight = FXPassiveCentre; LZRight = FZPassiveCentre; } else { LXRight = FXActiveCentre; LZRight = FZActiveCentre; LXLeft = FXPassiveCentre; LZLeft = FZPassiveCentre; } LNDim = FInterfaceData.ForbiddenLines.Length; if ((LNDim > 0)) { i = 0; while ((i < LNDim) && !LIntersects) { TForbiddenLine _wvar1 = FInterfaceData.ForbiddenLines[i]; // halve cirkel aan linkerkant LRadius = LZLeft - FZTangent; if (!LIntersects) { LIntersects = _wvar1.CircleIntersectsForbiddenLine(LXLeft, LZLeft, LRadius, FxLeftSurface, LXLeft); } if (!LIntersects) { LIntersects = _wvar1.LineIntersectsForbiddenLine(LXLeft, FZTangent, LXLeft, FZTangent); } LRadius = LZRight - FZTangent; if (!LIntersects) { LIntersects = _wvar1.CircleIntersectsForbiddenLine(LXRight, LZRight, LRadius, LXRight, FXRightSurface); } i++; } } if (LIntersects) { FErrorNote = FErrorNote + Constants.cIntersectsForbiddenLine; } } // ======================================================================================================================= // Date ID Modification // 2010-04-26 Created // ======================================================================================================================= protected bool IsUpliftSlipPlaneValid() { bool result; string lStr; try { FindUpliftEntreepoints(); FindUpliftExitpoints(); // See if intersection with forbidden lines are found if ((FErrorNote == "")) { CheckOnForbiddenLines(); } result = (FErrorNote == ""); } catch (Exception E) { lStr = E.Message; throw E; } return result; } protected bool IsSpencerSlipPlaneValid() //============================================== { FindSpencerEntreePoints(); CheckSpencerForbiddenLine(); return (FErrorNote == ""); } // ======================================================================================================================= // Date ID Modification // 1999-06-09 Best Created // ======================================================================================================================= private void DetermineShearForces() { int i; double LTotpore; double LTanPhiBot; double LNormalstress; double LCohBottom; double LLeftCombination; double LRightCombination; double LTotalPore; double LCosDif; double LCosBot; double LMass; for (i = 0; i < FNSlices; i++) { Tslice lSlice = FSlices[i]; if ((lSlice.ArcLength != 0)) { LLeftCombination = lSlice.LeftForce*Math.Sin(-lSlice.BottomAngle - lSlice.LeftForceAngle); LRightCombination = lSlice.RightForce*Math.Sin(-lSlice.BottomAngle - lSlice.RightForceAngle); LTotpore = Math.Sqrt(lSlice.VPoreOnSurface*lSlice.VPoreOnSurface + lSlice.HPoreOnSurface*lSlice.HPoreOnSurface); LTotalPore = lSlice.TotalPore; LCosDif = Math.Cos(lSlice.TopAngle - lSlice.BottomAngle); LCosBot = Math.Cos(-lSlice.BottomAngle); LMass = lSlice.Weight + lSlice.ExternalLoad + lSlice.SigmaSoilQuake*lSlice.Width; lSlice.NormalStress = (LTotpore*LCosDif - LTotalPore*lSlice.ArcLength + LMass*LCosBot - LLeftCombination + LRightCombination)/lSlice.ArcLength; } else { lSlice.NormalStress = 0; } LTanPhiBot = Math.Tan(lSlice.PhiBottom*Constants.CRad); LNormalstress = lSlice.NormalStress; LCohBottom = lSlice.CohBottom; lSlice.ShearStress = (LNormalstress*LTanPhiBot + LCohBottom)/FSafety; } } // ======================================================================================================================= // Date ID Modification // 2008-08-12 Best Created // ======================================================================================================================= private void CorrectForceAngles() { int i; for (i = 0; i < FNSlices; i++) { Tslice _wvar1 = FSlices[i]; _wvar1.LeftForceAngle = FDelta; _wvar1.RightForceAngle = FDelta; } // de start en end hoeken zijn nul FSlices[0].LeftForceAngle = 0.0; FSlices[FNSlices - 1].RightForceAngle = 0; } // ======================================================================================================================= // Date ID Modification // 2008-08-12 Best Created // ======================================================================================================================= private void CorrectSafetyFactor() { int i; for (i = 0; i < FNSlices; i++) { Tslice _wvar1 = FSlices[i]; _wvar1.TmpSafety = FSafety; } } // ======================================================================================================================= // Date ID Modification // 2008-08-12 Best Created // ======================================================================================================================= private void DeltaGrens(ref double AMinDelta, ref double AMaxDelta) { int i; double Ltmpphi; double LDeltaMin; double LDeltaMax; AMinDelta = -0.5*Math.PI; AMaxDelta = 0.5*Math.PI; for (i = 0; i < FNSlices; i++) { Tslice _wvar1 = FSlices[i]; if (_wvar1.PhiBottom*Constants.CRad < 1.0E-12) { Ltmpphi = 1.0E-12; } else { Ltmpphi = _wvar1.PhiBottom*Constants.CRad; } // minimale grenswaarde bepalen LDeltaMin = -FSlices[i].BottomAngle - (Math.PI + Math.Atan(-FSafety/Math.Tan(Ltmpphi))); AMinDelta = Math.Max(LDeltaMin, AMinDelta); // maximale grenswaarde bepalen LDeltaMax = LDeltaMin + Math.PI; AMaxDelta = Math.Min(LDeltaMax, AMaxDelta); } } // ======================================================================================================================= // Date ID Modification // 2008-08-12 Best Created // ======================================================================================================================= private double DetermineRestForce() { double result; int i; // de eerste lamel wordt een kracht nul aangenomen Tslice _wvar1 = FSlices[0]; _wvar1.LeftForce = 0; _wvar1.RightForce = _wvar1.GetRightForce(); // steeds wordt uitgaande van linker kracht // de rechter berekend for (i = 1; i < FNSlices; i++) { Tslice _wvar2 = FSlices[i]; FSlices[i].LeftForce = FSlices[i - 1].RightForce; FSlices[i].RightForce = _wvar2.GetRightForce(); } result = FSlices[FNSlices - 1].RightForce; return result; } // ======================================================================================================================= // Date ID Modification // 2008-08-12 Best Created // ======================================================================================================================= private double GetCriticalSafety() { double result; int i; double FMax; double Fk; FMax = 0.1; // minimale veiligheid per lamel bepalen for (i = 0; i < FNSlices; i++) { Tslice _wvar1 = FSlices[i]; Fk = -Math.Tan(-_wvar1.BottomAngle - _wvar1.LeftForceAngle)* Math.Tan(_wvar1.PhiBottom*Constants.CRad); FMax = Math.Max(Fk, FMax); } result = FMax; return result; } // ======================================================================================================================= // Date ID Modification // 2008-08-12 Best Created // ======================================================================================================================= private void DeleteRestMoments(ref double Mindelta, ref double MaxDelta, ref double FStart) { double LRestMoment; bool LReady; int LLocint; DeltaGrens(ref Mindelta, ref MaxDelta); LLocint = 0; do { LLocint++; LRestMoment = DetermineRestMoment(); LReady = (Math.Abs(LRestMoment) < Constants.CRestMomCriterium); if (!LReady) { ModifyDelta(LRestMoment); if ((FDelta >= MaxDelta) || (FDelta <= Mindelta)) { FStart = FStart + 0.1; FSafety = FStart; CorrectSafetyFactor(); DeltaGrens(ref Mindelta, ref MaxDelta); FDelta = (Mindelta + MaxDelta)/2.0; CorrectForceAngles(); } } } while (!(LReady || (LLocint > Constants.CMaxIter))); // writeln(lft,sRestmomentNa[GLO],locint:4,sIteratiesIs[GLO],RestMoment:12:5); } // ======================================================================================================================= // Date ID Modification // 2008-08-12 Best Created // ======================================================================================================================= private void ModifySafety(double ARestforce) { double LFold; double LDf; double LFRest1; double LFRest2; double LDFrest; double LFMin; double LFMax; LDf = 0.001; LFold = FSafety; LFRest1 = ARestforce; FSafety = LFold + LDf; CorrectSafetyFactor(); LFMax = 999.0; // Minimale veiligheids coefficient ivm het oplossings gebied bepalen LFMin = GetCriticalSafety(); // Bepaling restkracht bij factor nieuwe factor LFRest2 = DetermineRestForce(); LDFrest = LFRest2 - LFRest1; if ((Math.Abs(LDFrest) < 1E-12)) { if ((LDFrest >= 0)) { LDFrest = 1.0E-12; } else { LDFrest = -1.0E-12; } } // Nieuwe veiligheids coeficient bepalen met behulp van de gradient methode FSafety = LFold - (LDf/LDFrest)*LFRest1; // writeln(lft,FSafety:12:5,' ',Fold:12:5,' ',DF :12:5,' ',DFRest:12:5,' ',FRest1:12:5,' ',frest2:12:5); // Testen of veiligheidsCoefficient binen gestelde grenzen ligt if ((FSafety < LFMin)) { FDelta = FDelta + 0.01; CorrectForceAngles(); FSafety = LFold; CorrectSafetyFactor(); LFMin = GetCriticalSafety(); // Test of veiligheidscoefficient nog in oplosingsgebied valt if ((FSafety <= LFMin)) { FSafety = LFMin + 0.5; } } if ((FSafety > LFMax)) { FDelta = FDelta - 0.01; CorrectForceAngles(); FSafety = LFMax; } CorrectSafetyFactor(); } // ======================================================================================================================= // Date ID Modification // 2008-08-12 Best Created // ======================================================================================================================= private void DeleteRestforces() { double restforce; bool ready; int Locint; Locint = 0; do { Locint++; restforce = DetermineRestForce(); ready = (Math.Abs(restforce) < Constants.CRestForceCriterium); if (!ready) { ModifySafety(restforce); } } while (!(ready || (Locint > Constants.CMaxIter))); // writeln(lft,sRestforceNa[GLO],locint:4,sIteratiesIs[GLO],Restforce:12:5); } // ======================================================================================================================= // Date ID Modification // 2008-08-12 Best Created // ======================================================================================================================= private double DetermineRestMoment() { double result; int i; // initialiseer de eerst arm Tslice _wvar1 = FSlices[0]; FSlices[0].LeftForceZ = _wvar1.ZBottomMid; _wvar1.RightforceZ = _wvar1.GetRightForceZ(); if ((_wvar1.RightforceZ < _wvar1.ZBottomRight)) { _wvar1.RightforceZ = _wvar1.ZBottomRight; } if ((_wvar1.RightforceZ > _wvar1.ZTopRight)) { _wvar1.RightforceZ = _wvar1.ZTopRight; } // bereken vervolgens steeds de rechter arm uit for (i = 1; i <= (FNSlices - 2); i++) { Tslice _wvar2 = FSlices[i]; _wvar2.LeftForceZ = FSlices[i - 1].RightforceZ; _wvar2.RightforceZ = _wvar2.GetRightForceZ(); } // Bepaling rest moment laatste lamel Tslice _wvar3 = FSlices[FNSlices - 1]; _wvar3.LeftForceZ = FSlices[FNSlices - 2].RightforceZ; result = _wvar3.LeftForce*Math.Cos(_wvar3.LeftForceAngle)* (_wvar3.LeftForceZ - _wvar3.ZBottomMid + 0.5*_wvar3.Width*Math.Tan(-_wvar3.BottomAngle)) - _wvar3.LeftForce*Math.Sin(_wvar3.LeftForceAngle)*0.5*_wvar3.Width - _wvar3.HPoreOnSurface*(_wvar3.ZTopMid - _wvar3.ZBottomMid); return result; } // ======================================================================================================================= // Date ID Modification // 2008-08-12 Best Created // ======================================================================================================================= private void ModifyDelta(double restMoment) { double DDelta; double Delta1; double Restm1; double Restm2; double DRestMoment; DDelta = 0.01; Delta1 = FDelta; Restm1 = restMoment; // bepaal het restmoment met een andere delta FDelta = Delta1 + DDelta; CorrectForceAngles(); Restm2 = DetermineRestMoment(); // nieuwe delta bepalen met de gradient methode DRestMoment = Restm2 - Restm1; try { FDelta = Delta1 - (DDelta/DRestMoment)*Restm1; } catch { // indien drestmom = 0 dan ook geen variatie in delta } CorrectForceAngles(); } // ======================================================================================================================= // Date ID Modification // 2008-08-13 Best Created // ======================================================================================================================= private void ReverseData() { int i; int j; int Last; double Temp; Temp = FxLeftSurface; FxLeftSurface = -FXRightSurface; FXRightSurface = -Temp; // zet de coordinaten op een andere plaats Tslice lSlice; for (i = 0; i < FNSlices; i++) { lSlice = FSlices[i]; // xcoors Temp = lSlice.xLeft; lSlice.xLeft = -lSlice.xRight; lSlice.xRight = -Temp; // ZCoors Temp = lSlice.ZBottomLeft; lSlice.ZBottomLeft = lSlice.ZBottomRight; lSlice.ZBottomRight = Temp; Temp = lSlice.ZTopLeft; lSlice.ZTopLeft = lSlice.ZTopRight; lSlice.ZTopRight = Temp; Temp = lSlice.LeftForceZ; lSlice.LeftForceZ = lSlice.RightforceZ; lSlice.RightforceZ = Temp; Temp = lSlice.LeftForce; lSlice.LeftForce = lSlice.RightForce; lSlice.RightForce = Temp; // Loading lSlice.HPoreOnSurface = -lSlice.HPoreOnSurface; lSlice.LeftForceAngle = -lSlice.LeftForceAngle; lSlice.RightForceAngle = -lSlice.RightForceAngle; } // wissel de slices om (aleen de adressen) Last = Convert.ToInt32(FNSlices/2); for (i = 0; i < Last; i++) { j = FNSlices - 1 - i; swapslices(i, j); } for (i = 0; i < FNSlices; i++) { FSlices[i].addCalculatedProperties(); } } private void PrepareCalculation() { FGamWat = FInterfaceData.CurrentWaternet.GammaWater; } // ======================================================================================================================= // Date ID Modification // 2008-06-25 Best Created // ======================================================================================================================= private void ForceSolutionWithinRestraints() //============================================== { int LSliceIndex; double LInverseSafety; double LDelta; double LMaximumDelta; double LMinimumDelta; //Confine solution space; the constraints are found in equation (9) LInverseSafety = 1/FSafety + Constants.CInverseSafetyShift; if (LInverseSafety < Constants.CInverseSafetyShift) { // throw EDSpencerPlaneError.Create('New safety factor is negative; no proper solution is anymore expected'); } // The resultant slope range is sought for LInversSafety; then, also valid for LInversSafety - CInverseSafetyShift LMinimumDelta = double.MinValue; LMaximumDelta = double.MaxValue; // Find the Delta range belonging to LSafety for (LSliceIndex = 0; LSliceIndex < FNSlices; LSliceIndex++) { // BottomAngle is defined negative here LMinimumDelta = Math.Max(LMinimumDelta, (-Math.PI/2 - FSlices[LSliceIndex].BottomAngle)); LMaximumDelta = Math.Min(LMaximumDelta, (Math.PI/2 - Math.Atan( Math.Tan(FSlices[LSliceIndex].PhiBottom*Constants.CDegreeToRad)* LInverseSafety) - FSlices[LSliceIndex].BottomAngle)); } // Stay away CDeltaShift from the restraints; reserve space for LDelta - CDeltaShift if (LMinimumDelta > LMaximumDelta - 3*CDeltaShift) { // throw EDSpencerPlaneError.Create('No valid range for resultant slope; no proper solution is anymore expected'); } LDelta = Math.Max(LMinimumDelta + 2*CDeltaShift, Math.Min(FDelta, LMaximumDelta - CDeltaShift)); //{ Allow adapted proposed shift of resultant slope FDelta = LDelta; } private void DetermineInterSliceForces() //============================================== { int LSliceIndex; Tslice LSlice; LSlice = FSlices[0]; // The first slice has no left force } LSlice.LeftForce = 0; LSlice.LeftForceZ = LSlice.ZBottomLeft; LSlice.LeftForceAngle = FDelta; // Determine the right force and slope and pass them to the left values of the next slice } for (LSliceIndex = 0; LSliceIndex < FNSlices; LSliceIndex++) { LSlice = FSlices[LSliceIndex]; // Determine resultant } LSlice.TmpSafety = FSafety; LSlice.ResultantAngle = FDelta; LSlice.DetermineResultantNorm(); LSlice.DetermineResultantForce(); } for (LSliceIndex = 0; LSliceIndex < FNSlices - 1; LSliceIndex++) { LSlice = FSlices[LSliceIndex]; LSlice.DetermineInterSliceForces(); // Pass to next slice } FSlices[LSliceIndex + 1].LeftForce = LSlice.RightForce; FSlices[LSliceIndex + 1].LeftForceZ = LSlice.RightforceZ; FSlices[LSliceIndex + 1].LeftForceAngle = LSlice.RightForceAngle; } // It is neat to specify proper end values of the inter-slice slopes, though they are not relevant } LSlice = FSlices[FNSlices - 1]; LSlice.RightForce = 0; LSlice.RightforceZ = LSlice.ZBottomRight; LSlice.RightForceAngle = FDelta; } private double DetermineMomentImbalance() //============================================== { double result = 0; Tslice LSlice; // See equation (5); the right inter-slice force should vanish at the last slice } LSlice = FSlices[FNSlices - 1]; result = LSlice.ResultantMoment + LSlice.LeftForce*Math.Cos(LSlice.LeftForceAngle)* (LSlice.LeftForceZ - LSlice.ZBottomMid - Math.Tan(LSlice.LeftForceAngle)*LSlice.Width/2.0); return result; } private double DetermineForceImbalance() //============================================== { double result = 0; Tslice LSlice; //See equation (10) and (13); the right inter-slice force should vanish at the last slice. In (10) the ForceAngle's // are equal. We use this routine for both methods, Constant Resultant Slope and Equally Spaced Action ratio } LSlice = FSlices[FNSlices - 1]; result = LSlice.LeftForce*Math.Cos(LSlice.LeftForceAngle - LSlice.RightForceAngle) + LSlice.ResultantForce*Math.Cos(LSlice.ResultantAngle - LSlice.RightForceAngle); return result; } private bool IterateConstantResultantSlopes() //============================================== { const double CSafetyAccuracy = 0.0001; const double CDeltaAccuracy = 0.0001; bool result; double LRestForce; double LRestMoment; double LMomentPerSafety; double LForcePerSafety; double LMomentPerDelta; double LForcePerDelta; double LInverseSafetyShift; double LDeltaShift; // Determine present imbalance DetermineInterSliceForces(); LRestMoment = DetermineMomentImbalance(); LRestForce = DetermineForceImbalance(); // Determine imbalance for variation in resultant slope } FDelta = FDelta - CDeltaShift; DetermineInterSliceForces(); LMomentPerDelta = (LRestMoment - DetermineMomentImbalance())/CDeltaShift; LForcePerDelta = (LRestForce - DetermineForceImbalance())/CDeltaShift; FDelta = FDelta + CDeltaShift; // etermine imbalance for variation in safety factor } FSafety = FSafety/(1 + FSafety*Constants.CInverseSafetyShift); DetermineInterSliceForces(); LMomentPerSafety = (DetermineMomentImbalance() - LRestMoment)/Constants.CInverseSafetyShift; LForcePerSafety = (DetermineForceImbalance() - LRestForce)/Constants.CInverseSafetyShift; FSafety = FSafety/(1 - FSafety*Constants.CInverseSafetyShift); // ermine the adjustments } LDeltaShift = LForcePerSafety*LMomentPerDelta - LForcePerDelta*LMomentPerSafety; LInverseSafetyShift = (LMomentPerDelta*LRestForce - LForcePerDelta*LRestMoment)/LDeltaShift; LDeltaShift = (LForcePerSafety*LRestMoment - LMomentPerSafety*LRestForce)/LDeltaShift; //Check accuracy and adjust } result = (Math.Abs(LInverseSafetyShift) < CSafetyAccuracy) && (Math.Abs(LDeltaShift) < CDeltaAccuracy); // Result := ((Abs(LRestForce) < CEpsForce) and (Abs(LRestMoment) < CEpsMom)); FSafety = FSafety/(1 - FSafety*LInverseSafetyShift); FDelta = FDelta - LDeltaShift; return result; } // ======================================================================================================================= // Date ID Modification // 2008-08-12 Best Created // ======================================================================================================================= private double GetStabilityFactor() { int LIter = 0; bool lContinue; // double LFStart = 0; bool LIsAccurate = false; FGamWat = FInterfaceData.CurrentWaternet.GammaWater; FSafety = Constants.CFMinStart; if (FInterfaceData.ModelData.CalcType == TCalculationTypeSet.csLiftSpencer) { lContinue = IsUpliftSlipPlaneValid(); } else { lContinue = IsSpencerSlipPlaneValid(); } if ((lContinue)) { // assignfile(lft,sTestout[GLO]); // rewrite(lft); CreateSlices(); FillSlices(); // DumpSliceData; lContinue = (FNSlices > 0); } if ((lContinue)) { FReverseGeometry = FSlices[0].ZTopMid < FSlices[FNSlices - 1].ZTopMid; // Fill soil,stress data of slices ( in een apparte routine STSliceFill if (FReverseGeometry) { // om uit te kunnene rekenen moet de geometrie gespiegeld worden // dit gebeurd door alle xen van een - teken te voorzien en de lamelen // andersom te sorteren ReverseData(); } LIter = 0; FDelta = 0; FSafety = Constants.CFMinStart; DetermineResultantMoment(); // { For safety's sake check restraints } ForceSolutionWithinRestraints(); while (!LIsAccurate && LIter <= Constants.CMaxIterations) { LIsAccurate = IterateConstantResultantSlopes(); ForceSolutionWithinRestraints(); LIter++; if ((LIter == Constants.CMaxIterations) || (FSafety < 0)) { /*throw new Exception( "Gradient approach not successful; no proper solution is anymore expected");*/ FSafety = Constants.CFMinStart; break; } } DetermineInterSliceForces(); DetermineShearForces(); if (FReverseGeometry) { ReverseData(); } } return FSafety; } } }