CHAP 2 Nonlinear Finite Element Analysis Procedures

CHAP 2Nonlinear Finite Element AnalysisProceduresNam-Ho Kim Goals What is a nonlinear problem? How is a nonlinear problem different from a linear one? What types of nonlinearity exist? How to understand stresses and strains How to formulate nonlinear problems How to solve nonlinear problems When does nonlinear analysis experience difficulty?

Nonlinear Structural Problems What is a nonlinear structural problem?– Everything except for linear structural problems– Need to understand linear problems firstA(Du Ew) What is linearity?yInput x(load, heat)x1axOutput y(displ, temp)y1x2y2x1x1y2axDA(u) EA(w)x22y12x1 3x2Y2y1 3y2 Example: fatigue analysis What is a linear structural problem?ForceLinearLinearStressA0 AGlobalStrainLEFF?A0 A(F)HHA0EGLL0 Linearity is an approximation earV EHVF/2VLocalLocalGlobalA0EHInfinitesimal strain ( 0.2%)Infinitesimal displacementSmall rotationLinear stress-strain relationGL GL?L0 L

Observations in linear problems Which one will happen?MM Will this happen?TrussTrussF What types of nonlinearityexist?It is at every stage of analysis

Linear vs. Nonlinear Problems Linear Problem:1 § wui wujH – Infinitesimal deformation: ij2 wxj wxiD: – Linear stress-strain relation: – Constant displacement BCs· ¹8QGHIRUPHG FRRUG &RQVWDQW– Constant applied forces Nonlinear Problem:– Everything except for linear problems!– Geometric nonlinearity: nonlinear strain-displacement relation– Material nonlinearity: nonlinear constitutive relation– Kinematic nonlinearity: Non-constant displacement BCs, contact– Force nonlinearity: follow-up loads Nonlinearities in Structural nlinear displ. s-strainNonlinear force BCAppliedforce More than one nonlinearity can exist at the same time

Geometric Nonlinearity ϭdŝƉ ĚŝƐƉůĂĐĞŵĞŶƚ Relations among kinematic quantities (i.e., displacement,rotation and strains) are nonlinear Ϯϴ͘ϲ͘ϰ͘Ϯ͘ ϯ Ϭ͘Ϭ͘Ϭ Ϭ͘Ϯ Ϭ͘ϰ Ϭ͘ϲ Ϭ͘ϴ ϭ͘ϬEŽƌŵĂůŝnjĞĚ ĐŽƵƉůĞWhen du/dx is small Displacement-strain relation– Linear:2H(x)dudxE(x)du 1 § du · dx 2 dx ¹du§ du · dx dx ¹2– Nonlinear:H.O.T. can be ignoredH(x) E(x)Geometric Nonlinearity cont. Displacement-strain relationϬ͘ϯϱϬ͘ϯϬ– E has a higher-order termϬ͘Ϯϱ– (du/dx) 1 Æ H(x) E(x).F ƚƌĂŝŶϬ͘ϮϬϬ͘ϭϱϬ͘ϭϬ Domain of integration– Undeformed domain :0Ϭ͘ϬϱϬϬ Ϭ͘Ϭϱ – Deformed domain :xa( u, u )Ϭ͘ϭϬ Ϭ͘ϭϱ ĚƵͬĚdžϬ͘ϮϬ Ϭ͘Ϯϱ Ϭ͘ϯϬ³³: H(u ) : V(u) d:Deformed domain is unknown

Material Nonlinearity Linear (elastic) material{ V } [D]{ H }– Only for infinitesimal deformation Nonlinear (elastic) materialMore generally, {V} {f(H)}– [D] is not a constant but depends on deformation– Stress by differentiating strain energy density U– Linear material:dU1U2EH2VdHVdUdHEH– Stress is a function of strain (deformation): potential, pathindependentVEVHϭNonlinearspringE Linear springHLinear and nonlinear elastic spring models Material Nonlinearity cont. Elasto-plastic material (energy dissipation occurs)– Friction plate only support stress up to Vy– Stress cannot be determined from stress alone– History of loading path is required: path-dependentEVTYVVYEHElasto-plastic spring model Visco-elastic material– Time-dependent behavior– Creep, relaxationEVKVEHtimeVisco-elastic spring model

Boundary and Force Nonlinearities Nonlinear displacement BC (kinematic nonlinearity)– Contact problems, displacement dependent conditions&ŽƌĐĞEŵĂdž ŽŶƚĂĐƚ ďŽƵŶĚĂƌLJ ŝƐƉůĂĐĞŵĞŶƚ Nonlinear force BC (Kinetic nonlinearity) Mild vs. Rough Nonlinearity Mild Nonlinear Problems– Continuous, history-independent nonlinear relations betweenstress and strain– Nonlinear elasticity, Geometric nonlinearity, and deformationdependent loads Rough Nonlinear Problems– Equality and/or inequality constraints in constitutive relations– History-dependent nonlinear relations between stress and strain– Elastoplasticity and contact problems

Nonlinear Finite Element Equations Equilibrium between internal and external forcesP(d)StressStrainF(d)Linear problems[K]{ d}{F }Loads Kinetic and kinematic nonlinearities– Appears on the boundary– Handled by displacements and forces (global, explicit)– Relatively easy to understand (Not easy to implement though) Material & geometric nonlinearities– Appears in the domain– Depends on stresses and strains (local, implicit) Solution ProcedureWe can only solve for linear problems

Example – Nonlinear SpringsŬϭ Spring constantsŬϮ– k1 50 500uk2 100 200uƵϭƵϮ& Governing equation 300u12 400u1u2 200u22 150u1 100u2 ̄ 200u12 400u1u2 200u22 100u1 100u20P1100P2– Solution is in the intersection between two zero contours– Multiple solutions may exist– No solution exists in a certain situation Solution Procedure Linear ProblemsK dForP(d)F– Stiffness matrix K is constantP(d1 d2 ) P(d1 ) P(d2 )P(Dd) DP(d) DF FKFd d– If the load is doubled, displacement is doubled, too– Superposition is possible Nonlinear ProblemsP(d)F,P(2d) z 2FFKT– How to find d for a given F?Incremental Solution Proceduredid

Newton-Raphson Method Most popular method Assume di at i-th iteration is known Looking for di 1 from first-order Taylor series expansionP(di 1 ) P(di ) KTi (di ) 'diFi–KTi (di )§ wP ·{ wd ¹: Jacobian matrix or Tangent stiffness matrix Solve for incremental solutionKTi 'diF P(di )& Update solutiondi 1W;ĚŝнϭͿdi 'diW;ĚͿKTi �нϮĚŶĚ N-R Method cont. Observations:– Second-order convergence near the solution (Fastest method!)– Tangent stiffnessKTi (di )limis not constantn ofuexact un 1uexact un2c– The matrix equation solves for incremental displacement 'di– RHS is not a force but a residual force Ri { F P(di )– Iteration stops when conv tolerance j 1 (Rji 1 )2n1 j 1 (Fj )2nconv j 1 ('uij 1 )2n1 j 1 ( 'uj0 )2nOr,conv

N-R Algorithm1. Set tolerance 0.001, k 0, max iter 20, and initialestimate u u02. Calculate residual R f – P(u)3. Calculate conv. If conv tolerance, stop4. If k max iter, stop with error message5. Calculate Jacobian matrix KT6. If the determinant of KT is zero, stop with error message7. Calculate solution increment 'u8. Update solution by u u 'u9. Set k k 110. Go to Step 2 Example – N-R Method d d2 ½P(d) { 212¾ d1 d2 ¿KTwPwd 3 ½ ¾{F 9 ¿1 ºª 1« 2d 2d » 12¼d0 1 ½ ¾ 5 ¿R0F P( d 0 )P( d 0 ) 6 ½ ¾ 26 ¿ 3 ½ ¾ 17 ¿ Iteration 1ª 1 1 º 'd10 ½ ¾«» 2 10 ¼ 'd20 ¿ 3 ½ ¾ 17 ¿d1d0 'd0 0.625 ½ ¾ 3.625 ¿R1F P(d1 ) 0 ½ ¾ 4.531 ¿0 'd1 ½ 0¾ 'd2 ¿ 1.625 ½ ¾ 1.375 ¿

Example – N-R Method cont. Iteration 21 º 'd11 ½ ª 1« 1.25 7.25 » 1 ¾ ¼ 'd2 ¿d2d1 'd11 'd1 ½ 1¾ 'd2 ¿ 0 ½ ¾ 4.531 ¿ 0.092 ½ ¾ 3.092 ¿R2F P(d2 ) 0.533 ½ ¾ 0.533 ¿ 0 ½ ¾ 0.568 ¿ Iteration 31 º 'd12 ½ª 1« 0.184 6.184 » 2 ¾ ¼ 'd2 ¿d3d2 'd22 'd1 ½ 2¾ 'd2 ¿ 0 ½ ¾ 0.568 ¿ 0.003 ½ ¾ 3.003 ¿R3F P(d3 ) 0.089 ½ ¾ 0.089 ¿ 0 ½ ¾ 0.016 ¿ Example – N-R Method cont. Iteration 411 º 'd13 ½ ª« 0.005 6.005 » 3 ¾ ¼ 'd2 ¿d4d3 'd3 'd13 ½ 3¾ 'd2 ¿ 0 ½ ¾ 0.016 ¿ 0.000 ½ ¾ 3.000 ¿R4F P( d 4 ) 0.003 ½ ¾ 0.003 ¿ 0 ½ ¾ 0 ¿ResidualϮϬ,WHU5 ϴ ϰ Ϭ ϭϲϭϮϬϭϮIterationϯϰQuadratic convergence

When N-R Method Does Not Converge Difficulties– Convergence is not always guaranteed– Automatic load step control and/or line search techniques areoften used– Difficult/expensive to calculate KTi (di )&wPwd ŽůƵƚŝŽŶW;ĚͿwPwdĚŝнϮ ĚŝĚŶĚŝнϭĚ When N-R Method Does Not Converge cont. Convergence difficulty occurs when– Jacobian matrix is not positive-definiteP.D. Jacobian: in order to increase displ., force must be increased– Bifurcation & snap-through require a special algorithm & & &ŽƌĐĞ& ŝƐƉůĂĐĞŵĞŶƚ

Modified N-R Method Constructing KTi (di ) and solving KTi 'diRi is expensive Computational Costs (Let the matrix size be N x N)– L-U factorization N3– Forward/backward substitution N Use L-U factorized KTi (di ) repeatedly More iteration is required, butW;ĚͿ&each iteration is fastwPwd More stable than N-R method ŽůƵƚŝŽŶ Hybrid N-R methodĚŝĚŝнϭĚŶĚ Example – Modified N-R Method Solve the same problem using modified N-R method d d2 ½P(d) { 212¾ d1 d2 ¿KTwPwd 3 ½ ¾{F 9 ¿1 ºª 1« 2d 2d » 12¼d0R0 1 ½ ¾ 5 ¿P( d 0 )F P( d 0 ) 6 ½ ¾ 26 ¿ 3 ½ ¾ 17 ¿ Iteration 1ª 1 1 º 'd10 ½ ¾«» 2 10 ¼ 'd20 ¿ d1d0 'd0 'd10 ½ 0¾ 'd2 ¿ 3 ½ ¾ 17 ¿ 0.625 ½ ¾ 3.625 ¿R1 1.625 ½ ¾ 1.375 ¿F P(d1 ) 0 ½ ¾ 4.531 ¿

Example – Modified N-R Method cont. Iteration 2ª 1 1 º 'd11 ½ « 2 10 » 1 ¾ ¼ 'd2 ¿ 0.059 ½ ¾ 3.059 ¿d1 'd1d2 'd11 ½ 1¾ 'd2 ¿ 0 ½ ¾ 4.531 ¿ 0.566 ½ ¾ 0.566 ¿F P(d2 )R2Residual 0 ½ ¾ 0.358 ¿,WHU5 ϭϲ ϭϮ ϴ ϰ ϮϬϬϬϭϮϯϰϱϲϳIteration Incremental Secant Method Secant matrix– Instead of using tangent stiffness, approximate it using thesolution from the previous iteration– At i-th iterationKsi 'diF P(di )wPwd– The secant matrix satisfiesKsii 1 (d d )iW;ĚͿ& ŽůƵƚŝŽŶ ĞĐĂŶƚ ƐƚŝĨĨŶĞƐƐi 1P(d ) P(d )i– Not a unique process in high dimensionĚ ĚϭĚϮ ĚϯĚŶĚ Start from initial KT matrix, iteratively update it– Rank-1 or rank-2 update– The textbook has Broyden’s algorithm (Rank-1 update)– Here we will discuss BFGS method (Rank-2 update)

Incremental Secant Method cont. BFGS (Broyden, Fletcher, Goldfarb and Shanno) method– Stiffness matrix must be symmetric and positive-definite'di[Ksi ] 1 {F P(di )} { [Hsi ]{F P(di )}– Instead of updating K, update H (saving computational time)Hsi( I wi viT )Hsi 1 ( I wi viT )vi§( 'di 1 )T (Ri 1 Ri ) ·iRi 1 1 Ri T i 1( 'd ) R ¹wi'di 1( 'di 1 )T (Ri 1 Ri ) Become unstable when the No. of iterations is increased Incremental Force Method N-R method converges fast if the initial estimate is closeto the solution Solid mechanics: initial estimate undeformed shape Convergence difficultyoccurs when the appliedload is large(deformation is large) IFM: apply loads inincrements. Use thesolution from theprevious incrementas an initial estimate Commercial programscall it “Load Increment”or “Time Increment”

Incremental Force Method cont. Load increment does not have to be uniform– Critical part has smaller increment size Solutions in the intermediate load increments– History of the response can provide insight into the problem– Estimating the bifurcation point or the critical load– Load increments greatly affect the accuracy in path-dependentproblems Load Increment in Commercial Software Use “Time” to represent load level– In a static problem, “Time” means a pseudo-time– Required Starting time, (Tstart), Ending time (Tend) and increment– Load is gradually increased from zero at Tstart and full load at Tend– Load magnitude at load increment Tn:nFT n TstartFTend TstartTnn u 'T d Tend Automatic time stepping– Increase/decrease next load increment based on the number ofconvergence iteration at the current load– User provide initial load increment, minimum increment, andmaximum increment– Bisection of load increment when not converged

Force Control vs. Displacement Control Force control: gradually increase applied forces and findequilibrium configuration Displ. control: gradually increase prescribed displacements– Applied load can be calculated as a reaction– More stable than force control.– Useful for softening, contact, snap-through, etc.&&W;ƵͿ&Ŷ& &ϯW;ƵͿ& &Ϯ&ϭ& Ƶϭ ƵϮƵϯƵƵŶƵ Ƶ Ƶ Ƶ Ƶ Nonlinear Solution Steps1. Initialization: d00; i2. Residual CalculationRi0F P(di )3. Convergence Check (If converged, stop)4. Linearization–Calculate tangent stiffness KTi (di )5. Incremental Solution:–Solve KTi (di )'diRi6. State Determination–Update displacement and stress7. Go To Step 2di 1di 'diVi 1Vi 'Vi