BriefFiniteElementNet / BriefFiniteElement.Net

BriefFiniteElementDotNET (BFE.NET) is a library for linear-static Finite Element Method (FEM) analysis of solids and structures in .NET
GNU Lesser General Public License v3.0
148 stars 57 forks source link

problem with PartialNonUniformLoad on BarElement #111

Closed epsi1on closed 3 years ago

epsi1on commented 3 years ago

Describe the bug Adding a PartialNonUniformLoad to an BarElement do not cause any internal force or support reaction.

To Reproduce This is the code We've used and caused wrong result/ runtime error

public void Run3()
{
var model = new Model();

        //adding nodes
        model.Nodes.Add(new Node(0, 0, 0) { Label = "n0" });
        model.Nodes.Add(new Node(5, 5, 5) { Label = "n1" });

        //adding elements
        model.Elements.Add(new BarElement(model.Nodes["n0"], model.Nodes["n1"]) { Label = "e0" });

        //assign constraint to nodes
        model.Nodes["n0"].Constraints = model.Nodes["n1"].Constraints = Constraints.MovementFixed;

        //define sections and material
        var secAA = new Sections.UniformGeometric1DSection(SectionGenerator.GetISetion(0.24, 0.67, 0.01, 0.006));
        var mat = Materials.UniformIsotropicMaterial.CreateFromYoungPoisson(210e9, 0.3);

        //assign materials
        (model.Elements["e0"] as BarElement).Material = mat;

        //assign sections
        (model.Elements["e0"] as BarElement).Section = secAA;

        //creating loads
        var u1 = new Loads.PartialNonUniformLoad() { Direction = Vector.K, CoordinationSystem = CoordinationSystem.Global };
        u1.SeverityFunction = Mathh.SingleVariablePolynomial.FromPoints(-1, -6, 1, -4);
        u1.StartLocation = new IsoPoint(-0.5);      //set locations of trapezoidal load
        u1.EndLocation = new IsoPoint(0.5);         //set locations of trapezoidal load

        //assign loads
        model.Elements["e0"].Loads.Add(u1);

        //solve model
        model.Solve_MPC();

        //retrieve solve result
        var n0reaction = model.Nodes["n0"].GetSupportReaction();
        var n1reaction = model.Nodes["n1"].GetSupportReaction();

        Console.WriteLine("Support reaction of n0: {0}", n0reaction);
        Console.WriteLine("Support reaction of n1: {0}", n1reaction);

        var x = 1.0; //need to find internal force at x = 1.0 m
        var iso = (model.Elements["e0"] as BarElement).LocalCoordsToIsoCoords(x);//find the location of 1m in iso coordination system
        var e4Force = (model.Elements["e0"] as BarElement).GetInternalForceAt(iso[0]);//find internal force
        Console.WriteLine("internal force at x={0} is {1}", x, e4Force);
    }

Expected behavior There should be some non-zero support reaction.

Additional context From discussion #109

epsi1on commented 3 years ago

It gives nonzero values now, can you please do a check again?

cmoha commented 3 years ago

I will check this new commit. I did not receive a notification for your comment, sorry.

cmoha commented 3 years ago

I did some tests but results still don't match unfortunately

epsi1on commented 3 years ago

Can you please send the test code here and say what is expected value? I'll make it into a unit test... Thanks

cmoha commented 3 years ago

I have used the previous example, but changed frame element coordinates to be x-z plane for simplicity

         model.Nodes.Add(new Node(5, 0, 5) { Label = "n1" });

Case 1: load is along the bar element

            //creating loads
            var u1 = new Loads.PartialNonUniformLoad() { Direction = Vector.K, CoordinationSystem = CoordinationSystem.Local };
            double l = ((model.Elements["e0"] as BarElement).EndNode.Location - (model.Elements["e0"] as BarElement).StartNode.Location).Length;

            u1.SeverityFunction = Mathh.SingleVariablePolynomial.FromPoints(0, -6, l, -4);
            u1.StartLocation = new IsoPoint(-1);      //set locations of trapezoidal load
            u1.EndLocation = new IsoPoint(1);         //set locations of trapezoidal load
            //assign loads
            model.Elements["e0"].Loads.Add(u1);

            //solve model
            model.Solve_MPC();

            //retrieve solve result
            var n0reaction = model.Nodes["n0"].GetSupportReaction();
            var n1reaction = model.Nodes["n1"].GetSupportReaction();

            Debug.Assert(n0reaction == new Force(fx: -13.3333, fy: 0, fz: 13.3333, 0, 0, 0));
            Debug.Assert(n1reaction == new Force(fx: -11.6666, fy: 0, fz: 11.6666, 0, 0, 0));

Case 2: inverse the load (same code)

            u1.SeverityFunction = Mathh.SingleVariablePolynomial.FromPoints(0, -4, l, -6);

             Debug.Assert(n0reaction == new Force(fx: -11.6666, fy: 0, fz: 11.6666, 0, 0, 0));
             Debug.Assert(n1reaction == new Force(fx: -13.3333, fy: 0, fz: 13.3333, 0, 0, 0));
epsi1on commented 3 years ago

Are you sure about the expected support reactions? first node is on 0,0,0 and second node is on 5,5,5 so length of beam is about 8.6m. When i replace the distributed load with approximate series of concentrated loads (using Util.AreaLoad2ConcentratedLoads()) , i'll get same results, so result seems right to me.

image

https://github.com/BriefFiniteElementNet/BriefFiniteElement.Net/blob/master/BriefFiniteElementNet.Validation/GithubIssues/Issue111.cs

Thanks

cmoha commented 3 years ago

Kindly note that I Changed second node to 5,0,5 and changed load definition

I have used the previous example, but changed frame element coordinates to be x-z plane for simplicity

         model.Nodes.Add(new Node(5, 0, 5) { Label = "n1" });

image

. . also running same example in issue 111 (old one) gives not expected results image

epsi1on commented 3 years ago

also note that xi from element's coordinate is fed to SingleVariablePolynomial directly. So for your example it should be something like this:

u1.SeverityFunction = Mathh.SingleVariablePolynomial.FromPoints(-1, -6, 1, -4);
u1.StartLocation = new IsoPoint(-1);      //set locations of trapezoidal load
u1.EndLocation = new IsoPoint(1);         //set locations of trapezoidal load

but you used this:

u1.SeverityFunction = Mathh.SingleVariablePolynomial.FromPoints(0, -6, l, -4);
u1.StartLocation = new IsoPoint(-1);      //set locations of trapezoidal load
u1.EndLocation = new IsoPoint(1);         //set locations of trapezoidal load

i.e. parameters of Mathh.SingleVariablePolynomial.FromPoints(xi1,magnitude1,xi2,magnitude2) are in element's iso parametric coordination system...

epsi1on commented 3 years ago

I've checked again, sorry for dirty code this is result for first model and results are same as expected image

But didn't know the length of beam of second model to check that out. Can you please give dimensions of second model? beam length and start and end locations that distributed load is applied . I mean this:

image

Thanks

cmoha commented 3 years ago

Good news, the results are getting as expected element's coordinates are model.Nodes.Add(new Node(0, 0, 0) { Label = "n0" }); model.Nodes.Add(new Node(5, 5, 5) { Label = "n1" }); load is var u1 = new Loads.PartialNonUniformLoad() { Direction = Vector.K, CoordinationSystem = CoordinationSystem.Global }; u1.SeverityFunction = Mathh.SingleVariablePolynomial.FromPoints(-0.5, -6, 0.5, -4); u1.StartLocation = new IsoPoint(-0.5); u1.EndLocation = new IsoPoint(0.5);

i'll try to make more test cases to check expected results

epsi1on commented 3 years ago

OK, reopen this issue if problem still exists