Metoda elementelor finite (FEM) este o metodă de aproximare a variabilelor care sunt funcții complicate de dimensiuni spațiale folosind funcții mai simple pe bucăți, fiecare valabilă într-un element finit mic. Figura 2 prezintă un domeniu de interes V și unele regiuni ale acestuia discretizate în elemente finite. Un element finit este de obicei un poliedru convex cu fețe, muchii și vârfuri (numite noduri) care îi constituie limitele. Distribuția unei variabile în cadrul unui element finit, la un moment dat, se obține prin interpolarea valorilor nodale la acel moment.
The finite element method (FEM) is a method of approximating variables that are complicated functions of spatial dimensions using simpler piecewise functions each valid within a small finite element. Figure 2 shows a domain of interest V, and some regions of it discretized into finite elements. A finite element is typically a convex polyhedron with faces, edges and vertices (called nodes) constituting its boundaries. The distribution of a variable within a finite element, at a given time, is obtained by interpolating the nodal values at that time.

Fig. 2 Descrierea deformării unui domeniu și discretizarea domeniului și deformarea acestuia. Este de asemenea prezentată relația dintre deplasare și deformare (și dintre variațiile acestora). Vectorii sunt indicați folosind caractere aldine atunci când un index este omis. Indicii inferiori indică componente, iar indicele special n indică nodul cu care este asociată o cantitate. Indicii repetați implică însumarea tuturor direcțiilor de coordonate sau peste toate nodurile, iar o virgulă în indicii inferiori indică diferențierea față de direcția (direcțiile) indicată de indexul (indici) ulterior.

Pentru elementele izoparametrice, locațiile originale (materiale) ale celor n noduri Xn prescriu geometria originală a elementului X(ξ) = XnNn(ξ) [Rețineți că ecuațiile sunt scrise folosind notația indicială, Belytschko și colab. (2000). În timp ce alți indici variază în funcție de dimensiunile spațiale, indicele inferior n identifică cantitățile nodale și variază în funcție de numărul de noduri ale unui element sau ale întregii structuri, în funcție de context, iar e variază în funcție de numărul de elemente. Simbolurile cu caractere aldine reprezintă vectori/tensori.] ca o mapare de la un element părinte (o instanță geometrică simplă a tipului de element) prin funcții de interpolare Nn în același mod ca și valorile nodale ale unei variabile (deplasările un în Fig. 1) guvernează variația sa spațială asupra elementului u(ξ,t) = un(t)Nn(ξ). Locația spațială curentă a materialului în locația inițială X(ξ) este dată de x(ξ, t) = u(ξ, t) + X(ξ) = xn(t)Nn(ξ). Funcțiile de interpolare Nn(ξ), numite și funcții de formă sau coeficienți de influență, exprimă influența relativă a fiecăruia dintre cele n noduri ale unui element asupra variabilelor din fiecare punct din cadrul acestuia ca funcție de coordonatele sale locale ξ, dar pot fi exprimate și în termeni de coordonatele materialelor X ca Nn(ξ(X)). Rețineți că coeficienții de influență într-un punct ξ sunt aceiași pentru toate variabilele (deplasări, viteze, accelerații etc.) la un nod. În analiza termomecanică cuplată, variabilele nodale includ și temperatura.
For isoparametric elements, the original (material) locations of the n nodes Xn prescribe the original geometry of the element X(ξ) = XnNn(ξ) [Note that the equations are written using indicial notation, Belytschko et al. (2000). While other indices vary over spatial dimensions, subscript index n identifies nodal quantities and varies over the number of nodes of an element or the whole structure, depending on the context, and e varies over the number of elements. Boldface symbols represent vectors/tensors.] as a mapping from a parent element (a geometrically simple instance of the type of element) via interpolation functions Nn in the same manner as the nodal values of a variable (displacements un in Fig. 1) govern its spatial variation over the element u(ξ, t) = un(t)Nn(ξ). The current spatial location of the material at original location X(ξ) is given by x(ξ, t) = u(ξ, t) + X(ξ) = xn(t)Nn(ξ). The interpolation functions Nn(ξ), also called shape functions or influence coefficients, express the relative influence of each of the n nodes of an element on the variables at each point within it as a function of its local coordinates ξ but can also be expressed in terms of the material coordinates X as Nn(ξ(X)). Note that the influence coefficients at a point ξ are the same for all the variables (displacements, velocities, accelerations, etc.) at a node. In coupled thermomechanical analysis, the nodal variables also include temperature.

Dacă se aplică sarcini externe și/sau condiții de limită peste unele puncte, linii, arii sau volume ale domeniului de interes, discretizarea domeniului se realizează într-o manieră care să simplifice descrierea acestora; de exemplu, dacă se aplică sarcini concentrate în anumite puncte, nodurile ar fi situate în acele puncte. Dacă acest lucru este exclus, de exemplu, de locația acestor sarcini care se deplasează în raport cu materialul în timp, sarcinile pot fi aproximative sau ochiurile de rețea ar fi modificate în timp pentru a urma sarcinile, iar fluxul de material prin elemente ar fi urmărite.
If external loads and/or boundary conditions are applied over some points, lines, areas, or volumes of the domain of interest, the discretization of the domain is accomplished in a manner as to simplify the description of these; for instance, if concentrated loads are applied at some points, nodes would be situated at those points. If this is precluded, for instance, by the location of these loads moving with respect to the material over time, the loads may be approximated or the mesh would be modified over time to follow the loads, and the flow of material across elements would be tracked.

Funcțiile de interpolare utilizate pentru aproximarea distribuției spațiale a variabilelor în cadrul fiecărui element finit sunt de obicei polinoame. Funcțiile de interpolare sunt de obicei parametrizate în termeni de valori nodale ale acestor variabile așa cum s-a menționat mai sus. Pentru un set dat de valori nodale (adică, parametrii), polinomul este complet definit peste granițe, precum și interiorul elementului (adică, toate gradele de libertate ale polinomului sunt constrânse). Din ordinul polinomului se poate deduce ordinul continuității de-a lungul granițelor. În schimb, ordinul polinomului, ales pentru aproximarea variabilei depinde de cerințele de continuitate (netezime) peste granițele elementelor iar numărul de parametri, adică de grade de libertate, va crește odată cu ordinul polinomului – numărul de noduri și/sau numărul de variabile definite la noduri va crește (se pot adăuga panta, curbura etc.). Cerințele de continuitate sunt impuse prin recunoașterea faptului că valorile variabilelor de guvernare la un nod sunt aceleași pentru fiecare dintre elementele care împart nodul, în timpul asamblării ecuațiilor care guvernează comportamentul elementelor individuale într-un sistem de ecuații reprezentând comportamentul întregul domeniu de interes.
The interpolation functions used for approximating the spatial distribution of the variables within each finite element are typically polynomials. The interpolation functions are typically parameterized in terms of the nodal values of these variables as mentioned above. For a given set of nodal values (i.e., the parameters), the polynomial is fully defined over the boundaries as well as the interior of the element (i.e., all the degrees of freedom of the polynomial are constrained). From the order of the polynomial, the order of continuity along the boundaries can be inferred. Conversely, the order of the polynomial chosen for approximating the variable depends on the continuity (smoothness) requirements across element boundaries, and the number of parameters, i.e., degrees of freedom, will increase with the order of the polynomial – the number of nodes and/or the number of variables defined at the nodes will increase (slope, curvature, etc., may be added). The continuity requirements are enforced by recognizing that the values of the governing variables at a node are the same for each of the elements sharing the node, during assembly of the equations governing the behavior of the individual elements into a system of equations representing the behavior of the entire domain of interest.

Deformarea mecanică

Ideea fundamentală din spatele FEM este că valorile nodale care caracterizează variația spațială a cantităților necunoscute pot fi ușor de ajuns prin exprimarea principiilor care guvernează comportamentul domeniului într-o formă „concentrată”, adică ca principiu integral echivalent, că este satisfăcută pentru fiecare element finit. În mod obișnuit, principiile sunt exprimate ca ecuații diferențiale parțiale (PDE) valabile în fiecare punct din domeniu, numite formă „puternică”. Ele pot fi, de asemenea, exprimate ca o formă integrală echivalentă, care se numește forma „slabă”, deoarece satisfacerea formei integrale pentru un set de funcții de testare nu implică neapărat că PDE este satisfăcută peste tot în cadrul elementului. Pentru deformarea pur mecanică, ecuația impulsului (sau ecuația mișcării, Ec. 1) care trebuie să fie satisfăcută în fiecare punct din domeniu este forma puternică, iar principiul lucrului virtual (Ec. 2) este o formă slabă. Belytschko și colab. (2000, Sect. 2.3) arată derivarea acestei forme slabe din Ec. 1 prin integrarea produsului ecuației de echilibru cu o variație arbitrară a deplasării, funcția de test. Pașii cheie includ integrarea prin părți a termenului gradient de solicitare pentru a reduce ordinul derivatelor și recunoașterea faptului că produsul scalar al solicitării și variației gradientului de deplasare este același cu produsul scalar al solicitării și variației deformării, datorită simetriei tensorului de solicitare.
The fundamental idea behind FEM is that the nodal values that characterize the spatial variation of unknown quantities can be easily arrived at by expressing the principles that govern the behavior of the domain in a “lumped” form, i.e., as an equivalent integral principle, that is satisfied over each finite element. Typically, the principles are expressed as partial differential equations (PDEs) valid at every point within the domain, called the “strong” form. They can also be expressed as an equivalent integral form, which is called the “weak” form because satisfaction of the integral form for a set of test functions does not necessarily imply that the PDE is satisfied everywhere within the element. For purely mechanical deformation, the momentum equation (or the equation of motion, Eq. 1) that needs to be satisfied at every point within the domain is the strong form, and the principle of virtual work (Eq. 2) is a weak form. Belytschko et al. (2000, Sect. 2.3) show the derivation of this weak form from Eq. 1 by integrating the product of the equilibrium equation with an arbitrary variation in the displacement, the test function. Key steps include integration by parts of the stress gradient term to reduce the order of the derivatives and recognizing that the inner product of stress and variation in displacement gradient is the same as the inner product of the stress and variation in strain, due to the symmetry of the stress tensor.

Ecuația 1 afirmă că forța inerțială datorată ratei de schimbare a impulsului este egală cu suma forței corpului și a forței datorate gradientului de solicitare. Ecuația 2 arată că suma modificării energiei interne și a modificării energiei cinetice a structurii, datorită unei variații arbitrare a deplasării, ar trebui să fie egală cu suma lucrului extern efectuat de vectorii de tracțiune distribuiți ({τ} = τi), forțele corpului distribuite ({β} = βi) și vectorii de sarcină punctuală la n noduri (Pni).

Metoda diferențelor finite impune diferențiala parțială Ec. 1 în puncte de pe o grilă uniform distanțată, ceea ce impune unele limitări asupra geometriilor și condițiilor la limită. Pe de altă parte, FEM satisface principiul lucrului virtual pe o bază element cu element pentru geometrii arbitrare și condiții la limită, în timp ce satisface aproximativ ecuația mișcării.
Equation 1 states that the inertial force due to the rate of change of momentum equals the sum of the body force and force due to the stress gradient. Equation 2 states that the sum of the change in internal energy and the change in kinetic energy of the structure, due to an arbitrary variation in displacement, should equal the sum of the external work done by distributed traction vectors ({τ} = τi), distributed body forces ({β} = βi), and the point load vectors at the n nodes (Pni).

The finite difference method enforces the partial differential Eq. 1 at points on a uniformly spaced grid, which imposes some limitations on the geometries and boundary conditions. On the other hand, FEM satisfies the principle of virtual work on an element-by-element basis for arbitrary geometries and boundary conditions, while approximately satisfying the equation of motion.

După cum se arată în Fig. 2, deplasarea și deformarea în cadrul elementului pot fi exprimate ca funcții de deplasările nodale. Măsura deformării prezentată, εij, este tensorul de deformare infinitezimal lagrangian, care este partea simetrică a lui ui,j, gradientul de deplasare față de configurația originală (sau cadrul de referință al materialului). Variația deformării δεij poate fi exprimată ca parte simetrică a gradientului variației deplasării δui. Deoarece funcțiile de formă sunt sursa gradientului spațial de deplasare, gradientul de deplasare ui,j este un produs simplu al deplasărilor nodale cu gradientul funcțiilor de formă, ui,j= uîn Bnj, unde
As shown in Fig. 2, the displacement and strain within the element can be expressed as functions of the nodal displacements. The strain measure shown, εij, is the Lagrangian infinitesimal strain tensor, which is the symmetric part of ui,j, the displacement gradient with respect to the original configuration (or material reference frame). The variation in strain δεij can be expressed as the symmetric part of the gradient of the variation in displacement δui. Since the shape functions are the source of the spatial gradient of displacement, the displacement gradient ui,j is a simple product of the nodal displacements with the gradient of the shape functions, ui,j = uinBnj, where

Pentru problemele 3D, cele șase componente independente ale solicitării sunt scrise de obicei ca un vector {σ} = [σ11, σ22, σ33, σ12, σ23, σ31]^T (în forma Voigt). Pentru un element cu n noduri, dacă deplasările nodale sunt reformate ca un vector ^e{U} = [U1n, U2n, U3n]^T, de lungime 3n, matricea funcției de formă [N] pentru element este
For 3D problems, the six independent components of the stress are typically written as a vector {σ} = [σ11, σ22, σ33, σ12, σ23, σ31]^T (in Voigt form). For an element with n nodes, if the nodal displacements are recast as a vector ^e{U} = [U1n, U2n, U3n]^T, of length 3n, the shape function matrix [N] for the element is

de mărime 3 x 3n. Gradientul funcțiilor de formă [B] este dat de
of size 3 x 3n. The gradient of the shape functions [B] is given by

de mărime 6 x 3n care, înmulțită cu deplasările nodale, dă vectorul deformare {ε} = [ε11, ε22, ε33, γ12, γ23, γ31]^T. Rețineți că γij = εij + εji este utilizat în vectorul deformare, astfel încât {δε}^T{σ} oferă incrementul de lucru efectuat pe unitatea de volum. După cum s-a menționat mai sus, gradienții sunt de obicei calculați în raport cu coordonatele locale și transformați în cadrul de referință material folosind inversul jacobianului Jjk = δXj/δξk. Integrând termenii ecuației de lucru virtual peste volumul fiecărui element și însumând peste toate elementele, ecuațiile care trebuie îndeplinite pentru variația arbitrară în deplasările nodale sunt obținute așa cum se arată (folosind notația matriceală) în Ec. 3.
of size 6 x 3n that when multiplied with the nodal displacements gives the strain vector {ε} = [ε11, ε22, ε33, γ12, γ23, γ31]^T. Note that γij = εij + εji is used in the strain vector so that {δε}^T{σ} gives the increment of work done per unit volume. As noted above, gradients are usually calculated with respect to the local coordinates and transformed to the material reference frame using the inverse of the Jacobian Jjk = δXj/δξk. Integrating the terms of the virtual work equation over the volume of each element and summing over all elements, the equations to be satisfied for arbitrary variation in nodal displacements are obtained as shown (using matrix notation) in Eq. 3.

unde {U} este vectorul de deplasare globală, {P} este vectorul de încărcare globală, iar prefixul superior e indică cantități sau grade de libertate specifice elementului.
where {U} is the global displacement vector, {P} is the global load vector, and the superscript prefix e indicates element-specific quantities or degrees of freedom.

Solicitările sunt legate de deformații prin ecuații constitutive, de exemplu, elasticitatea liniară unde ^e{σ} = [D]^e{ε} = [D][B]^e{U}, unde [D] este matricea elasticității. Pentru modelele constitutive neliniare, de exemplu, plasticitatea, solicitarea trebuie urmărită separat și actualizată pentru a corespunde fiecărui increment în procesul de soluție (vezi mai jos). Deoarece vectorii deplasărilor nodale, accelerațiile și variațiile deplasării sunt constante, ei pot fi factorizați din integrale. Corpul distribuit și forțele de tracțiune sunt repartizate în mod corespunzător nodurilor prin funcțiile de formă care acționează ca funcții de ponderare, rezultând vectorii forțe nodale corespunzători ^e{b} și ^e{t}.
The stresses are related to the strains via constitutive equations, e.g., linear elasticity where ^e{σ} = [D]^e{ε} = [D][B]^e{U}, where [D] is the elasticity matrix. For nonlinear constitutive models, e.g., plasticity, the stress needs to be tracked separately and updated to correspond to each increment in the solution process (see below). Since the vectors of nodal displacements, accelerations, and variations in displacement are constants, they can be factored out of the integrals. Distributed body and traction forces are appropriately apportioned to the nodes by the shape functions acting as weighting functions, resulting in the corresponding nodal force vectors ^e{b} and ^e{t}.

Prima integrală dintre paranteze dă naștere matricei de rigiditate a elementului ^e[K] care, atunci când este înmulțită cu vectorul deplasare, dă vectorul forță internă ^e{f}^int necesar pentru a efectua deformația corespunzătoare deplasării ^e{U}. A doua integrală dă matricea masei elementului ^e[M] care, atunci când este înmulțită cu vectorul accelerație, dă forța de inerție necesară pentru a provoca accelerația. Integrările sunt efectuate numeric, de obicei folosind cuadratura Gauss care evaluează integrala ca o sumă ponderată a valorilor la funcțiile la unul sau mai multe puncte de integrare în cadrul elementului. Numărul de puncte de integrare și ponderile depind de ordinul total al polinoamelor din integrand. Integrarea redusă se referă la integrarea numerică folosind un număr mai mic de puncte de integrare decât este necesar pentru integrarea exactă a polinomului pentru a accelera calculele, precum și pentru a atenua răspunsul prea rigid (blocarea elementelor complet integrate, cap. 8 din Belytschko et al. ( 2000)). Insă, acest lucru face ca deformarea să fie calculată incorect ca fiind zero pentru unele moduri de deformare (moduri de deformare cu energie zero, numite și moduri de clepsidră deoarece aceasta este forma acestor moduri pentru elementele patrulatere) și ar trebui completată cu tehnici de control al clepsidrei care atribuie o forță de penalizare pentru a le atenua. Dacă este raportată o cantitate semnificativă de energie de clepsidră, aceasta indică necesitatea de a rafina rețeaua pentru a lua în considerare corect energia de deformare în regiunile cu gradienți mari de deformare. The first integral within the brackets gives rise to the element stiffness matrix ^e[K] that when multiplied by the displacement vector gives the internal force vector ^e{f}^int required to effect the deformation corresponding to the displacement ^e{U}. The second integral gives the element mass matrix ^e[M] that when multiplied by the acceleration vector gives the inertial force required to cause the acceleration. The integrations are performed numerically, typically using Gauss quadrature that evaluates the integral as a weighted sum of the values at the functions at one or more integration points within the element. The number of integration points and the weights depends upon the total order of the polynomials in the integrand. Reduced integration refers to numerical integration using a smaller number of integration points than required for accurate integration of the polynomial to speed up the calculations as well as to mitigate overly stiff response (locking of fully integrated elements, Chap. 8 of Belytschko et al. (2000)).
However, this causes the strain to be incorrectly calculated as being zero for some modes of deformation (zero-energy modes of deformation, also called hourglass modes because this is the shape of these modes for quadrilateral elements) and should be complemented by hourglass control techniques that assign a penalty force to mitigate these. If significant amount of hourglass energy is reported, it indicates the need to refine the mesh to correctly account for the deformation energy in regions of high strain gradients.

Pentru un element cu n grade de libertate, ecuația 5 dă naștere la n termeni ai vectorilor de tracțiune și forței corpului și la nxn termeni ai matricelor de rigiditate și masă care reprezintă interacțiunile dintre cele n grade de libertate. Termenii corespunzători unui grad de libertate la un nod vor figura în ecuațiile pentru fiecare dintre elementele din care face parte nodul. Procesul de adunare a vectorilor de deplasare a elementului într-un vector de deplasare global {U} se numește asamblare. În timpul asamblării, termenii forță de tracțiune și forța corpului corespunzători unui grad de libertate, care decurg din diferitele elemente din care face parte un nod, sunt însumați pentru a obține forțele globale ale corpului și de tracțiune de-a lungul acelui grad de libertate. În mod similar, elementele care împart granițe dau naștere la diferite matrici de rigiditate și de masă ale elementelor, fiecare conținând termeni reprezentând interacțiunea dintre aceleași perechi de grade de libertate, care trebuie însumate pentru a obține matricele globale de rigiditate și masă. Rezultatul final al asamblării este următoarea ecuație globală:
For an element with n degrees of freedom, equation 5 gives rise to n terms of the traction and body force vectors, and to nxn terms of the stiffness and mass matrices that represent the interactions between the n degrees of freedom. Terms corresponding to a degree of freedom at a node will figure in the equations for each of the elements that the node is part of. The process of collating the element displacement vectors into a global displacement vector {U} is called assembly. During assembly, the traction force and body force terms corresponding to a degree of freedom, arising from the different elements that a node is part of, are summed to get the global body and traction forces along that degree of freedom. Similarly, elements that share boundaries give rise to different element stiffness and mass matrices, each containing terms representing interaction between the same pairs of degrees of freedom, that need to be summed together to obtain the global stiffness and mass matrices. The end result of assembly is the following global equation.

Mărimea dintre paranteze pătrate este vectorul forță reziduală. Recunoscând că produsul său scalar cu variația arbitrară în deplasările nodale va fi zero numai dacă fiecare dintre termeni este zero, ecuațiile sistemului sunt scrise ca:
The quantity within square brackets is the residual force vector. Recognizing that its inner product with arbitrary variation in nodal displacements will be zero only if each of the terms is zero, the system equations are written as

Ecuația 7 de mai sus este o semi-discretizare cu elemente finite, deoarece, deși dependența spațială a fost discretizată prin utilizarea elementelor finite, dependența de timp este încă sub forma unei ecuații diferențiale. Rețineți că dacă forțele vâscoase ar fi fost incluse și în ecuația mișcării, ar exista un termen suplimentar care depinde de vitezele nodale. Condițiile inițiale și la limită trebuie aplicate, iar acest sistem de ecuații diferențiale ordinare (ODE) de ordinul doi trebuie rezolvat.
Equation 7 above is a finite element semi-discretization, because though the spatial dependence has been discretized through the use of finite elements, the time dependence is still in the form of a differential equation. Note that If viscous forces had also been included in the equation of motion, there would be an additional term dependent on the nodal velocities. The initial and boundary conditions need to be applied, and this system of second order ordinary differential equation (ODEs) needs to be solved.

Pentru încărcarea cvasistatică, cum ar fi atunci când se analizează tăierea unui material perfect plastic în scopul dezvoltării modelelor de câmp cu linii de alunecare, forța de rigiditate domină și ecuația sistemului poate fi scrisă ca [K] {U} = {f}. Dacă translațiile și rotațiile corpului rigid nu sunt constrânse, structurile sunt capabile de mișcări ale corpului rigid care necesită energie internă zero. Deoarece acestea sunt moduri proprii cu valoare proprie zero (deplasare finită care necesită forță zero), iar determinantul unei matrice este produsul valorilor sale proprii, matricea completă de rigiditate este singulară sau insuficientă ca rang. Deoarece un grad de libertate este constrâns, termenii de forță corespunzători acestuia pot fi mutați la termenul de forță externă din partea dreaptă, eliminând coloanele matricei de rigiditate corespunzătoare acelui grad de libertate. Rândul poate fi, de asemenea, eliminat deoarece este util să se calculeze doar forța de reacție odată ce celelalte deplasări sunt cunoscute. Când un număr suficient de aceste grade de libertate este constrâns, matricea de rigiditate devine nesingulară și se poate rezolva și problema cvasistatică (problema dinamică poate fi întotdeauna rezolvată).
For quasistatic loading, such as when analyzing the cutting of a perfectly plastic material for the purpose of developing slip-line field models, the stiffness force dominates and the system equation can be written as [K] {U} = {f}. If rigid body translations and rotations are not constrained, structures are capable of rigid-body motions requiring zero internal energy. Since these are eigenmodes with zero eigenvalue (finite displacement requiring zero force), and the determinant of a matrix is the product of its eigenvalues, the full stiffness matrix is singular or rank deficient. As a degree of freedom is constrained, the force terms corresponding to it can be moved to the external force term on the right side, eliminating the columns of the stiffness matrix corresponding to that degree of freedom. The row can also be eliminated as it is only useful to compute the reaction force once the other displacements are known. When a sufficient number of these degrees of freedom are constrained, the stiffness matrix becomes non-singular and the quasistatic problem can also be solved (the dynamic problem can always be solved).

Integrare explicită și implicită în timp

Cele două metode principale de soluție sunt (i) metodele de integrare directă și (ii) metodele modale. Metodele de integrare directă folosesc diferența finită pentru a rezolva ec. 7 și urmăresc deplasările nodale în timp. Există două clase de aproximări cu diferențe finite la i+1{U} în termeni de deplasări i{U}, viteze și accelerații în momente succesive de timp, it.
The two main methods of solution are (i) direct integration methods and (ii) modal methods. Direct integration methods use finite difference to solve Eq. 7 and track the nodal displacements over time. There are two classes of finite difference approximations to i+1{U} in terms of the displacements i{U}, velocities, and accelerations at successive instants of time, it.

Principala diferență dintre metodele implicite și cele explicite este că, în metodele explicite, deplasarea i+1{U} în următorul moment de timp i+1t depinde numai de deplasarea, viteza și accelerația la momentele curente și anterioare și poate fie pur și simplu calculată dacă acestea sunt cunoscute. Acest lucru se realizează prin utilizarea ecuațiilor sistemului la momentul it pentru a calcula viteza i+1/2{Ú} la i+1/2t și deplasarea i+1{U} la i+1t, așa cum se arată mai jos.
The main difference between implicit and explicit methods is that in explicit methods, the displacement i+1{U} at the next instant of time i+1t depends only upon the displacement, velocity, and acceleration at the current and previous instants and can be simply calculated if these are known. This is achieved by using the system equations at time it to calculate the velocity i+1/2{Ú} at i+1/2t and the displacement i+1{U} at i+1t, as shown below.

Rețineți că deplasarea i+1{U} este o extrapolare; dacă vitezele și accelerațiile continuă să rămână apropiate de ceea ce sunt în momentul curent, atunci deplasarea la următorul moment de timp poate fi bine aproximată prin aplicarea ecuației mișcării la momentul curent. Astfel, dimensiunea pasului de timp Δt = (i+1t - it) trebuie să fie mică în comparație cu perioada de timp chiar și a celor mai înalte oscilații de frecvență din sistem pentru ca deplasarea estimată i+1{U} să nu înceapă să se abată de la deplasarea efectivă. Mai exact, Δt ≤ 2/ωmax, care este echivalent cu Δt ≤ L/c unde L este cea mai scurtă latură a celui mai mic element și c este viteza sunetului (viteza undei de dilatație) în element. Din acest motiv, partea cea mai mică a elementului controlează dimensiunea pasului de timp, iar elementele foarte mici sau foarte distorsionate duc la timpi mari de soluție. Pentru modelarea exactă a deformației din jurul muchiei de tăiere, de exemplu, pentru a modela cu acuratețe tensiunile reziduale, delaminarea subsuprafeței etc., sunt necesare elemente foarte mici, ceea ce face ca pașii de timp să fie mici și timpul de simulare să fie mare. Subciclarea, prin care forțele nodale datorate fiecăruia dintre elemente sunt actualizate la diferiți multipli ai pasului de timp critic al celui mai mic element este uneori utilizată, dar nu este foarte comună. În afară de aceasta, o modalitate eficientă de a verifica dacă este utilizat un set consistent de unități pentru deplasări, sarcini, solicitări etc., este de a verifica dacă viteza undelor de dilatație (Vc = √E/ρ, unde E este modulul elastic și ρ este densitatea) se obține corect din valorile lui E și ρ.
Note that the displacement i+1{U} is an extrapolation; if the velocities and accelerations continue to remain close to what they are at the current time, then the displacement at the next instance of time can be approximated well by enforcing the equation of motion at the current time. Thus, the time step size Δt = (i+1t - it) needs to be small compared to the time period of even the highest frequency oscillations in the system for the estimated displacement i+1{U} to not begin to diverge from the actual displacement. Specifically, Δt ≤ 2/ωmax, which is equivalent to Δt ≤ L/c where L is the shortest side of the smallest element and c is the velocity of sound (dilatational wave speed) in the element. For this reason, the smallest element side controls the time step size, and very small or highly distorted elements result in large solution times. For modeling the deformation around the cutting edge accurately, for instance, to accurately model residual stresses, subsurface delamination, etc., very small elements are required causing the time steps to be small and the simulation time to be large. Subcycling, whereby the nodal forces due to each of the elements is updated at different multiples of the critical time step of the smallest element is sometimes used, but is not very common. As an aside, an effective way to verify that a consistent set of units is being used for displacement, loads, stresses, etc., is to verify that the speed of dilatational waves (Vc = √ E/ρ), where E is the elastic modulus and ρ is the density) is obtained correctly from the values of E and ρ.

Rețineți că atunci când se utilizează integrarea explicită, chiar și o problemă pur cvasistatică trebuie rezolvată ca o problemă dinamică în care încărcările i{f} și accelerațiile prescrise sunt netede (nu au salturi) și forțele inerțiale sunt neglijabile. Dacă sarcinile sunt aplicate brusc, accelerațiile devin prea mari și conduc la rezultate false (mai ales atunci când se folosește scalarea masei, discutată în continuare). Într-un caz în care se rezolvă o problemă aproximativ cvasistatică, adică în care forța de inerție este mult mai mică decât forța de rigiditate (sau energia cinetică este mult mai mică decât energia internă), creșterea densității materialului, cunoscută sub numele de scalare a masei, este o strategie eficientă pentru a reduce viteza sunetului, a mări dimensiunea pasului de timp și a accelera simularea. Limita de scalare a masei ar trebui aleasă astfel încât modificarea energiei cinetice de la pas la pas (pe parcursul unei simulări) să fie întotdeauna o fracțiune neglijabilă din modificarea energiei interne pentru a se asigura că energia cinetică stocată are un efect neglijabil asupra deformației. În majoritatea pachetelor software, este mai ușor de urmărit energia cinetică totală și de asigurat că este o mică parte din energia totală internă. Cu toate acestea, pentru structurile mari cu regiuni mici de deformare activă (cum ar fi o zonă mică de forfecare primară (PSZ) la un capăt al unei așchii mari foarte deformate în care energia internă este mare), efectele inerțiale locale pot rămâne nedetectate dacă este utilizat criteriul energiei totală, mai degrabă decât criteriul de modificare incrementală.
Note that when using explicit integration, even a purely quasistatic problem has to be solved as a dynamic problem where loads i{f} and prescribed accelerations are smooth (do not have jumps) and inertial forces are negligible. If the loads are suddenly applied, the accelerations become too high and lead to spurious results (especially when mass scaling, discussed next, is used). In a case where an approximately quasistatic problem is being solved, i.e., where the inertial force is much smaller than the stiffness force (or the kinetic energy is much smaller than the internal energy), increasing the density of the material, known as mass scaling, is an effective strategy to decrease the speed of sound, increase the time step size, and speed up the simulation. The limit to mass scaling should be chosen such that the change in kinetic energy from step to step (throughout a simulation) is always a negligible fraction of the change in internal energy to ensure that stored kinetic energy has negligible effect on the deformation. In most software packages, it is easier to track the total kinetic energy and ensure that it is a small fraction of the total internal energy. However, for large structures with small regions of active deformation (such as a small primary shear zone (PSZ) at one end of a large highly deformed chip within which the internal energy is high), local inertial effects may go undetected if the total energy criterion, rather than the incremental change criterion, is used.

Rezolvând în mod repetat ecuația de mai sus pentru creșteri mici în timp, se pot obține deplasările în orice moment de timp. Observați că rezolvarea ecuației de mai sus este banală dacă [M] este o matrice diagonală; în caz contrar, această metodă este imposibilă din punct de vedere computațional. Din acest motiv, matricea de masă este „grupată” într-o matrice diagonală prin mutarea (însumarea) tuturor maselor din fiecare rând în termenul diagonal, reducând la zero ceilalți termeni. De fapt, masa și forța internă corespunzătoare fiecărui grad de libertate pot fi calculate cu ușurință, evitând necesitatea formulării acestor matrici. Astfel, problemele foarte mari pot fi rezolvate cu ușurință.
By repeatedly solving the above equation for small increments in time, the displacements at any instant of time can be obtained. Notice that solving the above equation is trivial if [M] is a diagonal matrix; otherwise, this method is computationally infeasible. For this reason, the mass matrix is “lumped” into a diagonal matrix by moving (summing) all the masses in each row into the diagonal term, zeroing out the other terms. In fact, mass and internal force corresponding to each degree of freedom can be readily calculated, obviating the need to formulate these matrices. Thus, very large problems can be handled easily.

În metodele implicite, deplasările i+1{U} în următorul moment de timp i+1t depind și de viteza i+1{Ú} și de accelerația i+1{Ü} la i+1t, în plus față de deplasarea, viteza și accelerația în momentele curente și anterioare. Acest lucru se întâmplă atunci când ecuațiile sistemului la momentul i+1t sunt folosite pentru a calcula deplasarea i+1{U}, așa cum se arată mai jos:
In implicit methods the displacements i+1{U} at the next instant of time i+1t also depend upon the velocity i+1{U} and acceleration i+1{U} at i+1t, in addition to the displacement, velocity, and acceleration at the current and previous instants. This occurs when the system equations at time i+1t are used to calculate displacement i+1{U}, as shown below

Rezolvarea ecuațiilor de mai sus pentru i+1{Ú} și i+1{Ü}
Solving the above equations for i+1 _U and i+1 €U

Folosind ecuația mișcării la i+1t,
Using the equation of motion at i+1t,

Ultimul termen dintre parantezele pătrate din dreapta este din nou accelerația i+1{Ű} iar ultimii doi termeni împreună estimează modificarea vitezei între timpul it și i+1t, datorită accelerației medii i+1/2{Ű}. Deoarece deplasarea în următorul moment de timp i+1t este obținută prin forțarea ecuației de mișcare în acel moment, metodele implicite sunt necondiționat stabile. Însă, este necesară iterația prin diferite valori de probă ale i+1{U} atunci când [M], [K] sau {f} depind de i+1{U}, așa cum se întâmplă în cazul structurilor cu răspuns, contact neliniar etc. De fapt, pentru multe modele constitutive, solicitarea nu poate fi direct legată de deplasare, ci doar de solicitarea curentă și de creșterea deformării. Creșterea forței interne este proporțională cu deplasarea incrementală printr-o matrice de rigiditate tangentă. O formă a matricei de rigiditate tangentă este calculată ca ∰ve =[ B]^T[C][B]dV, unde [C] raportează creșterea solicitării cu creșterea deformării ca [σ́]= C[ε ́] (Belytschko et al. 2000, Sect. 6.4). Procesul de actualizare a matricei solicitărilor și rigidității tangente necesită de obicei iterare folosind scheme predictor-corector pentru a se asigura că constrângerile impuse solicitărilor și incrementelor de deformare de către modelul constitutiv sunt satisfăcute.
The last term within the square brackets on the right is again the acceleration i+1{Ű} and the last two terms together estimate the change in velocity between time it and i+1t, due to the average acceleration i+1/2{Ű}. Since the displacement at the next instant of time i+1t is obtained by enforcing the equation of motion at that time, implicit methods are unconditionally stable. However, iteration through various trial values of i+1{U} is needed when [M], [K], or {f} depend on i+1{U}, as occurs for structures with nonlinear response, contact, etc. In fact, for many constitutive models, the stress cannot be directly related to the displacement, but only to the current stress and strain increment. The increment of internal force is proportional to the incremental displacement via a tangent stiffness matrix. One form of the tangent stiffness matrix is calculated as ∰ve =[B]^T[C][B]dV, where [C] relates the increment of stress to increment of strain as [σ́]=C[έ] (Belytschko et al. 2000, Sect. 6.4). The process of updating the stress and tangent stiffness matrix typically requires iteration using predictor–corrector schemes to ensure that the constraints imposed on the stresses and strain increments by the constitutive model are satisfied.

Plasticitatea, contactul, efectele temperaturii, evoluția geometriei etc., fac din prelucrare o problemă extrem de neliniară. Dar, dimensiunea pasului de timp poate fi de 20-100 de ori dimensiunea pasului de timp critic în metodele explicite, limitată doar de necesitatea de a urmări evoluția [M], [K] și {f}. Totuși, de observat că inversarea matricei este necesară pentru rezolvarea ec. 11, deoarece [K] are termeni ne-diagonali. Acest lucru face ca metodele implicite să fie sensibile la dimensiunea problemei, iar timpul de rezolvare crește ca a treia putere a numărului de grade de libertate. În timp ce rafinarea locală a ochiurilor, condensarea regiunilor liniare ale structurii în superelemente etc., sunt strategii eficiente pentru a menține dimensiunea problemei gestionabilă, dacă dimensiunea pasului de timp scade, de asemenea, odată cu rafinarea locală a ochiurilor (de exemplu, cu deformare de contact sau localizată), s timpul de soluție poate crește și mai mult. În plus, dificultățile de convergență în prezența gradienților foarte mari în deformare și a interacțiunilor neliniare foarte cuplate (cum ar fi în benzi de forfecare adiabatică) pot determina eșecul soluției. Unele pachete software de uz general sunt predispuse la acest lucru, iar depășirea dificultăților de convergență este uneori o experiență frustrantă.
Plasticity, contact, temperature effects, geometry evolution, etc., make machining a highly nonlinear problem. Still, the time step size can be 20–100 times the critical time step size in explicit methods, limited only by the need to track the evolution of [M], [K], and {f}. However, notice that matrix inversion is required for solving Eq. 11, because [K] has non-diagonal terms. This makes implicit methods sensitive to problem size, and solution time increases as the third power of the number of degrees of freedom. While local mesh refinement, condensation of linear regions of the structure into superelements, etc., are effective strategies to keep the problem size manageable, if the time step size also decreases with local mesh refinement (for instance with contact or localized deformation) the solution time can increase even more. Moreover, convergence difficulties in the presence of very high gradients in strain and highly coupled nonlinear interactions (such as in adiabatic shear banding) can cause the solution to fail. Some general purpose software packages are prone to this, and overcoming convergence difficulties is sometimes a frustrating experience.