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

An calculation error about a simple beam element model #158

Closed GIfengxin closed 12 months ago

GIfengxin commented 1 year ago

Hello, i have build a simple beam element model with this brilliant library, but when I want to solve the model, I got an error said "Matrix must be symmetric positive definite", I can not figure out what's wrong with the model. The model has 3 nodes with label named n1, n2 and n3, the nodes' coordinate is (0,0,0), (1000,0,0) and (2000,0,0) repectively. And the model has 2 elements which type is BarElement, the first one's nodes is n1 and n2, and the second one's nodes is n2 and n3. The two elements is under element uniform load which is directed to the negative Y axis in global coordinate system. When the first and the last node is under fiexed constraints, the model will output an error named "Matrix must be symmetric positive definite" as i said before. I can not figure out why this happend. But When all three nodes is under fixed constraints, the model can be solved and I can get the support reaction force,.

The whole code of model is showned below:

        UniformParametric1DSection uniformParametric1DSection = new UniformParametric1DSection(
            6400, 1706670, 1706670);
        UniformIsotropicMaterial uniformIsotropicMaterial = new UniformIsotropicMaterial(
            9350, 0.3);

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

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

        Element element1 = new BarElement(node1, node2)
        {
            Section = uniformParametric1DSection,
            Material = uniformIsotropicMaterial,
        };
        Element element2 = new BarElement(node2, node3)
        {
            Section = uniformParametric1DSection,
            Material = uniformIsotropicMaterial,
        };

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

        BriefFiniteElementNet.Model model = new BriefFiniteElementNet.Model();
        model.Nodes.Add(node1, node2, node3);
        model.Elements.Add(element1, element2);
        model.Solve_MPC();
epsi1on commented 1 year ago

the problem is DoF.Rx of node 2 which is not fixed to ground neither connected to other nodes. Note that for uniformParametric1DSection you have not defined J property, so the BarElement will not have Shaft behaviour thus Rx of node2 is the problem. you can either fix the Rx of node 2, or simply add an extra value for j like this (it will be 1.0 for following code):

UniformParametric1DSection uniformParametric1DSection = new UniformParametric1DSection( 6400, 1706670, 06670,1);

or

UniformParametric1DSection uniformParametric1DSection = new UniformParametric1DSection( 6400, 1706670, 06670);
uniformParametric1DSection.J = 1;
epsi1on commented 1 year ago

Also pleas have a look at here, you can find out which node is the problem of PosDef exception https://github.com/BriefFiniteElementNet/BriefFiniteElement.Net/wiki/How-to-fix-NotPosDef-error

GIfengxin commented 1 year ago

Thank you for your timely response, I have solve the problem after assigning the J property of the section object. But I encounter another error when I want return the internal force of the BarElement object. What I want is maximum bending moment and shear along the element, but I don't find a ready-made function, so I have to use a for loop function to tranverse the entire element object. The function used is GetExactInternalForceAt() by assigning IsoPoint parameter. When IsoPoint is equal to -1 or 1, the function will return an InvalidInternalForceLocationException error, I can move the point a tiny distance just like -0.9999 or 0.9999, but when IsoPoint is equal to -0.6666666(-2/3), the same error returned again. When I want to know is that how many descrete points there are and how can I avoid the points. Thank you.

epsi1on commented 1 year ago

Thank you for your timely response, I have solve the problem after assigning the J property of the section object. But I encounter another error when I want return the internal force of the BarElement object. What I want is maximum bending moment and shear along the element, but I don't find a ready-made function, so I have to use a for loop function to tranverse the entire element object. The function used is GetExactInternalForceAt() by assigning IsoPoint parameter. When IsoPoint is equal to -1 or 1, the function will return an InvalidInternalForceLocationException error, I can move the point a tiny distance just like -0.9999 or 0.9999, but when IsoPoint is equal to -0.6666666(-2/3), the same error returned again. When I want to know is that how many descrete points there are and how can I avoid the points. Thank you.

There is a method called BarElement.GetInternalForceDiscretationPoints() which give you ISO points which you should bypass. actually i think you should have two points at each discrete point, one a little before that point, and one a little after that point, because some loads, like ConcentratedLoad can affect the internal force, so you will have different shears a little before and a little after the concentrated load apply point.

GIfengxin commented 1 year ago

Thank you for your reply. I use the method as you mentioned as below: IsoPoint[] discreteIsoPoints = (element1 as BarElement).GetInternalForceDiscretationPoints().. And this method return two IsoPoint which are {-1,0,0} and {1,0,0}, respectively. But when I set one IsoPoint with the Xi property is equal to -2/3(As I mentioned in the first question), the program return an error telling me that this IsoPoint is a discrete point. Is than I used a wrong method? Why the GetInternalForceDiscretationPoints didn't return the IsoPoint who's Xi property is equal to -2/3? Thank you.

epsi1on commented 1 year ago

I'll check it, there could be a bug. I'll reply here later

epsi1on commented 1 year ago

I use the method as you mentioned as below: IsoPoint[] discreteIsoPoints = (element1 as BarElement).GetInternalForceDiscretationPoints().. And this method return two IsoPoint which are {-1,0,0} and {1,0,0}, respectively. But when I set one IsoPoint with the Xi property is equal to -2/3(As I mentioned in the first question), the program return an error telling me that this IsoPoint is a discrete point.

Is the error for same model as code above?

epsi1on commented 12 months ago

inactivity