Elements

13 Pages·2015·4.33 MB·English

Hindawi Publishing Corporation
Mathematical Problems in Engineering
Volume 2015, Article ID 457046, 12 pages
http://dx.doi.org/10.1155/2015/457046
Research Article
Stress Field Gradient Analysis Technique Using
𝐶0
Lower-Order Elements
JianweiXingandGangtieZheng
SchoolofAerospaceEngineering,TsinghuaUniversity,Beijing100084,China
CorrespondenceshouldbeaddressedtoGangtieZheng;[email protected]
Received1September2014;Accepted30September2014
AcademicEditor:SongCen
Copyright©2015J.XingandG.Zheng.ThisisanopenaccessarticledistributedundertheCreativeCommonsAttributionLicense,
whichpermitsunrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalworkisproperlycited.
Forevaluatingthestressgradient,amathematicaltechniquebasedonthestressfieldoflower-order𝐶0elementsisdevelopedin
thispaper.Withnodalstressresultsandlocationinformation,anoverdeterminedandinconsistentequationofstressgradientis
establishedandtheminimumnormleastsquaressolutionisobtainedbytheMoore-Penrosepseudoinverse.Thistechniquecanbe
appliedtoanyelementtypeincomparisonwiththesuperconvergentpatch(SCP)recoveryforthestressgradient,whichrequires
thequadraticelementsatleastandhastoinverttheJacobiandHessianmatrices.Theaccuracyandvalidityofthepresentedmethod
aredemonstratedbytwoexamples,especiallyitsmeritofachievinghighaccuracywithlower-orderlinear𝐶0elements.Thismethod
canbeconvenientlyintroducedintothegeneralfiniteelementanalysisprogramsasapostprocessingmodule.
1.Introduction the𝐶0-1 enhancedpatchtestofaconvergencecondition.In
addition to the reason that only a few types of 𝐶1 elements
Effects and evaluations of the stress gradient have always are successful and effective, the limitation in the element
been important to the design of an engineering structure. shape and the number of DOFs all impede the practical
As nowadays the stress is usually calculated with a finite applications.
element method (FEM), the calculation technique of the In order to circumvent problems of the 𝐶1 element,
stressgradienthasbeentheeffortofmanyresearchworks. a mixed or coupled formulation using lower-order ele-
The size effect on the stress/strain gradient has already ments (𝐶0) is implemented. Zhao et al. [13] proposed a
beenmanifestedinmanyexperiments[1–4].Toovercomethe mixedquadrilateralelement(CQ12+RDKQ)forthecouple
problem of the conventional elasticity theory in describing stress/strain gradient elasticity, in which one satisfies 𝐶0
theeffectofsize,thecouplestress[5–7]andstraingradient continuitytocomputestrainsandtheotheronesatisfiesweak
plasticity [8, 9] theories have been developed and proved 𝐶1 continuitytocalculatestraingradients.Amanatidouand
effectiveinintroducingstraingradienttermsintoconstitutive Aravas [14] developed an alternative mixed finite element
equations and yield functions. However, according to these formulation and studied the Mindlin [15] strain gradient
theories, the element should be at least 𝐶1 continuity. This elasticity theories in detail. The 𝐶1 continuous Hermitian
complicates the element construction. Dasgupta and Sen- shape functions for the plastic and damage multipliers and
gupta[10]developedan18-DOF(degree-of-freedom)trian- the 𝐶0 continuousinterpolation functions for displacement
gular plate bending element, whose displacement function field were employed by Dorgan and Voyiadjis [16]. For
is fifth-order polynomial in the area coordinates. Zervos et reasonssuchastheadditionalandinevitablecomputational
al. [11] proposed a 𝐶1 continuity elastoplasticity triangular cost and only a few available element types, the mixed
element with assumption that the stress increment consists elementsareinconvenienttopracticalapplications.
of both the strain increment and its Laplacian. Chen and Anothermotivationforthestressgradientanalysisisthe
Li [12] put forward a 24-DOF quadrilateral spline element fatigue life calculation and the fatigue behavior prediction.
forthecouplestress/straingradientelasticity,whichcanpass Neuber[17]developedawidelyadoptedmethodinpredicting
2 MathematicalProblemsinEngineering
thenotchfatiguelife;nevertheless,themethodonlyfocuses TheMoore-Penrosepseudoinverse,whichiscomputedwith
on the stress of dangerous point but without considering thesingularvaluedecomposition,isemployedtoobtainthe
the influence of stress gradient and multiaxial effects. Even minimum norm least squares solution [32, 33] of the stress
if the stresses of dangerous points for various structures gradient.Apopular2Dexampledemonstratesitssatisfactory
are the same, the fatigue life and behavior are different accuracyincomparisonwiththeSCPrecovery.Thepresented
becauseofdifferentstressgradientdistribution[18].Taking methodespeciallycanalsoachievehighaccuracywithlower-
into account the stress gradient for assessing the fatigue orderlinear𝐶0elements.Another3Dexampleshowsthatthe
strength [19], the stress field intensity method [20] and the method can be easily introduced into general FEAP (such
advancedvolumetricmethod[21]aretwoofthemostpopular as ANSYS, PATRAN, and ABUQUS) as a postprocessing
techniques. Although the stress field intensity method is module.
simple, it needs experiment to provide the size of notch This paper is organized as follows. First, the super-
damaged zone. The advanced volumetric method requires convergent patch recovery method is briefly introduced in
theaccuratecalculationofstressgradientsnearnotchroots; Section2. Then, the stress filed gradient analysis technique
however, it is usually difficult to guarantee the calculation and the error estimation are presented in Section3. In
accuracy. Section4,twonumericalexamplesareemployedtoevaluate
Thetechniqueoffiniteelementpostprocessingisanother theperformanceofthepresentedmethod.Thelastsectionof
approach to calculating continuous stress and its gradient. thepaperistheconclusions.
The superconvergent patch (SCP) recovery is known as
the most effective technique [22], which can obtain higher
accuracy than any other points at the Gauss integration 2.SCPStressGradientRecovery
point within an element. Utilizing the patches of elements
2.1.SecondDerivativesinFEM. Takingtheisotropicnormal
in FEM, Zienkiewicz and Zhu [23, 24] obtained continu-
ous stress field with data at the superconvergent points. A stress 𝜎𝑥, for example, the gradient vector G(𝜎𝑥) can be
writtenas
detailed analysis of the SCP recovery technique and error
estimation was presented in [25, 26]. Since then, Wiberg
andAbdulwahab[27],Wibergetal.[28],andParkandShin G(𝜎𝑥)=[𝑔𝑥(𝜎𝑥),𝑔𝑦(𝜎𝑥)]
[29] improved the SCP approach using the method of least
square for fitting the residuals of the equilibrium equation 𝐸 𝑛 𝜕2𝑁 𝜕2𝑁
andapplyingtheprincipleofvirtualwork.Recently,theSCP = 0 [∑( 𝑖𝑢 +𝜇 𝑖V),
technique has been extended to recover the higher-order (1−𝜇0) 𝑖=1 𝜕𝑥2 𝑖 0𝜕𝑥𝜕𝑦 𝑖 (1)
derivatives.UsingQ8(8-nodequadratic,9-pointquadrature)
𝑛 𝜕2𝑁 𝜕2𝑁
quadrilateralelement,T6(6-nodequadratic,4-pointquadra- ∑( 𝑖𝑢 +𝜇 𝑖V)],
ture) triangular element, and the mixed mesh of Q8 and 𝜕𝑦𝜕𝑥 𝑖 0 𝜕𝑦2 𝑖
𝑖=1
T6, Gan and Akin [30] developed an element patch based
superconvergentsecondderivativerecoverytechniquefora
lower-orderstraingradientplasticitymodel.Banketal.[31] where𝑔𝑥 and𝑔𝑦 arethecomponentsofthestressgradient,
proposedasuperconvergentderivativerecoveryschemewith 𝑁𝑖 istheshapefunction,𝑛isthenumberofelementnodes,
the Lagrange triangular element and the superconvergence and(𝑢𝑖,V𝑖)isthenodaldisplacementvector;𝐸0 = 𝐸,𝜇0 = 𝜇
estimation of the 𝑝th derivatives. Nevertheless, the higher- fortheplanestress,𝐸0 = 𝐸/(1−𝜇2),𝜇0 = 𝜇/(1−𝜇)forthe
order derivatives such as stress gradients also need higher- planestrain;𝐸and𝜇areYoung’smodulusandPoisson’sratio,
order element shape functions, which lead to additional respectively.
DOFsandcalculationcost,especiallyforcomplexstructures Firstly,thesecondderivativesoftheshapefunctionsare
orlargesystems. necessary for evaluating the stress gradient, but they are all
Therefore, the effort is still required to enhance tech- equal to zero for lower-order linear 𝐶0 elements, such as
niques for calculating stress gradient in such aspects as Q4(4-nodequadrilateralelement)andT3(3-nodetriangular
reducingcomputationcost,developingmoreelementtypes, element).
and providing codes for engineering applications. In the Assumingusingparametricelements,thefirstderivative
rest of this paper, a mathematical technique based on the transformationsfromtheparametriccoordinatetothephys-
stress field results of general finite element analysis pro- icalcoordinateare
grams (FEAP) is proposed to calculate the stress gradient.
Because it works with the stress field results, this method
𝜕 𝜕𝑥 𝜕𝑦 𝜕
cttoeacnhpnhbyieqsuiccoea.nlAscisdoneororedtdirnaaanstsefaoirssmtrreaeqtsisuoinfireefdlrdo,mtghripasadrtieeacnmhten(tiSrqiFcuGeco)doaornedasilnynasotiest {{{{{{𝜕𝜕𝜉}}}}}}=[[[[𝜕𝜕𝑥𝜉 𝜕𝜕𝑦𝜉]]]]{{{{{{𝜕𝜕𝑥}}}}}}, (2)
involve Jacobi and Hessian matrices and hence does not {𝜕𝜂} [𝜕𝜂 𝜕𝜂]{𝜕𝑦}
need to calculate their inverses. With the location and
stress information of the central node and its surrounding
samplenodes,anoverdeterminedandinconsistentequation where𝜉and𝜂aretheparametriccoordinates,andthematrix
of the stress gradient at the central node is established. istheso-calledJacobimatrix.
MathematicalProblemsinEngineering 3
y
x
Global coordinates
(a) Q4 (b) Q8
Gauss quadrature points Gauss quadrature points
Stress sample points Stress sample points
Central points Central points
(c) T3 (d) T6
Figure1:Constructionof2Dpatcheswithquadrilateralandtriangularelements.
Therelationshipbetweenthesecondparametricderiva- wherethesecondsquarematrixontherighthandsideisthe
tivesandthesecondphysicalderivativescanbewrittenas Hessianmatrix.From(3),thephysicalsecondderivativescan
be expressed by the first and second parametric derivatives
andtheinverseoftheJacobiandHessianmatrices.
𝜕2 𝜕2𝑥 𝜕2𝑦
{{{{{{{{{ 𝜕𝜕𝜉22 }}}}}}}}} [[[[ 𝜕𝜕2𝜉𝑥2 𝜕𝜕2𝜉𝑦2 ]]]]{{{𝜕𝜕𝑥}}} t2i.o2n.SpCoPinRtsecaorevearlysoanthdeEsrurpoerrEcostnimveartgieonnt.pThoinetGsafourssthinetsetgrreass-
=[ ]
{{{{{{{{{ 𝜕𝜕𝜂22 }}}}}}}}} [[[[ 𝜕𝜕2𝜂𝑥2 𝜕𝜕2𝜂𝑦2 ]]]]{{{{𝜕𝜕𝑦}}}} gsmhruaodswtinesnaitnt,isFafniygduqrtuehased1r(sabut)picaenracdto1nl(evdae)sr.tgTheinnettephlaeetmcSheCenPstscshatanrpesebsefugbnruacditliitoennasst
{𝜕𝜉𝜕𝜂} [𝜕𝜉𝜕𝜂 𝜕𝜉𝜕𝜂] recoverybecauseoftherequirementontheexistenceofthe
secondderivatives.
[(𝜕𝑥)2 (𝜕𝑦)2 2𝜕𝑥𝜕𝑦 ] WiththeSCPrecovery,therecoveredstressgradient𝑔𝑥∗,𝑦
[[ 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 ]] canbeexpressedas
+[[[[[(𝜕𝜕𝑥𝜂)2 (𝜕𝜕𝑥𝜂)2 2𝜕𝜕𝑥𝜂𝜕𝜕𝑦𝜂 ]]]]] (3) 𝑔𝑥∗,𝑦 =Pa, (4)
[𝜕𝑥𝜕𝑥 𝜕𝑦𝜕𝑦 𝜕𝑥𝜕𝑦 𝜕𝑥𝜕𝑦]
+
[ 𝜕𝜉 𝜕𝜂 𝜕𝜉𝜕𝜂 𝜕𝜉𝜕𝜂 𝜕𝜂𝜕𝜉] where P is the polynomial vector and a is the coefficient
columnvectortobedetermined.
𝜕2
×{{{{{{{{{ 𝜕𝜕𝑥22 }}}}}}}}}, ForquadraticPel=em[1ent𝑥sli𝑦ke𝑥Q28𝑥an𝑦d𝑦T26],,Pisassumeda(s5)
{{{{{{{{{ 𝜕𝜕𝑦22 }}}}}}}}} where the recovered stress gradient has the same order of
{𝜕𝑥𝜕𝑦} displacementfiled[29],andthena=[𝑎1 𝑎2 ⋅⋅⋅ 𝑎6]𝑇.
4 MathematicalProblemsinEngineering
The functional integral of the difference between the
stressgradientcalculatedfromtheSCPrecovery𝑔∗ andthe
𝑥,𝑦
stressgradientobtainedbytheFEM𝑔̃𝑥,𝑦canbeexpressedas
y
𝑚
𝐹(𝑔∗ ,𝑔̃ )= ∑∫ (𝑔∗ −𝑔̃ )2𝑑𝑆
𝑥,𝑦 𝑥,𝑦 𝑥,𝑦 𝑥,𝑦
𝑒=1 𝑆𝑒 x
(6) z
𝑚
2
= ∑∫ (Pa−𝑔̃𝑥,𝑦) 𝑑𝑆,
𝑒=1 𝑆𝑒
where𝑚isthenumberofpatchelementsand𝑆𝑒istheelement
area.
Stress sample points
The assumed coefficient vector a can be obtained by
minimizingthefunctionalintegralin(6),whichyields Central points
𝑚 𝑚 Figure2:Constructionof3Dpatchwithhexahedralelement.
∑∫ P𝑇P𝑑𝑆a= ∑∫ P𝑇𝑔̃𝑥,𝑦𝑑𝑆, (7)
𝑒=1 𝑆𝑒 𝑒=1 𝑆𝑒
wheresuperscript𝑇denotesthetranspose,andbothsidesof In this section, a simple and convenient mathematical
technique for calculating the stress gradient will be devel-
(7)canbecalculatedbytheGaussintegral:
oped.Thismethodissuitableforanyelementtypesandonly
𝑛 𝑛 needstheinformationofnodallocationandstress.
∑P𝑇(𝑥𝑖,𝑦𝑖)P(𝑥𝑖,𝑦𝑖)𝑤𝑖a=∑P𝑇(𝑥𝑖,𝑦𝑖)𝑔̃𝑥,𝑦(𝑥𝑖,𝑦𝑖)𝑤𝑖,
𝑖=1 𝑖=1
(8) 3.1.GeneralFormulation. AsshowninFigure1for2Dprob-
lem, the central point 𝑃𝑐(𝑥0,𝑦0) is surrounded by sample
where 𝑛 = 𝑚 ⋅ 𝑛𝑔, 𝑛𝑔 is the number of Gauss points per points 𝑃𝑗(𝑥𝑗,𝑦𝑗), 𝑗 = 1,2,...,𝑛𝑠 (𝑛𝑠 is number of sample
element,𝑤𝑖istheweightcoefficientofeachintegrationpoint, points).Thevectorfrom𝑃𝑐 to𝑃𝑗 isn𝑐𝑗 = (𝑥𝑗 −𝑥0,𝑦𝑗 −𝑦0);
and (𝑥𝑖,𝑦𝑖) is the global coordinates of each Gauss point. thus,theunitvectorcanbewrittenasn𝑐𝑗 = (1/‖𝑃𝑐𝑃𝑗‖)(𝑥𝑗 −
Finally,thereis 𝑥0,𝑦𝑗 −𝑦0),where‖𝑃𝑐𝑃𝑗‖ = √2(𝑥𝑗−𝑥0)2+(𝑦𝑗−𝑦0)2 isthe
a=C−1D, (9) vectorlength.
Thedirectionalderivativealongn𝑐𝑗atthecentralpoint𝑃𝑐
inwhich canbeexpressedas
𝑛
C=∑𝑖=1P𝑇(𝑥𝑖,𝑦𝑖)P(𝑥𝑖,𝑦𝑖)𝑤𝑖, (10) 𝜕𝜕n𝜎̃𝑐𝑥𝑗 (𝑃𝑐)=𝑃𝑃l→i∈mn𝑃𝑐𝜎̃󵄩󵄩󵄩󵄩𝑥𝑃𝑃−𝑐𝑃𝜎̃󵄩󵄩󵄩󵄩𝑥𝑃𝑐
𝑛 𝑐𝑗
D= ∑P𝑇(𝑥𝑖,𝑦𝑖)𝑔̃𝑥,𝑦(𝑥𝑖,𝑦𝑖)𝑤𝑖. 1 𝜕𝜎̃
𝑒=1 = 󵄩󵄩󵄩󵄩󵄩𝑃𝑐𝑃𝑗󵄩󵄩󵄩󵄩󵄩[(𝑥𝑃𝑗 −𝑥𝑃𝑐) 𝜕𝑥𝑥 (𝑃𝑐) (11)
In[25,26],utilizingthequadraticcompletenesspolyno-
mial, the SCP recovery for stress gradient has the order of +(𝑦 −𝑦 )𝜕𝜎̃𝑥 (𝑃)],
accuracy𝑂(ℎ3)(ℎistheelementcharacteristiclength).And 𝑃𝑗 𝑃𝑐 𝜕𝑦 𝑐
then, it is easy to verifythat the errorofthe stress gradient
recoveryisoforder𝑂(ℎ2)bythesameprocedure.Ithasbeen where 𝜎̃𝑥 is the stress result of commercial finite element
found that the recovered variable field is more accurate at software.Thelimitcanbeapproximatelywrittenas
nodepoints.
3.StressFieldGradientAnalysis 𝑃𝑃l→i∈mn𝑃𝑐𝜎̃󵄩󵄩󵄩󵄩𝑥𝑃𝑃−𝑐𝑃𝜎̃󵄩󵄩󵄩󵄩𝑥𝑃𝑐 ≈ 𝜎̃󵄩󵄩󵄩󵄩󵄩𝑥𝑃𝑃𝑗𝑐−𝑃𝑗𝜎̃󵄩󵄩󵄩󵄩󵄩𝑥𝑃𝑐. (12)
𝑐𝑗
It is known that the SCP recovery for the stress gradient
needs quadratic elements at least, and the shape functions Substituting(12)into(11)yields
should be determined beforehand to compute the second
derivatives of the displacement field. All of this leads to n𝑐𝑗⋅G𝑃(𝜎̃𝑥)𝑇 =𝜎̃𝑥𝑃𝑗 −𝜎̃𝑥𝑃𝑐. (13)
𝑐
an additional calculation cost and difficulty, especially to
thoseengineerswhoareusedtoworkwithcommercialfinite Forallsamplepoints,thestressgradientequationcanbe
elementsoftware,andprovidingthepreliminaryanalysisby writtenas
lower-orderlinearelements,wheretheSCPtechniquecannot
beimplementedyet. A⋅G𝑃(𝜎̃𝑥)𝑇 =b, (14)
𝑐
MathematicalProblemsinEngineering 5
y
y
C
B
B C
P0 A P0 5r0 P0
r0 E D x
A
O E D x
r0 10r0
(a) (b)
Figure3:Aninfiniteplatewithacentralholesubjectedtoaremotelyuniformtension:(a)theorymodel,(b)finiteelementmodel.
(a) Quadrilateralcoarsemesh400elements: (b) Quadrilateral fine mesh 1600 elements: (c) Quadrilateraldouble-finemesh6400ele-
441nodes(Q4),1281nodes(Q8) 1681nodes(Q4),4961nodes(Q8) ments:6561nodes(Q4),19521nodes(Q8)
(d) Trianglecoarsemesh564elements:323 (e) Triangle fine mesh 2220 elements: 1191 (f) Triangle double-fine mesh 9040 ele-
nodes(T3),1209nodes(T6) nodes(T3),4601nodes(T6) ments:4681nodes(T3),18401nodes(T6)
Figure4:Meshesofquadrilateralandtriangularelements.
where In general, 𝑛𝑠 > 2; hence, the coefficient matrix A is
notsquareand(14)isoverdeterminedandinconsistent.An
n𝑐1 [𝜎̃𝑥𝑃1 −𝜎̃𝑥𝑃𝑐] ePffenecrtoisveepwseauydtooinsovlevresethoefAeq𝑛𝑠u×a2.tiThoneipssteoudemoipnlvoeyrstehMeM2×𝑛o𝑠ocraen-
[ ] [ ] 𝑛×2
A=[[[n...𝑐2]]], b=[[[[𝜎̃𝑥𝑃2 −... 𝜎̃𝑥𝑃𝑐]]]]. (15) tbheactaislc,ulated by the singular value decomposition of A 𝑠 ;
[n𝑐𝑛𝑠] [𝜎̃𝑥𝑃𝑛𝑠 −𝜎̃𝑥𝑃𝑐] A=UΣV𝑇, (16)
6 MathematicalProblemsinEngineering
200 200 250
200
150 150
150
Error (%) 10500 CoLairnsee- m45e∘sh Error (%) 10500 FLininee m-4e5s∘h Error (%) 100 D-Lfiinnee -m45e∘sh
50
0 0 0
−50 −50 −50
0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6
r/r r/r r/r
0 0 0
(a) (b) (c)
150 75 150
100 50
100
Line-Y
50 25
%) %) %) D-fine mesh
or ( or ( or ( 50
Err 0 Err 0 Err
−50 Line-Y −25 Line-Y 0
Coarse mesh
Fine mesh
−100 −50 −50
0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6
r/r0 r/r0 r/r0
SFG-Q4 SFG-Q4 SFG-Q4
SFG-Q8 SFG-Q8 SFG-Q8
SCP-Q8 SCP-Q8 SCP-Q8
(d) (e) (f)
Figure5:Stressgradientsofsamplelineswithquadrilateralelements.
whereU ∈ 𝑅𝑛𝑠×𝑛𝑠 andV ∈ 𝑅2×2areorthogonalmatrices,Σ ∈ 3.2. Error Estimation. The Moore-Penrose pseudoinverse
𝑅𝑛𝑠×2isadiagonalmatrixwithnonnegativerealnumbers,and M𝑛0×𝑛𝑠 (𝑛0 = 2 or 3 for 2D or 3D case) of A𝑛𝑠×𝑛0 ∈ 𝑅𝑛𝑠×𝑛0
thediagonalentriesareknownasthesingularvalues.Then, isuniqueandsatisfiesallofthefollowingfourcriteria:
M=VΣ†U𝑇, (17)
AMA=A,
whereΣ†isthetransposeofΣanditsdiagonalentriesarethe MAM=M,
inversesofthosenonzerodiagonalelementsofΣ.
(19)
Finally,thestressgradientatthecentralpoint𝑃𝑐 canbe (AM)𝑇 =AM,
expressedby
(MA)𝑇 =MA,
G𝑃(𝜎̃𝑥)𝑇 =M⋅b. (18)
𝑐
Asthesameprocedureas2Dproblem,thestressgradient whichyields
equationof3Dproblemcanalsobewrittenas(14),inwhich
𝑥th0e,𝑦d𝑗ir−ec𝑦t0io,𝑧n𝑗v−ec𝑧t0o)raannddGst𝑃re(s𝜎̃s𝑥g)r=ad(i𝜕e𝜎ñ𝑥t/v𝜕e𝑥ct,o𝜕r𝜎̃a𝑥r/e𝜕𝑦n,𝑐𝑗𝜕𝜎=̃𝑥(/𝑥𝜕𝑗𝑧−), (AMA)𝑇 =A𝑇(AM)𝑇 =A𝑇AM=A𝑇,
𝑐 (20)
respectively.Asanexample,the3Delementpatchisshown
(MAM)𝑇 =M𝑇(MA)𝑇 =M𝑇MA=M𝑇.
inFigure2with8-nodelinearhexahedralelements.
MathematicalProblemsinEngineering 7
200 220000 115500
150 115500
110000
Line-45∘
%) 100 Line-45∘ %)%) 110000 Fine mesh %)%) Line-45∘
or ( Coarse mesh or (or ( or (or ( D-fine mesh
Err 50 ErrErr 5500 ErrErr 5500
0 00
00
−50 −−5500
0 1 2 3 4 5 6 00 11 22 33 44 55 66 00 11 22 33 44 55 66
r/r rr//rr rr//rr
0 00 00
(a) (b) (c)
150 100 50
100
50
50
%) %) %)
Error ( 0 Error ( 0 Error ( 0
−50 Line-Y
−50 Line-Y
Line-Y D-fine mesh
Coarse mesh
Fine mesh
−100 −100 −50
0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6
r/r r/r r/r
0 0 0
SFG-T3 SFG-T3 SFG-T3
SFG-T6 SFG-T6 SFG-T6
SCP-T6 SCP-T6 SCP-T6
(d) (e) (f)
Figure6:Stressgradientsofsamplelineswithtriangularelements.
y
z
o x
h
d W=20mm
H=20mmmm
L
=
100 1692 nodes
m
m
1220 elements
F=1600N
(a) (b)
Figure7:Geometry,loadingcondition,andmeshof3Dproblem.
8 MathematicalProblemsinEngineering
Table1:Stressgradientsatpoint𝑃withdifferentmeshesandmethods.
Numericalsolution Exactsolution Error
Mesh-Method-Element
𝑆𝑥𝑥/108 𝑆𝑥𝑦/108 𝑆𝑥𝑥/108 𝑆𝑥𝑦/108 𝑆𝑥𝑥 𝑆𝑥𝑦 Average
C-SFG-Q4 −10.3667 4.0309 2.3524% 8.2482% 5.3003%
C-SCP-Q8 −8.7140 3.3071 −13.9646% −11.1893% 12.5769%
C-SFG-Q8 −1.0908 4.1159 −10.1284 3.7238 7.6936% 10.5311% 9.1123%
C-SFG-T3 −9.2718 4.8073 −8.4575% 29.0992% 18.7783%
C-SCP-T6 −10.1010 3.9277 −0.2813% 5.4765% 2.8789%
C-SFG-T6 −10.2858 3.6292 1.5533% −2.5402% 2.0467%
F-SFG-Q4 −10.1920 3.6648 2.2362% 0.5960% 1.4161%
F-SCP-Q8 −9.54347 3.2959 −4.2689% −9.5305% 6.8997%
F-SFG-Q8 −10.4925 3.8314 −9.9690 3.6431 5.2506% 5.1685% 5.2096%
F-SFG-T3 −10.1635 3.7822 1.9511% 3.8174% 2.8842%
F-SCP-T6 −10.2152 3.6521 2.4697% 0.2479% 1.3588%
F-SFG-T6 −10.1343 3.5049 1.6575% −3.7928% 2.7251%
DF-SFG-Q4 −10.1587 3.6648 2.6837% 1.6590% 2.1714%
DF-SCP-Q8 −9.8485 3.3578 −0.4519% −6.8582% 3.6551%
DF-SFG-Q8 −10.2297 3.5939 −9.8932 3.6050 3.4010% −0.3094% 1.8552%
DF-SFG-T3 −9.9612 3.7442 0.6875% 3.8603% 2.2739%
DF-SCP-T6 −10.1757 3.5357 2.8555% −1.9237% 2.3896%
DF-SFG-T6 −10.1554 3.5479 2.6497% −1.5834% 2.1166%
Meshes:C:coarse;F:fine;DF:double-fine.
Methods:SFG:stressfieldgradient;SCP:superconvergentpatch.
Forallg∈𝑅𝑛0×1,thereis Apparently,theconstraintconditionsof(25)and(26)are
󵄩󵄩󵄩󵄩Ag−b󵄩󵄩󵄩󵄩=󵄩󵄩󵄩󵄩Ag−b+AMb−AMb󵄩󵄩󵄩󵄩 thesameasthefourcriterialistedin(19),sothatthematrixM̂
mustbetheMoore-PenrosepseudoinverseM.Inthisway,it
=‖AMb−b+Aw‖ (21) isshownthattheMoore-Penrosepseudoinversecanprovide
theuniqueminimumnormleastsquaressolutionof (14).
=‖(AM−I)b+Aw‖, The error of (14) consists of two parts, one is the
where w=g−Mb ∈ 𝑅𝑛0×1. With the first equation of (20), limit approximation of order 𝑂(ℎ) in (12) and the other
A𝑇(AM−I)=Oandthusthereisarelationship.Consider one is the stress field 𝜎̃𝑥. The stress calculated with FEM
is not continuous across elements, and the general FEAP
⟨(AM−I)b,Aw⟩=0, (22) commonlycalculatesthecontinuousstressbythemethodof
nodalstressweightedaverage:
in which ⟨,⟩ denotes the inner product. With this relation-
ship,wecanhave
‖(AM−I)b+Aw‖=‖(AM−I)b‖+‖Aw‖
(23) 𝑤 𝜎 +𝑤 𝜎 +⋅⋅⋅+𝑤 𝜎
≥‖(AM−I)b‖, 𝜎̃𝑥 = 1 𝑥(1𝑤 +2𝑤𝑥2+⋅⋅⋅+𝑤 𝑚) 𝑥𝑚, (27)
1 2 𝑚
or
‖AMb−b‖≤󵄩󵄩󵄩󵄩Ag−b󵄩󵄩󵄩󵄩, (24)
whichmeansthatMoore-Penrosepseudoinverseprovidesthe
in which 𝑚 isthenumberofsurroundingelementsforone
linearleastsquaressolution.
Theleastsquaressolutionisusuallynotunique,andthe node, 𝑤𝑖 (𝑖 = 1,2,...,𝑚) is the weight (element area or
generalsolutioncanbeexpressedas volume),and 𝜎𝑥𝑖 (𝑖 = 1,2,...,𝑚)isthestresscalculatedby
the node displacements of each element. The stress field 𝜎̃𝑥
ĝ=M̂b+(I−M̂A)z, ∀z∈R𝑛0×1, M̂∈R𝑛0×𝑛𝑠 (25) evaluated by (27) has the order of accuracy 𝑂(ℎ) at least.
Therefore, the order of the minimum norm least squares
inwhichAM̂A=Aand(AM̂)𝑇 =AM̂. solution of (14) is of 𝑂(ℎ) at least. Moreover, it should not
Similarly,itcanbeverifiedthatifM̂AM̂ = M̂,(M̂A)𝑇 = be ignored that the least squares processing has the strong
̂
MA,then abilitytoimprovetheaccuracy.Thefollowingexampleshows
󵄩󵄩󵄩󵄩󵄩M̂b󵄩󵄩󵄩󵄩󵄩≤󵄩󵄩󵄩󵄩󵄩M̂b+(I−M̂A)z󵄩󵄩󵄩󵄩󵄩. (26) thattheaccuracyofthesuggestedmethod,evenemploying
thelower-orderlinearelements,isequaltoorhigherthanthe
Itindicatesthissolutionhastheminimumnorm. SCPrecoveryusingquadraticelements.
MathematicalProblemsinEngineering 9
areshowninFigure4,includingthenumbersofnodesand
elements.
Start
Firstly, a checking point 𝑃 around (2.2𝑟0,2.2𝑟0) is used
to compare the results of different methods as shown
in Figure4. The exact locations are (2.212𝑟0,2.212𝑟0),
(2.225𝑟0,2.225𝑟0),and(2.232𝑟0,2.232𝑟0)forcoarse,fine,and
ANSYS APDL files input
double-finemeshes,respectively.Resultsfromtheproposed
(Including model, mesh, load, and solve)
stress field gradient (SFG) analysis technique and the SCP
recovery method with various element types and meshes
are compared in Table1. The average error is the root sum
square of 𝑥-direction gradient (𝑆𝑥𝑥) error and 𝑦-direction
Output gradient (𝑆𝑥𝑦) error. The accuracy of these methods is all
improvedwiththemeshrefinementnomatterwhichelement
(1) Nodes location and nodal Von-Mises stress
is employed, especially for quadrilateral elements. Except
(2) Elements information (types and nodes)
fine mesh T6 element (F-T6) only, the presented method
is always more accurate than the SCP recovery technique.
For lower-order linear elements Q4, it is clear that the SFG
technique can still provide a better result than the SCP
Stress field gradient analysis (Matlab)
methodinallthemeshtypes.Ontheotherhand,usingT3
(1) Import ANSYS output files element, results of the SFG get closer to those of the SCP
(2) Confirm nodes and elements information withthemeshrefinement,andfinallyahigheraccuratestress
(3) Find the surrounding elements for each node gradient can be obtained in the double-fine mesh. In the
(4) Identify the sample nodes for each node stress field gradient (SFG) analysis, higher-order elements
do not always correspond to higher accuracy because of
(5) Solve the stress gradient equation
the treatment of minimum norm least squares solution. In
(6) Output the nodal location and stress gradient
otherwords,alower-orderelement,whichgenerallymeans
lowcomputationcost,canprovidealmostthesameaccuracy
as a higher-order element. This is particularly important to
practicalapplications.
Contour figures (Tecplot)
Furthermore, results from the two methods along two
(import the Matlab output files and plot) checking lines, 𝜃 = 45∘ and 𝜃 = 90∘, are shown in Figures
5and6,respectively.Bothfiguresindicatethattheyhavethe
uniformaccuracyinmostareas,includingtheSFGprocedure
employinglinearT3andQ4elements.Althoughtheresults
attheedgeoftheholehavethelowestprecisionforallcases,
End
the accuracy is very high once not at the edge. Specifically,
at the edge of the hole, the SFG method using linear 𝐶0
Figure 8: The flowchart of postprocessing using the presented elementsprovidestheworststressgradientforline𝜃 = 45∘
method. (Line-45∘) as well as the SCP method using quadratic Q8
element for line 𝜃 = 90∘ (Line-𝑌). Utilizing quadrilateral
elements,thetwoapproacheshavealmostidenticalstability
4.NumericalExampleandDiscussions foreachmesh.Ontheotherhand,theSCPrecoverymethod
with triangular elements works better in the stability and
4.1. Example 1: 2D Problem. For existing exact solution convergencecomparingwiththeSFGmethod.
andexhibitingstrongstressconcentrationphenomenon,the
infiniteplatewithacenterholeunderunidirectionaltensile
stress as shown in Figure3 becomes a classic and popular
two-dimensionalexampleusedinpreviouspapers[13–15,24, 4.2.Example2:3DProblem. Athree-dimensionalcantilever
29]. beam with a groove under vertical end load is modeled
In this case, the normal stress 𝜎𝑥 is dominant and the using8-nodehexahedralelements,anditsgeometry,loading
exactsolutionisgivenas condition,andmeshareshowninFigure7.AllDOFsonthe
𝑥=0surfacearefixed.
𝑟2 3 3𝑟4
𝜎𝑥 =𝑃0[1− 𝑟02 (2cos2𝜃+cos4𝜃)+ 2𝑟04 cos4𝜃], (28) of pBoasstepdroocensstihnegswoiftthwatrheepplarotfpoormsedAmNeStYhSo,dthise sflhoowwcnhairnt
Figure8. The code of the stress field gradient analysis is
where𝑟and𝜃arethepolarcoordinates. developed by Matlab language and the contour figures are
Due to the symmetry in the geometry, only a quarter generatedusingsoftwareTecplot.Thegeneralfiniteelement
areaismodeledusingquadrilateralandtriangularelements, analysis program (FEAP), which provides the preliminary
suchasT3,Q4,T6,andQ8elements.Threedifferentmeshes analysisofthestressfield,isnotjustlimitedtoANSYS,and
10 MathematicalProblemsinEngineering
y y y
z x z x z x
(a) (b) (c)
y
y y
z x
z x z x
(d) (e) (f)
Figure9:Von-Misesstressanditsgradientsof3Dproblem:(a)Von-Misesstress,(b)Von-Misesstressgradient-total,(c)stressgradient
slices-total,(d)stressgradient-𝑋,(e)stressgradient-𝑌,and(f)stressgradient-𝑍.
Table 2: Maximal Von-Mises stresses and gradients of different thewidth,andthemaximumofthestressgradientisreduced
groovewidthsanddepths. from 4.15𝐸 + 10 to 2.52𝐸 + 10. The distribution of stress
gradientbecomesmoreandmoreuniformandsmoothwith
ℎ/mm 𝑑/mm Max-Von-Mises/Pa Max Totalgradient
thewidthincreasing,whichisbeneficialtopreventstructural
5 5 1.4339E+08 4.1542E+10
failure.Asthedepthincreases,thetotalstressgradientgoes
10 5 1.2381E+08 2.4833E+10 upfrom4.0531𝐸+10to2.3719𝐸+11andthehigh-gradient
20 5 1.1546E+08 2.5172E+10 zone in the bottom of the groove becomes more and more
4 4 1.3163E+08 4.0531E+10 concentrated.MoredetailsareshowninFigure11,inwhich
4 8 2.3143E+08 9.2233E+10 themiddlelinescrossthebottomsurfacesofthegrooveswith
4 12 4.6672E+08 2.3719E+11 differentwidthsanddepthsalongthe𝑥-direction.
5.Conclusions
the proposed method can also be applied to other software
An FEAP-based mathematical technique is developed for
platformsasPATRANandABUQUS,andsoforth.
accurately extracting stress gradient. Comparing with the
Details of the stress gradients are shown in Figure9 for
SCP recovery method, which needs the quadratic elements
ℎ = 4mm and 𝑑 = 5mm, in which the total gradient
atleastandmustinverttheJacobiandHessianmatrices,this
𝐺(𝜎̃) =√2(𝜕𝜎̃/𝜕𝑥)2+(𝜕𝜎̃/𝜕𝑦)2+(𝜕𝜎̃/𝜕𝑧)2.Itcanbeseen methodonlyrequiresnodalstressresultsaswellaslocation
Total
that the stress gradient in Figure9(b) is more powerful in informationandcanbeimplementedtoanyelementtypes.
characterizingandidentifyingthestressconcentrationthan The classic plane example shows that the suggested
thestressitselfinFigure9(a). Thegradientslicesshowthat method can achieve more accurate stress gradient than the
thestressgradientmaximumappearsaroundthesurfacesand SCPrecoverytechniquewiththesameelements.Inaddition,
corners of the groove, which are the sensitive or dangerous its performance in calculating the stress gradient with the
areastothefatiguefailureandtheplasticdeformation. linear 𝐶0 elements is proved better than the SCP recovery
Resultsofdifferentgroovewidthsanddepthsobtainedby method.Itisalsoobservedthatwiththismethodthequadri-
theSFGtechniquearecomparedinFigure10andTable2.The lateral elements can provide better stability than the trian-
high-gradientzonemovestothefixedendwiththeincreaseof gularelements.Besides,basedontheproposedtechnique,a