1、15.FEM of material non-linearity5.0 非线性问题概述非线性问题概述5.1 Solution of non-linear algebraic equations2Non-linear behaviour of solids:material non-linearity Geometric non-linearityMaterial non-linearity:the loading and unloading response of the material is different.Typical here is the case of classical e
2、lasticplastic and viscoelastic behavior,non-linear elasticity.Small deformation5.0 非线性问题概述3 Finite deformation:Geometric non-linearity 大变形大变形位移与应变非线性关系位移与应变非线性关系与材料非线性一样重要。如：平板的弯曲问与材料非线性一样重要。如：平板的弯曲问题，大挠度理论分析结果更符合实际情况。题，大挠度理论分析结果更符合实际情况。大变形问题有限元分析的理论和方法存在不大变形问题有限元分析的理论和方法存在不同学派之间的争鸣。如解的稳定性问题、收同学派之间的
3、争鸣。如解的稳定性问题、收敛性问题以及效率等。敛性问题以及效率等。4Elastic-plastic Materialirrecoverable strains on load removal,showing a history dependent different path.ideal plastic behavior in a limiting yield stress.hardening/softening plastic material in which the yield stress depends on some parameter k5 Viscoelasticity his
4、tory dependence of deformationLinear models for viscoelasticityThe rate at which inelastic strains develop depends not only on the current state of stress and strain but,in general,on the full history of their development.Viscoplasticityelasto-viscoplastic material 65.1 Solution of non-linear algebr
5、aic equationsRuk0)(0)(RuuKuTnuuuuu.321iterative methodTnuuuuu00302010.00)(KuKRKu101RKunn11|to1nnuuthe solution process is path dependent75.1.2 Newtons method0)()(RuuKuNewtons method is the most rapidly convergent process for solutions of problems.0)()(RuuKu0)(1nnnnuuuunnnuuu1)(1nnnuuu)(k1Tnnnuu8nnuT
6、k)(k1Tnnnuunnnuuu1切线变刚度法切线变刚度法9The Newton process,despite its rapid convergence,has some negative features:a new KT matrix has to be computed at each iteration;if direct solution for Eq.is used the matrix needs to be factored at each iteration;1.on some occasions the tangent matrix is symmetric at a
7、 solution state but unsymmetric otherwise(e.g.in some schemes for integrating large rotation parameters or non-associated plasticity).In these cases an unsymmetric solver is needed in general.105.1.3 Modified Newtons methodthe variable jacobian matrix KT by a constant approximation0TTkknnnnuuu1)(k10
8、Tnnuu KT can be chosen as the matrix corresponding to the first iteration.Converge at a slower rate,but some of the difficulties mentioned above for the Newton process disappear.切线常刚度法115.1.4 Incremental-secant or quasi-Newton methodsthe first iteration of the preceding)(k111T1uu)()(k2112T1uuu)(k212
9、T2uu割线变刚度法12direct iterative methodsdynamic relaxation methods5.1.5 incremental and rate methodsIf the behaviour is path dependent(e.g.as in plasticity-type constitutive laws)the use of small increments is desirable to preserve accuracy in solution changes.135.2 Stress and strain theoryzzyzxyzyyxxzx
10、yxzzyzxyzyyxxzxyxijmijijsssssssssssss3iimijijmijmijijmijsxyxzzyzxyzyxxyznnxyzninvariants of the stress032213IIInnn3211zyxI1332212222zxyzxyxzzyyxI32122232xyzzxyyzxzxyzxyzyxI14Stress deviator invariants032213JsJsJs03211ssssssJzyx1332212222sssssssssssssssJzxyzxyxzzyyx3213ssssssssssssJzyzxzzyyxyzxyxx22/
11、1213232221321Jqe纯剪应力OQO213p3),(321Q偏平面等压线1rOO通过原点与三个坐标轴夹角相等的线-等压线。与OO垂直的平面称为偏平面，方程r 332115OQO213p3),(321Q偏平面等压线1rr=0的平面称为平面pIr3131)(31OO1321qJrOQ322)()()(31)(31QO22121323222121232123222122作用在平面上的应力分量：16OQO213p3),(321Q偏平面等压线1ro30113222323332M应力偏量与与 轴垂直线的夹角称为Lode角/23233sinJJ2332326632sin 221J3ssin 222
12、J3s 32sin 223J3sLode参数 tan317)34sin(sin)32sin(22321J3mmm185.3 review of Isotropic plasticity models5.3.1 Yield functionsThe yielding can occur only if the stress satisfies the general yield criterion0),(F an isotropic hardening parameter.This yield condition can be visualized as a surface in an n-di
13、mensional space of stress with the position and size of the surface dependent on the instantaneous value of the parameters k.0),(321F0),(321IIIF0),(321JJIF0),(qpF19 Elastic and ideal plastic material:It is mean destruction when the yielding occur,yielding surface is as same as destructing surface.Ha
14、rdening material:from yielding to destructing,following development of yielding surface.Assuming:the similar shape between yielding surface and critical one.For an isotropic material all functions can be represented in terms of the three stress invariants.0),(321JJIFOne useful form of these invarian
15、ts for use in yield functions is given by0),(mF2J纯剪应力201.Tresca0)(cos2YF2J当最大剪应力达到某一值时，材料发生屈服kkk222133221)34sin(sin)32sin(22321J3mmm04)(4)(4)(221322322221kkk123setting0Tresca yield criterion,Mises yield criterion2.Hubervon Mises:0)(322YF21For soils,concrete and other frictional materials the MohrCou
16、lomb or DruckerPrager surfaces is frequently used.3.MohrCoulomb:sin2cos23131c)34sin(sin)32sin(22321J3mmm0cossinsin31cossin3121cJI0cossinsin31cossincFm2Jcpqtansinsincos3sin3tansinsincos3cos3cc22c(k)and(k)are the cohesion and the angle of friction,respectively,which can depend on an isotropic strain h
17、ardening parameter k.criterion yield Tresca0 If0F3sintanPresent the minimum value of Freduce DruckerPrager yield functions.0)(3KFm2sin33sin2)sin3(3cos6K4.DruckerPrager yield functions 23 MohrCoulomb受拉破坏DruckerPrager受压破坏663sintan12324Rewrite the following formijijijmmijFFFFF00011131ijmijijijijijJJJJJ
18、JJJ2233332261913tanijijJJ2221ijijijmmJJFJJFFF3322Alternatively,we can always write:0),(mF255.3.2 Flow rule(normality principle)the plastic strain increments is related to the yield surface；the plastic strain rate components to be normal to the yieldsurface in the space of nine stress and strain dime
19、nsions.ijpijF)d(d is a proportionality constant,as yet undetermined,often referred to as the plastic consistency parameter.During sustained plastic deformation we must have0)(and0dFpijeijijpijeijijdddThe normality principle during elastic loading/unloading0dand0dF26specifying separately a plastic fl
20、ow rule potential),(ijQQ 0d,dijpijQ加载卸载准加载卸载准则用以判别从塑性状态出发是继续塑性加载塑性加载还是弹性卸载弹性卸载，这是计算过程中判定是否继续塑性变形塑性变形以及决定是采用弹塑性本构关系弹塑性本构关系还是弹性本构关系弹性本构关系所必须的。Q=F，associative plasticity.Q F，non-associative.The plastic flow rule potential of MohrCoulomb:constantsinsin31cossinmQ00无剪涨最大剪涨27ijijijijFAFd)d(ddLoadingijijFAd
21、1)d(A:coefficient related to hardening/softening rules5.3.3 Hardening/softening rulesDefine hardening(or softening)of the plastic behaviour of the material.Isotropic hardeningThe evolution of,governing the size of the yield surface,is commonly related to the rate of plastic work or directly to the c
22、onsistency parameter.If related to the rate of plastic work has dimensions of stress and a relation of the type.28ddTppijijdijpijQ dd0ddFFdFijijAssuming0dddddijijijijpijijijijQFFFFF0),(FijijFAd1d01dijijijijijijQFAFF0dijijijijijijQFFFAijijQFA295.3.4 Plastic stressstrain relationsThe strains are assum
23、ed to be divisible into elastic and plastic parts given aspeijijijpeijijijdddeijijDddijeDijdd1ijpijQ ddijijijQDddd1the elastic stress changesijijDeddThe plastic strain(rate)will occur only if the elastic stress changes，tend to put the stress outside the yield surface,that is,is in the plastic loadin
24、g direction.If,on the other hand,this stress change is such that unloading occurs then of course no plastic straining will be present.When plastic loading is occurring the stresses are on the yield surface given by Flow rule.30ddd1QDT312312332211dddddddTQQQQQQQ312312332211T312312332211ddddddddddQDDi
25、jijFAd1dd1dTFA d1ddQDFADTdddQDFADATddQDFADATddQDFADAT31eeepQDFADADTddepDdddpeDDRewritepeepDDDepepDDDeeeQDFADQFDTT该式对硬化、软化和理想塑性材料均适用，也适用于非耦合下相关联与非关联流动法则。The continuum elasto-plastic matrix is symmetric only when plasticity is associative and when kinematic hardening is symmetric.In general,non-associ
26、ative materials present stability difficulties,and special care is needed to use them effectively.Similar difficulties occur if the hardening moduli are negative which,in fact,leads to a softening behaviour.For ideal plasticity A are zero.325.4 Elastic-Plastic FE equations in increment formTraction
27、componentsDisplacement components)(tui)(tTibody force components)(tbi02dTudbudTTT0)()()(2dttTudttbudttTTTTime t+tuuuDD)()(ttt)()(btbttb)()(TtTttT)()(ututtu)()(ttt33dtTudtbudtdTudbudTTTTTT)()()(22epuDTTTuudTudbuduDu2TTepTTdtTudtbudtuTTTT)()()(2Dividing the domain into NE finite elements,NP nodes.eNod
28、euNueeNodeuNueTTTeNuueNodeNB34 eeeeeedTNudbNuduBDBu2TTNodeTTNodeNodeepTTNode)()()(2TTNodeTTNodeTTNodedtTNudtbNudtBu eeeeeedTNdbNudBDB2TTNodeepTdtTNdtbNdtB)()()(2TTTIf the components of stress at t satisfy the equilibrium equations,then02TTNodeepT eeeeeedTNdbNudBDB35else eeeeeedttTNdttbNudBDB)()(2TTN
29、odeepTdtB)(TeeeeeetQttQuK)()(Nodeep)()(NodeeptQttQuKThis simple and general description of solution scheme:For t+t,give loading Increments pCompute nodal u using FEM equationsu(t+t)=u(t)+uCompute compute1.adjust non-equilibrium forces)()(ttt36dttBdttTNdttbNttQttQttQttQeeeeeeeee)()()()()()()(TTTexter
30、nalinternalexternalinternal0)()(externalinternalttQttQIfReturn 1;Else,start next loading Increments.375.5 Computation of stress incrementsWe have emphasized that with the use of iterative procedures within a particular increment of loading,it is important always to compute the stresses askkttt)()(co
31、rresponding to the total change in displacement parameterskukkuB kkikudu0stepwhich has accumulated in all previous iterations within the step.This point is of considerable importance as constitutive models with path dependence have different responses for loading and unloading.The stresses have to b
32、e integrateddDtttkk0ep)()(incorporating into Dep the dependence on variables in a manner corresponding to a linear increase of or .kku385.5.1 Judge elastic or plastic stress stateFor each increment or iterative step,first compute strains increment after obtaining displacements increment.kkuB Second,
33、predicting stresses increment(the trial value)according to elastic relationskkDekkttt)()(Third,computing value of yielding function),()(ttFttFk0)(and0)(IftFttFElastic state0)(and0)(IftFttFElastic unloading state0)(and0)(IftFttFplastic loading state0)(and0)(IftFttFDividing stresses increment to two p
34、artskkmm)1(one accord to elastic relations and other accord to plastic relations.39For plastic statekdDk0epThe procedures for integration of equation can be classified into explicit and implicit categories.5.5.2 The procedures for integration of the constitutive equation Explicit methodsSubincrement
35、ation-A direct integration process is used Some form of the RungeKutta process is adopted.1.SubincrementationThe known increment is subdivided into m intervals and the integral of the equation is replaced by direct summation,k10ep1mjkjkDm denotes the tangent matrix computed for stresses and hardenin
36、gParameters updated from the previous increment in the sum.jD ep40Such integration will generally result in the stress change departing from the yield surface by some margin.2.RungeKutta method A more precise explicit procedureFirst an increment of is applied in a single-step explicit manner to obta
37、in using the initial elasto-plastic matrix.2/212/1DThis increment of stress is evaluated to compute D1/2 and finally we evaluate2/121DThis process has a second-order accuracy and,in addition,can give an estimate of errors incurred as2/12If such stress errors exceed a certain norm the size of the inc
38、rement can be reduced.41F.B.Hildebrand.Introduction to Numerical Analysis.Dover Publishers,2nd edition,1987.G.C.Nayak and O.C.Zienkiewicz.Note on the alpha-constant stiffness method for the analysis of nonlinear problems.International Journal for Numerical Methods in Engineering,4:579582,1972.42Rate
39、-dependent plasticity behaviors 一、一、Nonlinear model Abaqus1 Brittle cracking Isotropic elasticity and brittle shear 2 Modified Drucker-Prager/Cap plasticity Drucker-Prager/Cap plasticity hardening and isotropic elasticity or porous elasticity 3 Cast iron plasticity Cast iron compression hardening,ca
40、st iron tension hardening,and isotropic elasticity 5.6 Abaqus 43This option is used to specify the plastic part of the material behavior for elastic-plastic materials that use the extended Cam-clay plasticity model.#Logarithmic plastic bulk modulus#Stress ratio at critical state M#the initial yield
41、surface size#the parameter defining the size of the yield surface on the“wet”side of critical state 0a#the ratio of the flow stress in triaxial tension to the flow stress in triaxial compression K 4 Cam-clay plasticity Elasticity or porous elasticity 445 Concrete Isotropic elasticity to define the p
42、roperties of plain concrete outside the elastic range,field variable dependencies included in the definition of the compressive yield stress,in addition to temperature.Concrete damaged plasticity#Concrete compression hardening,concrete tension stiffening,and isotropic elasticity#define flow potentia
43、l,yield surface,and viscosity parameters for the concrete damaged plasticity material model.$Dilation angle,$Flow potential eccentricity,defines the rate at which the hyperbolic flow potential approaches its asymptote.45$the ratio of initial equibiaxial compressive yield stress to initial uniaxial c
44、ompressive yield stress$Kc,the ratio of the second stress invariant on the tensile meridian,q(TM),to that on the compressive meridian,q(CM),at initial yield for any given value of the pressure invariant p such that the maximum principal stress is negative.Viscosity parameter,used for the visco-plast
45、ic regularization of the concrete constitutive equations in Abaqus/Standard analyses.7 Crushable foam plasticity Crushable foam hardening and isotropic elasticity the material behavior for elastic-plastic materials that use the crushable foam plasticity model.468 Drucker-Prager plasticity Drucker-Pr
46、ager hardening and isotropic elasticity or porous elasticity Material angle of friction,in the pt plane.K,the ratio of the flow stress in triaxial tension to the flow stress in triaxial compression.Dilation angle,in the pt plane.9 Mohr-Coulomb plasticity Mohr-Coulomb hardening and isotropic elastici
47、ty.define the yield surface and flow potential parameters for elastic-plastic materials that use the Mohr-Coulomb plasticity model.Friction angle,Dilation angle 4710 Metal plasticity Elasticity or hyperelasticity.Isotropic elasticity,orthotropic elasticity(requires anisotropic yield),hyperelasticity
48、,or the combination of an equation of state and isotropic linear elastic shear behavior for an equation of state.use the Mises or Hill yield surface Rate-dependent plasticity behaviors Cap creep:Elasticity,modified Drucker-Prager/Cap plasticity,and Drucker-Prager/Cap plasticity hardening.Metal creep
49、 Drucker-Prager creep,Metal plasticity,Rate-dependent viscoplasticity,Time-dependent volumetric swelling 48二、二、Characterizing elementsFive aspects of an element characterize its behavior#Family#Degrees of freedom#Number of nodes#Formulation#Integration 49Degrees of freedomThe degrees of freedom are
50、the fundamental variables calculated during the analysis.For a stress/displacement simulation the degrees of freedom are the translations and,for shell and beam elements,the rotations at each node.For a heat transfer simulation the degrees of freedom are the temperatures at each node;for a coupled t
