BriefFiniteElementNet / BriefFiniteElement.Net

BriefFiniteElementDotNET (BFE.NET) is a library for linear-static Finite Element Method (FEM) analysis of solids and structures in .NET
https://www.bfe-framework.net
GNU Lesser General Public License v3.0
154 stars 58 forks source link

An Error when get InterForce on a simple beam element model #183

Closed GIfengxin closed 4 months ago

GIfengxin commented 4 months ago

Sorry for the late reply. I have raise the issure #158 for the bug when returning InterForce at IsoPoint equals to -0.6666666(-2/3) last year. I have update the BriefFiniteElement.Net NuGet package to 2.1.2, but the issue remains. But this point doesn't contained in the IsoPoint List return by GetInternalForceDiscretationPoints method. Is it a wrong method been used and how can I avoid this Point, thank you.

epsi1on commented 4 months ago

Could you please put a little code that reproduces the error?

GIfengxin commented 4 months ago

Thank you for reply. I have 2 beam elements, the first one is (-300,0) to (0,0) and the second one is (0,0) to (300,0). For the first element, I load a concentrated force at (-250,0) in the Global -Z direction. And when I want to get the moment internal force at the same point, the error returned as I said before. If it's illegal to get the internal force at the concentrated force loaded point?

epsi1on commented 4 months ago

You are welcome, and thanks for information, if you put the simple code that throws the exception, i will try to help you. Could you please code what you described? I need to trace/debug the code in order to see what is going on and tell you if there is any problem, either with this library or your code

GIfengxin commented 4 months ago

Thank you for your reply. The example code is showed below, you can test it on simple project.

//Section
UniformParametric1DSection uniformParametric1DSection = new UniformParametric1DSection(6400, 1706670, 1706670, 1000);
//Material
UniformIsotropicMaterial uniformIsotropicMaterial = new UniformIsotropicMaterial(9350, 0.3);

BriefFiniteElementNet.Node node1 = new BriefFiniteElementNet.Node(-300, 0, 0){ Label = "n1" };
BriefFiniteElementNet.Node node2 = new BriefFiniteElementNet.Node(0, 0, 0){ Label = "n2" };
BriefFiniteElementNet.Node node3 = new BriefFiniteElementNet.Node(300, 0, 0){ Label = "n3" };

//Constraint
node1.Constraints = new Constraint(DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Released , DofConstraint.Fixed);
node2.Constraints = new Constraint(DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Released , DofConstraint.Fixed);
node3.Constraints = new Constraint(DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Released , DofConstraint.Fixed);

//Element
double eleLength = 300 - 0.0;
Element element1 = new BarElement(node1, node2)
{
    //Behavior = BarElementBehaviour.BeamZEulerBernoulli,
    Section = uniformParametric1DSection,
    Material = uniformIsotropicMaterial,
};
Element element2 = new BarElement(node2, node3)
{
    //Behavior = BarElementBehaviour.BeamZEulerBernoulli,
    Section = uniformParametric1DSection,
    Material = uniformIsotropicMaterial,
};

//element uniform load
element1.Loads.Add(new UniformLoad(LoadCase.DefaultLoadCase, new BriefFiniteElementNet.Vector(0, 0, 1), -0.0208, CoordinationSystem.Global));
element2.Loads.Add(new UniformLoad(LoadCase.DefaultLoadCase, new BriefFiniteElementNet.Vector(0, 0, 1), -0.0208, CoordinationSystem.Global));

//concentrated load
List<double> ptLoadCoords = new List<double>(){-250,-150,-50,50,150,250};

for (int i = 0; i < ptLoadCoords.Count(); i++)
{
    if (ptLoadCoords[i] < 0)  //left element
    {
        double isoPointCoord = (2.0 * (ptLoadCoords[i] - element1.Nodes[0].Location.X)- eleLength) / eleLength;
        IsoPoint isoPoint = new IsoPoint(isoPointCoord);
        element1.Loads.Add(new ConcentratedLoad(new Force(0, 0, -287.71875, 0, 0, 0) , isoPoint, CoordinationSystem.Global));
    }
    else //right element
    {
        double isoPointCoord = (2.0 * (ptLoadCoords[i] - element2.Nodes[0].Location.X) - eleLength) / eleLength;
        IsoPoint isoPoint = new IsoPoint(isoPointCoord);
        element2.Loads.Add(new ConcentratedLoad(new Force(0, 0, -287.71875, 0, 0, 0) , isoPoint, CoordinationSystem.Global));
    }
}

//Model
BriefFiniteElementNet.Model model = new BriefFiniteElementNet.Model();
model.Nodes.Add(node1, node2, node3);
model.Elements.Add(element1, element2);

//model.Solve();
model.Solve_MPC();

//Return Result

//Return IsoPoint, only (-1,0,0) and (1,0,0) returned
IsoPoint[] discreteIsoPoints = (element1 as BarElement).GetInternalForceDiscretationPoints();

//Use -250 to calculate a IsoPoint, which is equal to -2/3.
double currIsoPoint11 = (2.0 * (-250 - element1.Nodes[0].Location.X)- eleLength) / eleLength;
//Exception thrown: 'BriefFiniteElementNet.InvalidInternalForceLocationException' in BriefFiniteElementNet.dll
//The reason for the exception may be that a concentrated load is added on the same location
Force force = (element1 as BarElement).GetExactInternalForceAt(currIsoPoint11);

Looking forward to your reply.

epsi1on commented 4 months ago

Sorry for late reply, yes you are right. it is illegal to get internal force exactly at concentrated load point.

Imagine this case:

image

at center of beam, the moment is PL/4, but shear force do not have a specific value thus cannot be reported to user. Library will throw exception if user ask internal force at load point (center in this case).

for more information take a look at this example: LINK

GIfengxin commented 4 months ago

Thank you for your reply. I have use two points before and after the concentrated loaded points to get the internal force I want. Another question is that why the GetInternalForceDiscretationPoints() function could not return all the discretion points as I have comment ted early in this issue?

epsi1on commented 4 months ago

why the GetInternalForceDiscretationPoints() function could not return all the discretion points as I have comment ted early in this issue?

Not sure, but if you look at the source, will see that the code uses GetInternalForceDiscretationPoints() method to get the problematic points as a array, then check if given point is in the array and if exists, throws the exception. Could you please give more information about the case that BarElement.GetInternalForceDiscretationPoints() gives noncomplete list of points? preferably with code please...

Thanks

GIfengxin commented 4 months ago

Thank you for your reply. The code is the same as I sent before, which I pasted below, too. In the statement IsoPoint[] discreteIsoPoints = (element1 as BarElement).GetInternalForceDiscretationPoints();, the array discreteIsoPoints does not contain the IsoPoint -2/3 on which a concentrated force is loaded.

The code:

//Section
UniformParametric1DSection uniformParametric1DSection = new UniformParametric1DSection(6400, 1706670, 1706670, 1000);
//Material
UniformIsotropicMaterial uniformIsotropicMaterial = new UniformIsotropicMaterial(9350, 0.3);

BriefFiniteElementNet.Node node1 = new BriefFiniteElementNet.Node(-300, 0, 0){ Label = "n1" };
BriefFiniteElementNet.Node node2 = new BriefFiniteElementNet.Node(0, 0, 0){ Label = "n2" };
BriefFiniteElementNet.Node node3 = new BriefFiniteElementNet.Node(300, 0, 0){ Label = "n3" };

//Constraint
node1.Constraints = new Constraint(DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Released , DofConstraint.Fixed);
node2.Constraints = new Constraint(DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Released , DofConstraint.Fixed);
node3.Constraints = new Constraint(DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Fixed, DofConstraint.Released , DofConstraint.Fixed);

//Element
double eleLength = 300 - 0.0;
Element element1 = new BarElement(node1, node2)
{
    //Behavior = BarElementBehaviour.BeamZEulerBernoulli,
    Section = uniformParametric1DSection,
    Material = uniformIsotropicMaterial,
};
Element element2 = new BarElement(node2, node3)
{
    //Behavior = BarElementBehaviour.BeamZEulerBernoulli,
    Section = uniformParametric1DSection,
    Material = uniformIsotropicMaterial,
};

//element uniform load
element1.Loads.Add(new UniformLoad(LoadCase.DefaultLoadCase, new BriefFiniteElementNet.Vector(0, 0, 1), -0.0208, CoordinationSystem.Global));
element2.Loads.Add(new UniformLoad(LoadCase.DefaultLoadCase, new BriefFiniteElementNet.Vector(0, 0, 1), -0.0208, CoordinationSystem.Global));

//concentrated load
List<double> ptLoadCoords = new List<double>(){-250,-150,-50,50,150,250};

for (int i = 0; i < ptLoadCoords.Count(); i++)
{
    if (ptLoadCoords[i] < 0)  //left element
    {
        double isoPointCoord = (2.0 * (ptLoadCoords[i] - element1.Nodes[0].Location.X)- eleLength) / eleLength;
        IsoPoint isoPoint = new IsoPoint(isoPointCoord);
        element1.Loads.Add(new ConcentratedLoad(new Force(0, 0, -287.71875, 0, 0, 0) , isoPoint, CoordinationSystem.Global));
    }
    else //right element
    {
        double isoPointCoord = (2.0 * (ptLoadCoords[i] - element2.Nodes[0].Location.X) - eleLength) / eleLength;
        IsoPoint isoPoint = new IsoPoint(isoPointCoord);
        element2.Loads.Add(new ConcentratedLoad(new Force(0, 0, -287.71875, 0, 0, 0) , isoPoint, CoordinationSystem.Global));
    }
}

//Model
BriefFiniteElementNet.Model model = new BriefFiniteElementNet.Model();
model.Nodes.Add(node1, node2, node3);
model.Elements.Add(element1, element2);

//model.Solve();
model.Solve_MPC();

//Return Result

//Return IsoPoint, only (-1,0,0) and (1,0,0) returned
IsoPoint[] discreteIsoPoints = (element1 as BarElement).GetInternalForceDiscretationPoints();

//Use -250 to calculate a IsoPoint, which is equal to -2/3.
double currIsoPoint11 = (2.0 * (-250 - element1.Nodes[0].Location.X)- eleLength) / eleLength;
//Exception thrown: 'BriefFiniteElementNet.InvalidInternalForceLocationException' in BriefFiniteElementNet.dll
//The reason for the exception may be that a concentrated load is added on the same location
Force force = (element1 as BarElement).GetExactInternalForceAt(currIsoPoint11);`
epsi1on commented 4 months ago

Yes you are right. Thanks for reporting this issue. there were two versions of GetInternalForceDiscretationPoints:

I did removed the first one (usually anyone knows the node iso locations, which are -1 and +1, now it will throw exception on call) and only second one works now.

epsi1on commented 4 months ago

please note that code in github repo is updated (in 1f3744f), but nuget package not updated yet. Changes on ~github~ nuget package will be published in next nuget publish (soon). Thank you again for reporting the issue.

GIfengxin commented 4 months ago

Thank you for your work. For the current lastest nuget package(2.1.2), it I use the function GetInternalForceDiscretationPoints(LoadCase) function, can I get the correct result, that is all the discrete IsoPoints including concentrated load points will be returned.

epsi1on commented 4 months ago

good, then will close the issue for now. :)