CHAP 2 Nonlinear Finite Element Analysis Procedures - PowerPoint Presentation

1. CHAP 2Nonlinear Finite Element Analysis ProceduresNam-Ho Kim1

2. GoalsWhat is a nonlinear problem?How is a nonlinear problem different from a linear one?What types of nonlinearity exist?How to understand stresses and strainsHow to formulate nonlinear problemsHow to solve nonlinear problemsWhen does nonlinear analysis experience difficulty?

3. Nonlinear Structural ProblemsWhat is a nonlinear structural problem?Everything except for linear structural problemsNeed to understand linear problems firstWhat is linearity?Example: fatigue analysisInput x(load, heat)Output y(displ, temp)x1y1x2y22x12y12x1 +3x22y1+3y2x1x2Y

4. What is a linear structural problem?Linearity is an approximationAssumptions: Infinitesimal strain (<0.2%)Infinitesimal displacementSmall rotationLinear stress-strain relationF/2F/2A0AL0dLLEses = EeForceStressStrainDisplacementLinearLinearLinearGlobalGlobalLocalLocal

5. Observations in linear problemsWhich one will happen?Will this happen?MMFTrussTruss

6. What types of nonlinearity exist?It is at every stage of analysis

7. Linear vs. Nonlinear ProblemsLinear Problem:Infinitesimal deformation: Linear stress-strain relation:Constant displacement BCsConstant applied forcesNonlinear Problem:Everything except for linear problems!Geometric nonlinearity: nonlinear strain-displacement relationMaterial nonlinearity: nonlinear constitutive relationKinematic nonlinearity: Non-constant displacement BCs, contactForce nonlinearity: follow-up loadsUndeformed coord.Constant

8. Nonlinearities in Structural ProblemsMore than one nonlinearity can exist at the same timeDisplacementStrainStressPrescribeddisplacementApplied forceNonlinear displacement-strainNonlinear stress-strainNonlinear displ. BCNonlinear force BC

9. Geometric NonlinearityRelations among kinematic quantities (i.e., displacement, rotation and strains) are nonlinearDisplacement-strain relationLinear:Nonlinear: 0.0 0.2 0.4 0.6 0.8 1.0Normalized couple8.6.4.2.0.Tip displacementC0C1C2C3When du/dx is smallH.O.T. can be ignored

10. Geometric Nonlinearity cont.Displacement-strain relationE has a higher-order term(du/dx) << 1  e(x) ~ E(x).Domain of integrationUndeformed domain W0Deformed domain Wx0 0.05 0.10 0.15 0.20 0.25 0.30du/dx0.350.300.250.200.150.100.050StraineEDeformed domain is unknown

11. Material NonlinearityLinear (elastic) materialOnly for infinitesimal deformationNonlinear (elastic) material[D] is not a constant but depends on deformationStress by differentiating strain energy density ULinear material:Stress is a function of strain (deformation): potential, path independent1esEesELinear springNonlinear springLinear and nonlinear elastic spring modelsMore generally, {s} = {f(e)}

12. Material Nonlinearity cont.Elasto-plastic material (energy dissipation occurs)Friction plate only support stress up to syStress cannot be determined from stress aloneHistory of loading path is required: path-dependentVisco-elastic materialTime-dependent behaviorCreep, relaxationVisco-elastic spring modelsEhtimesEeElasto-plastic spring modelsEsYesEsY

13. Boundary and Force NonlinearitiesNonlinear displacement BC (kinematic nonlinearity)Contact problems, displacement dependent conditionsNonlinear force BC (Kinetic nonlinearity)Contact boundarydmaxDisplacementForce

14. Mild vs. Rough NonlinearityMild Nonlinear ProblemsContinuous, history-independent nonlinear relations between stress and strainNonlinear elasticity, Geometric nonlinearity, and deformation-dependent loadsRough Nonlinear ProblemsEquality and/or inequality constraints in constitutive relationsHistory-dependent nonlinear relations between stress and strainElastoplasticity and contact problems

15. Nonlinear Finite Element EquationsEquilibrium between internal and external forcesKinetic and kinematic nonlinearitiesAppears on the boundaryHandled by displacements and forces (global, explicit)Relatively easy to understand (Not easy to implement though)Material & geometric nonlinearitiesAppears in the domainDepends on stresses and strains (local, implicit)Linear problemsStressStrainLoads

16. Solution ProcedureWe can only solve for linear problems …

17. Example – Nonlinear SpringsSpring constantsk1 = 50 + 500u k2 = 100 + 200u Governing equationSolution is in the intersection between two zero contoursMultiple solutions may existNo solution exists in a certain situationk1k2u1u2FP1P2

18. Solution ProcedureLinear ProblemsStiffness matrix K is constantIf the load is doubled, displacement is doubled, tooSuperposition is possibleNonlinear ProblemsHow to find d for a given F?FddiKT2F2ddKFIncremental Solution Procedure

19. Newton-Raphson MethodMost popular methodAssume di at i-th iteration is knownLooking for di+1 from first-order Taylor series expansion : Jacobian matrix or Tangent stiffness matrixSolve for incremental solutionUpdate solutiondidi+1dDdi+1FP(d)di+2dnSolutionP(di)P(di+1)

20. N-R Method cont.Observations:Second-order convergence near the solution (Fastest method!)Tangent stiffness is not constantThe matrix equation solves for incremental displacementRHS is not a force but a residual forceIteration stops when conv < tolerance Or,

21. N-R AlgorithmSet tolerance = 0.001, k = 0, max_iter = 20, and initial estimate u = u0Calculate residual R = f – P(u)Calculate conv. If conv < tolerance, stop If k > max_iter, stop with error messageCalculate Jacobian matrix KTIf the determinant of KT is zero, stop with error messageCalculate solution increment uUpdate solution by u = u + uSet k = k + 1Go to Step 2

22. Example – N-R MethodIteration 1

23. Example – N-R Method cont.Iteration 2Iteration 3

24. Example – N-R Method cont.Iteration 4IterationResidualQuadratic convergenceIter||R||017.26314.53120.01630.0

25. When N-R Method Does Not ConvergeDifficultiesConvergence is not always guaranteedAutomatic load step control and/or line search techniques are often usedDifficult/expensive to calculatedidi+1dFP(d)di+2dnSolution

26. When N-R Method Does Not Converge cont.Convergence difficulty occurs whenJacobian matrix is not positive-definiteBifurcation & snap-through require a special algorithmABCDEABCDEDisplacementForceFFBFCP.D. Jacobian: in order to increase displ., force must be increased

27. Modified N-R MethodConstructing and solving is expensiveComputational Costs (Let the matrix size be N x N)L-U factorization ~ N3Forward/backward substitution ~ NUse L-U factorized repeatedlyMore iteration is required, buteach iteration is fastMore stable than N-R methodHybrid N-R methoddidi+1dFP(d)dnSolution

28. Example – Modified N-R MethodSolve the same problem using modified N-R methodIteration 1

29. Example – Modified N-R Method cont.Iteration 2IterationResidualIter||R||017.26314.531020.358430.083140.020450.005160.001370.0003

30. Incremental Secant MethodSecant matrixInstead of using tangent stiffness, approximate it using the solution from the previous iterationAt i-th iterationThe secant matrix satisfiesNot a unique process in high dimensionStart from initial KT matrix, iteratively update itRank-1 or rank-2 update The textbook has Broyden’s algorithm (Rank-1 update)Here we will discuss BFGS method (Rank-2 update)Fd0d1dP(d)dnSolutiond2d3Secant stiffness

31. Incremental Secant Method cont.BFGS (Broyden, Fletcher, Goldfarb and Shanno) methodStiffness matrix must be symmetric and positive-definiteInstead of updating K, update H (saving computational time)Become unstable when the No. of iterations is increased

32. Incremental Force MethodN-R method converges fast if the initial estimate is close to the solutionSolid mechanics: initial estimate = undeformed shape Convergence difficulty occurs when the applied load is large (deformation is large)IFM: apply loads in increments. Use the solution from the previous increment as an initial estimateCommercial programs call it “Load Increment” or “Time Increment”

33. Incremental Force Method cont.Load increment does not have to be uniformCritical part has smaller increment sizeSolutions in the intermediate load incrementsHistory of the response can provide insight into the problemEstimating the bifurcation point or the critical loadLoad increments greatly affect the accuracy in path-dependent problems

34. Load Increment in Commercial SoftwareUse “Time” to represent load levelIn a static problem, “Time” means a pseudo-timeRequired Starting time, (Tstart), Ending time (Tend) and incrementLoad is gradually increased from zero at Tstart and full load at TendLoad magnitude at load increment Tn:Automatic time steppingIncrease/decrease next load increment based on the number of convergence iteration at the current loadUser provide initial load increment, minimum increment, and maximum incrementBisection of load increment when not converged

35. Force Control vs. Displacement ControlForce control: gradually increase applied forces and find equilibrium configurationDispl. control: gradually increase prescribed displacementsApplied load can be calculated as a reactionMore stable than force control.Useful for softening, contact, snap-through, etc.u1 u2 u F P(u) un F2 F3 Fn u3 F1 uA uB u F P(u) uD FB FC uC FA 

36. Nonlinear Solution StepsInitialization:Residual CalculationConvergence Check (If converged, stop)LinearizationCalculate tangent stiffness Incremental Solution:SolveState DeterminationUpdate displacement and stressGo To Step 2

37. Nonlinear Solution Steps cont.State determinationFor a given displ dk, determine current state (strain, stress, etc)Sometimes, stress cannot be determined using strain aloneResidual calculationApplied nodal force − Nodal forces due to internal stressesWeak form:Discretization:Residual:

38. Example – Linear Elastic MaterialGoverning equation (Scalar equation)Collect ResidualLinear elastic materialdFKT

39. Example – Nonlinear BarRubber barDiscrete weak formScalar equation0 0.01 0.02 0.03 0.04 0.05Strain120100806040200StressF = 10kNL = 1m12x

40. Example – Nonlinear Bar cont.JacobianN-R equationIteration 1Iteration 2

41. N-R or Modified N-R?It is always recommended to use the Incremental Force MethodMild nonlinear: ~10 incrementsRough nonlinear: 20 ~ 100 incrementsFor rough nonlinear problems, analysis results depends on increment sizeWithin an increment, N-R or modified N-R can be usedN-R method calculates KT at every iterationModified N-R method calculates KT once at every incrementN-R is better when: mild nonlinear problem, tight convergence criterionModified N-R is better when: computation is expensive, small increment size, and when N-R does not converge wellMany FE programs provide automatic stiffness update optionDepending on convergence criteria used, material status change, etc

42. Accuracy vs. ConvergenceNonlinear solution procedure requiresInternal force P(d)Tangent stiffnessThey are often implemented in the same routineInternal force P(d) needs to be accurateWe solve equilibrium of P(d) = FTangent stiffness KT(d) contributes to convergenceAccurate KT(d) provides quadratic convergence near the solutionApproximate KT(d) requires more iteration to convergeWrong KT(d) causes lack of convergence

43. Convergence CriteriaMost analysis programs provide three convergence criteriaWork, displacement, load (residual)Work = displacement * loadAt least two criteria needs to be convergedTraditional convergence criterion is load (residual)Equilibrium between internal and external forcesUse displacement criterion for load insensitive systemForceDisplacementUse loadcriterionUse displacementcriterion

44. Load Step (subcase or step)Load step is a set of loading and boundary conditions to define an analysis problemMultiple load steps can be used to define a sequence of loading conditionsSolution StrategiesLS1LS2LoadTimeNASTRANSPC = 1SUBCASE 1 LOAD = 1SUBCASE 2 LOAD = 2

45. Solution StrategiesLoad Increment (substeps)Linear analysis concerns max load Nonlinear analysis depends on load path (history)Applied load is gradually increasedwithin a load stepFollow load path, improve accuracy, and easy to convergeConvergence IterationWithin a load increment, an iterative method (e.g., NR method) is used to find nonlinear solutionBisection, linear search, stabilization, etcFaFu1234500.20.40.60.81-100-50050100150200250300DisplacementForceLoadingUnloading

46. Solution Strategies cont.Automatic (Variable) Load IncrementAlso called Automatic Time SteppingLoad increment may not be uniformWhen convergence iteration diverges, the load increment is halvedIf a solution converges in less than 4 iterations, increase time increment by 25%If a solution converges in more than 8 iterations, decrease time increment by 25%Subincrement (or bisection)When iterations do not converge at a given increment, analysis goes back to previously converged increment and the load increment is reduced by halfThis process is repeated until max number of subincrements

47. When nonlinear analysis does not convergeNR method assumes a constant curvature locallyWhen a sign of curvature changes around the solution, NR method oscillates or divergesOften the residual changes sign between iterationsLine search can help to converge

48. When nonlinear analysis does not convergeDisplacement-controlled vs. force-controlled procedureAlmost all linear problems are force-controlledDisplacement-controlled procedure is more stable for nonlinear analysisUse reaction forces to calculate applied forcesABCDEFABCDEDisplacementForceFBFC

49. When nonlinear analysis does not convergeMesh distortionMost FE programs stop analysis when mesh is distorted too muchInitial good mesh may be distorted during a large deformationMany FE programs provide remeshing capability, but it is still inaccurate or inconvenientIt is best to make mesh in such a way that the mesh quality can be maintained after deformation (need experience)Initial mesh

50. MATLAB Code for Nonlinear FEA

51. NLFEA.mNonlinear finite element analysis programIncremental force method with N-R methodBisection method when N-R is failed to convergeCan solve for linear elastic, hyperelastic and elasto-plastic material nonlinearities with large deformationGlobal arraysNameDimensionContentsGKFNEQ x NEQTangent matrixFORCENEQ x 1Residual vectorDISPTDNEQ x 1Displacement vectorDISPDDNEQ x 1Displacement incrementSIGMA6 x 8 x NEStress at each integration pointXQ7 x 8 x NEHistory variable at each integration point

52. StopUpdate history variablesPrint stress & displacementT = T + DTBisection controlFinal time?Input dataIncrease load & BCITER = 0Calculate R & KDisplacement BCConverged?Max ITER?Solve DU = K\RU = U + DUITER=ITER + 1T = T – DTDT = DT/2YesNoYesNoYesNo

53. NLFEA.m cont.Nodal coordinates and element connectivitythe node numbers are in sequencenodal coordinates in XYZ(NNODE , 3)eight-node hexahedral elements LE(NELEN, 8)Applied forces and prescribed displacementsEXTFORCE(NFORCE, 3): [node, DOF, value] formatSDISPT(NDISPT, 3)Load steps and incrementsTIMS(NTIME,5): [Tstart, Tend, Tinc, LOADinit, LOADfinal] formatMaterial propertiesMooney-Rivlin hyperelasticity (MID = -1), PROP = [A10, A01, K]infinitesimal elastoplasticity (MID = 1), PROP = [LAMBDA, MU, BETA, H, Y0]

54. NLFEA.m cont.Control parametersITRA: maximum number of convergence iterationsif residual > ATOL, then solution diverges, bisection startsThe total number of bisections is limited by NTOLThe convergence iteration converges when residual < TOLProgram prints out results to NOUT after convergencefunction NLFEA(ITRA,TOL,ATOL,NTOL,TIMS,NOUT,MID,PROP,EXTFORCE,SDISPT,XYZ,LE)%***********************************************************************% MAIN PROGRAM FOR HYPERELASTIC/ELASTOPLASTIC ANALYSIS%***********************************************************************

55. %% Nodal coordinatesXYZ=[0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1];%% Element connectivityLE=[1 2 3 4 5 6 7 8];%% External forces [Node, DOF, Value]EXTFORCE=[5 3 10.0; 6 3 10.0; 7 3 10.0; 8 3 10.0];%% Prescribed displacements [Node, DOF, Value]SDISPT=[1 1 0;1 2 0;1 3 0;2 2 0;2 3 0;3 3 0;4 1 0;4 3 0];%% Load increments [Start End Increment InitialLoad FinalLoad]TIMS=[0.0 0.5 0.1 0.0 0.5; 0.5 1.0 0.1 0.5 1.0]';%% Material properties%PROP=[LAMBDA MU BETA H Y0]MID=1;PROP=[110.747, 80.1938, 0.0, 5., 35.0];%% Set program parametersITRA=20; ATOL=1.0E5; NTOL=5; TOL=1E-6;%% Calling main functionNOUT = fopen('output.txt','w');NLFEA(ITRA, TOL, ATOL, NTOL, TIMS, NOUT, MID, PROP, EXTFORCE, SDISPT, XYZ, LE);fclose(NOUT);Extension of a Single Element Example

56. Tension of Elastoplastic Bar ExamplelmsYH110.7 GPa80.2 GPa400 MPa100 MPa11kNx21568742x13x391011kN11kN11kN1211Material propertiesUniaxial stress condition s33 ≠ 0When Fi = 10 kN, s = 400 MPa (Elastic limit)Elastoplastic when Fi = 10 ~ 11kN%% Two-element example%% Nodal coordinatesXYZ=[0 0 0; 1 0 0; 1 1 0; 0 1 0; 0 0 1; 1 0 1; 1 1 1; 0 1 1; 0 0 2; 1 0 2; 1 1 2; 0 1 2]*0.01;%% Element connectivityLE=[1 2 3 4 5 6 7 8; 5 6 7 8 9 10 11 12];%% Prescribed displacements [Node, DOF, Value]SDISPT=[1 1 0;1 2 0;1 3 0;2 2 0;2 3 0;3 3 0;4 1 0;4 3 0];1432x1x2u=v=w=0v=w=0w=0u=w=0

57. Tension of Elastoplastic Bar Example%% External forces [Node, DOF, Value]EXTFORCE=[9 3 10.0E3; 10 3 10.0E3; 11 3 10.0E3; 12 3 10.0E3];%% Load increments [Start End Increment InitialFactor FinalFactor]TIMS=[0.0 0.8 0.4 0.0 0.8; 0.8 1.1 0.1 0.8 1.1]';%% Material properties PROP=[LAMDA MU BETA H Y0]MID=1;PROP=[110.747E9 80.1938E9 0.0 1.E8 4.0E8];%% Set program parametersITRA=70; ATOL=1.0E5; NTOL=6; TOL=1E-6;%% Calling main functionNOUT = fopen('output.txt','w');NLFEA(ITRA, TOL, ATOL, NTOL, TIMS, NOUT, MID, PROP, EXTFORCE, SDISPT, XYZ, LE);fclose(NOUT);01.10.40.90.814891011TimeForce (kN)10kN * 1.1 = 11kN

58. Tension of Elastoplastic Bar Example Time Time step Iter Residual 0.40000 4.000e-01 2 3.80851e-12   Time Time step Iter Residual 0.80000 4.000e-01 2 4.32010e-12   Time Time step Iter Residual 0.90000 1.000e-01 2 3.97904e-12   Time Time step Iter Residual 1.00000 1.000e-01 2 3.63798e-12   Time Time step Iter Residual 1.10000 1.000e-01 2 6.66390e+02 3 1.67060e-09Convergence iteration outputs (output.txt)Linear elastic regionElastoplastic regionLoad factoru5zu9zS33 Elem1S33 Elem2State0.47.73×10−61.55×10−5160 MPa160 MPaElastic0.81.55×10−53.09×10−5320 MPa320 MPaElastic0.91.74×10−53.48×10−5360 MPa360 MPaElastic1.01.93×10−53.87×10−5400 MPa400 MPaElastic1.14.02×10−38.04×10−3440 MPa440 MPaPlastic