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
154 stars 57 forks source link

Non linear Analysis #6

Closed Loceno closed 2 years ago

Loceno commented 6 years ago

Hello @epsi1on, thank you for this powerfull library! I study FEM for my phd programme so i found this great work. In particular I'm focusing on the fire resistance of the structures applying FSE and I think that I could use this library to perform a parametric analysis but I would know if it's possible to take in account second order effects and the material's non linearity. I developed in C# a tool which build the moment-curvature diagram for steel-concrete composite sections in fire condition, I would use it to evaluate a "Young modulus" to perform a step by step analysis. You think that is it possible ?

epsi1on commented 6 years ago

Hello Loceno, Unfortunately no nonlinear analysis API is implemented into this framework. I would recommend you to use OpenSEES as i believe you can take both geometrical and material non linearity into account. That is in C++, and there are some official tutorials in its source code on developing new classes for elements and materials if it already its code doesn't exists.

I developed in C# a tool which build the moment-curvature diagram for steel-concrete composite sections in fire condition, I would use it to evaluate a "Young modulus" to perform a step by step analysis. You think that is it possible ?

I am not familiar with nonlinear analysis and solving procedures of it, so i can not confirm or deny that step by step nonlinear analysis is possible. Can you please describe how it works (step by step analysis)?

Thanks

Loceno commented 6 years ago

Hello @epsi1on , I apologize for late answer but I had some trouble. I know that opensees let take in account non linearities, but I'm interested at this fem core because it's fast and more flexible. As regard my analysis, I pubblished a paper about this topic and if you want I'll send it to you. In summary, my procedure: • starting from the mesh of the section and the thermal analysis we can define a mechanic model based on the uniaxial material behaviour described by the constitutive law ( in this case the ramberg osgood constitutive law of the Eurocode 1992 p1-2 for materials in fire condition); • Fix timestep; • Fix Axial load; • Fix curvature; • Evaluating the stress field in plane section hypotesis. In case of beam It's nearly true also in fire condition. The mechanical stress filed is evaluated by subtracting the from the total strain the thermal strain (e_sigma = e_tot - e_th); • variation of the neutral axis position in order to obtain the equilibrium with the axial load; • Evaluation of the resultant moment of stress field; • repeat for all curvature; • Repeat for all axial load up to the limit axial load; • Repeat for all timestep. The results of the analysis is a set of moment curvature diagram, that I would use in order to evaluate the stiffness of each beam element, and change it for each step of the analysis (In this moment I think about the anisotropic material properties in order to change only the young modulus that refers to bending). In addition in these days i tried to add second order effects by an iterative way and I made it. I think that is really simple to add to the core of the software but my programming skill aren't good as yours. the code is here https://pastebin.com/bAgctA99 As regard the mechanical non linearities, I think that is possible to implement on the non linear analysis by adding the internal distorsions and (first of all) adding the solver LU instead of the Cholesky decomposition. I hope to make a positive contribution.

epsi1on commented 6 years ago

Thanks for detailed information. I have a question about geo & mat NL analysis. If we have a model with n nodes, then total number of DoFs will be N = 6*n. if we assume a trial displacement vector named Δtrial with size N x 1 then force vector (we name F vector) with same size as delta can be calculated;

Ftrial = K * Δtrial this equation is for linear analysis and K is constant. Can you please give a brief comment on how to calculate Ftrial if analysis is nonlinear (mat and geo nonlinearity)?

Kind Regards

Loceno commented 6 years ago

I think that for the nonlinear analysis you have to update the stiffness matrix in a step by step calculus. In case of Geometric non linearity you can apply the iterative procedure that vary the geometry of structure (Like I showed before), https://www.google.com/search?q=second+order+effects&client=firefox-b&source=lnms&tbm=isch&sa=X&ved=0ahUKEwjrpLjm8v3bAhXQ2KQKHYvpDJ4Q_AUICigB&biw=1638&bih=785#imgrc=XgKZwFfYLvy_5M: . In order to take in account the material's NL, You have to update the stiffness matrix component in order to take in account the plastic part of constitutive law. For example in structures with beam elements where the material is caracterized by Elastic-perfect plastic behaviour you can perform the analysis adding Hinges at model in nodes where the Mpl is attained. By this way you can perform elastic analysis at each step with another configuration with a ΔF and take in account material non linearity. I developed a class which perform a geometric non linear analysis (wrote better than last time) . https://pastebin.com/1f4MkmnT there is an example of implementation

public FrameElement2Node assign(FrameElement2Node el ) { var S1 = BriefFiniteElementNet.SectionGenerator.GetRectangularSection(0.2, 0.1); el.Geometry = S1; double E = 210E9, ni = 0.3; el.E = E; el.G = E / (2 * (1 + ni)); el.UseOverridedProperties = false; return el; }

private void example() { secondOrder solver = new secondOrder(); var model = new Model(); Node n1 = new Node(0, 0, 0); Node n2 = new Node(0, 0, 4); n1.Constraints = Constraints.Fixed; n2.Loads.Add(new NodalLoad(new Force(10000, 0, -2000000,0,0,0))); model.Nodes.Add(n1, n2); FrameElement2Node el = new FrameElement2Node(model.Nodes[0], model.Nodes[1]); el = assign(el); model.Elements.Add(el); model = solver.BeamStructuresMesh(model,0.2); model.Solve(); double my = (model.Elements[0] as FrameElement2Node).GetInternalForceAt(0).My; solver.SetNodeLog(1); model = solver.RunSecondOrderAnalysis(model, 500); }

Soon I'll build an example for a step by step analysis with elastoplastic material. I have one question about this topic: can you explain to me how can I use the uniform anisotropic material properties with beam elements ?

P.S.: BFE.net is fantastic!

epsi1on commented 6 years ago

I have one question about this topic: can you explain to me how can I use the uniform anisotropic material properties with beam elements ?

Can you please describe what you mean by anisotropic material? you mean different Young modulus in 3 direction? or you mean nonuniform material along beam's length? Anyways i think you should use BarElement instead of FrameElement2Node. Also please note that bar element is under development.

Thanks

Loceno commented 6 years ago

I mean different young modulus in 3 direction because I would assign different stiffness for bending moment and axial loads. Why should I use the bar element ?

epsi1on commented 6 years ago

Sorry for delayed response. As far as i know only young modulus in local X direction is used for both axial load and bending moment (young modulus in local X direction) as both bending and axial load cause stress in local X direction.

i suggested to use bar element because that supports both nonuniform material and geometric section of bar element. example of nonuniform geometric section is a tapered beam.

Loceno commented 6 years ago

young modulus in local X direction is used for both axial load and bending moment (young modulus in local X direction) as both bending and axial load cause stress in local X direction.

the tensor has the 3 Young modulus for each direction, I supposed that I can use Ez for elongation and Ex and Ey for bending stresses. why bar element isn't fuly implemented yet ? what is missing ?

epsi1on commented 6 years ago

A BarElement can act as a truss or single direction beam or a shaft, or a frame (combination of these all). BarElement.Behaviour defines this. (more info) I'm not sure but maybe elognation is taken into account in Timoshenko theory for beams, but currently BarElement only support Euler-Bernoulli theory for beam behavior that only uses Second area moment of section, and Ex for all calculations. BarElement is missing Timoshenko behavior also partial end release.

epsi1on commented 5 years ago

Hi @Loceno, I have a question. Please correct me if i say something wrong. In solving models, we can define F vector as a function of Δ like this: F = G(Δ) in linear system: G(Δ) = K . Δ but in nonlinear system (both material and geometrical) this is not true anymore and relationship between F and Δ is not simply possible to be written as a direct formula. how we can computer F vector for a given test Δ vector in simplest form? The reason i am asking this is this: imagine i have a very sophisticated magical program for finding root of ANY nonlinear function, say we have a function named H(x), where x is a vector, the magic root finder will find us a x0 vector that causes the H(x0) = 0 and x0 is root of our function. This is how magic root finder works: When this software starts, it does give me a test Δ0 vector and ask me what is F(Δ0), then it does some calculations with itself, then gives me a Δ1 vector and ask me the F(Δ1) and so on ... finally in step N where it gives me ΔN, F(ΔN) would be very near to zero and i accept the ΔN as root of my equation system. For linear system, H(Δ) is equal to K.Δ-F, but in nonlinear form i do not know what is that, but need to know how to calculate it...

Thanks for reading this comment

Loceno commented 5 years ago

I think that you are using an algorithm like line search or newton-raphson in multidimensional form. In order to calculate the function you have to perform a step by step analysis. Fixing the initial value of F near 0 in the first step, you can linearize the problem so you can find the zero of H(Δ) as you do in linear systems. At each step you can refresh the stiffness matrix introducing the material and geometrical non linearities. As reagard the material non linearity it can be taken in account by varying the carachteristics of materials, in case of monodimensional beams you can change the stiffness of element (EA, EIx or EIy or only the elastic modulus in one direction). The non linear geometry can be taken in account by adding some quantities to the elements of stiffness matrix. Applying this strategy you can raplace at the H(Δ) which has a non linear behaviour, a sum of h(Δ) with a linear behaviour. By the point of view of Forces the equation F = G(Δ) is an equilibrium between the external force and the internal force of system which you calculate step by step. By the energetic point of view you can write the equilibrium as a sum of the energy of external force Eext the internal energy (theorem of Clapeyron) Eint and the dissipated energy Edis There you can find how to change the matrix during the calculation for second order effects.

I hope that this help you.

Alberto.

epsi1on commented 5 years ago

Thanks, i have some little knowledge about finite element method, but my knowledge about finite element theory is not as good as you. The scenario is that if i have a solver that uses an unknown method, and only gives me a Δ vector each time and ask me H(Δ) of that, how can i directly calculate it. i think i still can use G(Δ) = K . Δ for nonlinear scenario, but K is not constant, but is also depended to Δ, in other words: G(Δ) = K(Δ) . Δ and K(Δ) = Kmat + Kgeo hope it be right :)

Loceno commented 5 years ago

Hello @epsi1on, I'm using the library in order to perform nonlinear Analysis. In order to Take in account the non linear geometry I'm changing the position of nodes at each iteration of the calculus. To do this, I'm using two model: one contains the undeformed structure and the other contains the deformed Structure. I noticed that if I change the position of one node in one model the library refresh also the undeformed model. Can you tell me how to lock the undeformed model ?

epsi1on commented 5 years ago

Hello, I think you should create a clone of model. Assume original model is origModel and the model with nodal deformation at itteration is defModel : var defModel = origModel.Clone(); instead of: var defModel = origModel;

Loceno commented 5 years ago

Thank You for the Istant answer, I'm declaring the model as Model model1 = new Model(); because I'm developing some classes. Does it work anyway ?

epsi1on commented 5 years ago

Thanks for information. Can you please give more info on how do you declare the second model? I think that on is more important.

Thanks

Loceno commented 5 years ago

I didn't know the Method Clone() so I wrote a method wich clone the element by element and node by node.

Thank you

epsi1on commented 5 years ago

As far as I can understand if you use Model.Clone() method, you will have a complete separated model without any dependency to original model. Can you please use Model.Clone() and check the result? Thanks

Loceno commented 5 years ago

It Return me a ArgomentNullException with the message "The graphic object can not be null."

epsi1on commented 5 years ago

OK, there should be something null in the model, and it causes the error "Object graph cannot be null". I cannot recreate this error on my side, so can you please give me a model that gives same error when you call .Clone() method of it?

Thanks

Loceno commented 5 years ago

Hello @epsi1on , I made a mistake, the method work fine! I only had to refresh my code. there's something new about the library ? Have you implemented the non linear analysis yet ? :D Btw I would be very gald to have you among my linkedin contacts. Alberto Compagnone

epsi1on commented 5 years ago

Hi @Loceno, I've done some tries for implementing geometric nonlinear analysis based on stress stiffening theory, (actually I've ported the geometric non linearity from this code: https://sourceforge.net/p/frame3dd/code/HEAD/tree/trunk/src/frame3dd.c). It seems works fine, but there is some problems in converging if structure is regular, it means nodes are exactly on a grid. But if model have a little irregularity it do converge better for first step. I think it needs more effort. Anyways Just wanted to kindly ask you please sending me the paper of yours about nonlinear analysis. Thanks

ahmadizm commented 5 years ago

Hi @epsi1on,

The solver behavior you mentioned can be expected, because when the elements are initially perfect, they can be theoretically loaded beyond their stability limits. In these situations, or even when the loads are close to their instability limits, any subsequent lateral movement (resulting from a load, connected element, or even numerical round-off errors) can trigger a sudden buckling-like behavior, which will then result in large changes of stiffness. Such significant changes of behavior cannot be handled by most solvers. This can also occur even when geometry is not updated if the geometric stiffness fully cancels out the initial stiffness in some degrees of freedom.

To address these issues, usually special solving procedures need to be used (such as arc-length or work-control), or elements must have some small initial out-of-straightness (like 0.001xLength for a beam with intermediate nodes).

Good luck to you, you are doing some interesting things.

M.

rubsy92 commented 4 years ago

Hi @epsi1on ,

Where did you put the nonlinear code? How should I perform a non-linear analysis? I might check it out.

Thanks!

epsi1on commented 4 years ago

Where did you put the nonlinear code? How should I perform a non-linear analysis? I might check it out.

This project is dedicated for linear analysis. Nonlinear would be another project if it goes for existence. I've just have stuck at implementing the arc-length method, any help would be appreciated.

epsi1on commented 4 years ago

I've found this project, that definitely is interesting, implements nonlinear solver like ArcLength... https://github.com/EduardBargues/NonLinearEquationsSolver

rubsy92 commented 4 years ago

I've found this project, that definitely is interesting, implements nonlinear solver like ArcLength...

Interesting library! Would it be a lot of work to implement it or does it need some adjustments first?

epsi1on commented 4 years ago

I've found this project, that definitely is interesting, implements nonlinear solver like ArcLength...

Interesting library! Would it be a lot of work to implement it or does it need some adjustments first?

Indeed lot of work. I've made new repository for nonlinear analysis here https://github.com/epsi1on/NonlinearBFE

epsi1on commented 2 years ago

continue this topic in discussion section: https://github.com/BriefFiniteElementNet/BriefFiniteElement.Net/discussions/104

epsi1on commented 2 years ago

Hello epsi1on, it is possible to perform a dynamic analysis. You can add an example where the use of Slab is included and how to discretize it using the tool.

Sorry can you please describe? I didn't get it...

epsi1on commented 2 years ago

Hi,

There are two type of 2D element in the library, TriangleElement and QuadrilaturalElement.

The TriangleElement is tested and verified, also more reliable than QuadrilaturalElement so i suggest you to use TriangleElement for now.

Splitting a rectangle into smaller triangles is not something hard, first you need to make the node grid, and then generate triangles. but for more complex shapes you can use libraries who implement Deluanay triangulation (for example Triangle.NET). for example pls. take a look at this:

https://github.com/BriefFiniteElementNet/BriefFiniteElement.Net/tree/master/BriefFiniteElementNet.Validation/Case_02

image

epsi1on commented 2 years ago

Can you please check the UniformLoad and see if it works? There is example for uniform load on beam here. Same is for triangle element