研究表明, 汽车质量每减轻10%, 汽车油耗和排放可以降低6%~8%[1].在节能减排和汽车轻量化的背景下, DP590高强钢因比强度高等优点在汽车车身中的使用越来越广泛.高强钢成形过程中的破裂是阻碍其在汽车车身中推广应用的主要问题之一, 因此研究韧性断裂准则在DP590高强钢破裂预测中的应用具有重要的意义.韧性断裂准则是一种建立在韧性损伤力学基础上, 以微观孔洞的成核、生长、联合引起的累积损伤为判据的破裂预测方法.韧性断裂准则主要分为两大类:耦合韧性断裂准则和非耦合韧性断裂准则.耦合断裂准则考虑了损伤对材料带来的软化作用, 而非耦合断裂准则没有考虑该因素.最著名的耦合断裂准则是由Gurson[2]、Tvergaard等[3]提出的GTN模型.耦合韧性断裂准则形式通常比较复杂, 而非耦合韧性断裂准则形式简洁, 计算简单, 在工业生产中应用更加广泛.在前人的不断探索中, 发展出了很多成熟的非耦合韧性断裂准则. Cockcroft等[4]实验观察到韧性断裂往往发生在最大拉应力区域, 提出了一种基于最大主应力的非耦合韧性断裂准则.随后, Oh、Clift和Ko等[5-7]分别在Cockcroft-Latham准则的基础上提出了一些新的修正模型. Rice等[8]通过应力三轴度的指数函数来近似描述单个球形孔洞的生长, 由于简化了分析结果, 使Rice-Tracey模型在板材成形中得到广泛的应用.随后, Oyane等[9]根据多孔材料的塑性理论方程推导出了新的韧性断裂准则.考虑到孔洞的成核、生长和聚结效应, Lou等[10]提出了应力三轴度截断值为-1/3的DF2012准则.为了描述金属在低应力三轴度下的韧性断裂, Cao等[11]基于Lemaitre模型提出了修正的断裂准则; 而Mohr等[12]提出了唯象Hosford-Coulomb模型.近年来, 为了适应新材料和新成形工艺的应用, 出现了一些新的非耦合韧性断裂准则[13-15].
随着韧性断裂准则理论的发展, 其已经被广泛应用于金属板材破裂的预测. Lou等[16]使用DF2012韧性断裂准则成功预测了DP980高强钢的破裂.以DP780高强钢为例, Wang等[17]比较了Brozzo、McClintock、Rice-Tracey和Oyane等4种韧性断裂准则对高强钢破裂预测的精度.陈劼实等[18]使用Clift等4种韧性断裂准则分别预测了HS钢、IF钢以及6111-T4铝合金的成形极限.王在林等[19]将常用的6种韧性断裂准则用于超高强钢辊弯成形破裂预测, 并验证了Brozzo韧性断裂准则在超高强钢破裂预测中的准确性.此外, 余心宏等[20]使用Oyane韧性断裂准则成功预测了铝合金A5182-O和SPCC的成形极限.
随着DP590高强钢在工业中应用越来越广泛, 探索一种更加灵活、准确的韧性断裂准则在其成形模拟中的应用具有重要的研究价值.本文将以DP590高强钢为研究对象, 设计6种单向加载试样获取材料的断裂参数, 运用实验-模拟混合法标定一种灵活的DF2012韧性断裂准则并嵌入ABAQUS软件中进行断裂预测模拟.
1 实验研究 1.1 实验材料和试样制备本研究以厚度为1.5 mm的DP590高强钢板材为研究对象.为了探究DP590高强钢的塑性和断裂性能, 设计了6种不同形状和尺寸的单向加载试样, 如图 1所示.其中, 图 1中试样Ⅰ是狗骨试样, 用于测定DP590的塑性性能; 其余分别是圆孔试样、R5缺口试样、R10缺口试样、R20缺口试样和剪切试样, 这些试样包含纯剪切到平面应变的应力状态, 用于测定DP590的材料断裂参数.
为研究DP590高强钢板的塑性各向异性, 每种试样均制备了与轧制方向成0°、45°和90°的3个试样.
1.2 实验过程本研究中, 6种试样的拉伸实验均可在同一台岛津AG-X Plus 100 kN万能试验机上完成, 拉伸速度设置为2 mm/min, 数据采样频率为100 Hz.为避免实验数据的偶然性, 每组实验均重复进行3次并取平均值.实验中材料的变形过程通过高清单反相机进行记录, 图像采样频率为10 Hz, 最后使用数字图像相关(DIC)软件Ncorr[21]计算试样的全场应变.以狗骨试样为例, 通过实验得到3个取向的工程应力-工程应变曲线如图 2所示.
本研究标定断裂准则参数采用的是一种实验-模拟混合方法, 选用由Lou等[22]提出的各向异性Drucker屈服函数, 该屈服函数已被验证可以准确地表征面心立方和体心立方晶体结构金属的各向异性塑性, 并且该屈服函数相对其他经典屈服函数可以大大提高仿真模拟的计算效率.各向异性Drucker屈服函数的基本形式为
$ \overline{\boldsymbol{\sigma}}=\left(J_{2}^{\prime 3}-c J_{3}^{\prime 2}\right)^{1 / 6}. $ |
式中:参数c是经验常数, 这里c取建议值1.226[22], J′2和J′3是应力张量s′的第2和第3应力不变量, 这两个应力不变量的计算表达式如下:
$ J_{2}^{\prime}=(1 / 2) \boldsymbol{s}^{\prime}: \boldsymbol{s}^{\prime}=-s_{11}^{\prime} s_{22}^{\prime}-s_{22}^{\prime} s_{33}^{\prime}-s_{11}^{\prime} s_{33}^{\prime}+\\s_{12}^{\prime 2}+s_{23}^{\prime 2}+s_{13}^{\prime 2}, $ | (1) |
$ J_{3}^{\prime}=\operatorname{det}\left(s^{\prime}\right)=-s_{11}^{\prime} s_{22}^{\prime} s_{33}^{\prime}+2 s_{12}^{\prime} s_{23}^{\prime} s_{13}^{\prime}-\\s_{11}^{\prime} s_{23}^{\prime 2}-s_{22}^{\prime} s_{13}^{\prime 2}-s_{33}^{\prime} s_{12}^{\prime 2}. $ | (2) |
式(1)、(2)中应力张量s′的应力分量计算如下:
$ {\mathit{\boldsymbol{s}}^\prime } \equiv \left[ {\begin{array}{*{20}{c}} {s_{11}^\prime }\\ {s_{22}^\prime }\\ {s_{33}^\prime }\\ {s_{23}^\prime }\\ {s_{13}^\prime }\\ {s_{12}^\prime } \end{array}} \right] = {\mathit{\boldsymbol{L}}^\prime }\left[ {\begin{array}{*{20}{l}} {{\mathit{\sigma }_{xx}}}\\ {{\mathit{\sigma }_{yy}}}\\ {{\mathit{\sigma }_{zz}}}\\ {{\mathit{\sigma }_{yz}}}\\ {{\mathit{\sigma }_{xz}}}\\ {{\mathit{\sigma }_{xy}}} \end{array}} \right], $ |
$ \boldsymbol{L}^{\prime}=\left[\begin{array}{cccccc} \frac{c_{2}^{\prime}+c_{3}^{\prime}}{3} & \frac{-c_{3}^{\prime}}{3} & \frac{-c_{2}^{\prime}}{3} & 0 & 0 & 0 \\ \frac{-c_{3}^{\prime}}{3} & \frac{c_{3}^{\prime}+c_{1}^{\prime}}{3} & \frac{-c_{1}^{\prime}}{3} & 0 & 0 & 0 \\ \frac{-c_{2}^{\prime}}{3} & \frac{-c_{1}^{\prime}}{3} & \frac{c_{1}^{\prime}+c_{2}^{\prime}}{3} & 0 & 0 & 0 \\ 0 & 0 & 0 & c_{4}^{\prime} & 0 & 0 \\ 0 & 0 & 0 & 0 & c_{5}^{\prime} & 0 \\ 0 & 0 & 0 & 0 & 0 & c_{6}^{\prime} \end{array}\right]. $ |
在各向异性Drucker屈服函数中, 共有6个各向异性参数.其中, c′1、c′2、c′3、c′6是面内各向异性塑性参数, 可用0°、45°和90°方向的拉伸屈服应力以及等双拉的屈服应力σ0、σ45、σ90、σb进行标定.因为厚度方向的材料属性难以通过实验测得, 所以厚向塑性参数c′4、c′5通常取值与c′6相等[22].
表 1为实验测定的DP590沿轧制方向各个角度的屈服应力和厚向异性系数r值.由于实验条件的限制, 本研究未进行等双向拉伸实验, 而采用一种等效的方法获得σb和rb.
结合表 1和图 2可知, DP590初始屈服应力的各向异性不明显, 因此等双拉屈服应力σb使用理论公式σb=(σ0+2σ45+σ90)/4进行求解, 屈服应力各向异性不明显时, 采用该方法可以同时兼顾参数求解的效率和理论预测值的精度; 因为厚向异性系数r值的各向异性较为明显, 采用各向异性屈服函数Yld1996对rb进行预测, 尽可能使得rb预测值更加准确.使用屈服应力标定各向异性Drucker屈服函数, 参数标定结果如表 2中Drucker_y对应数据所示.为了增加各向异性Drucker屈服准则的灵活性, 数值模拟中使用了非相关塑性流动法则, 材料的塑性势函数形式与屈服面的函数形式一致, 改用厚向异性系数r值对塑性面的参数进行标定, 标定结果如表 2中Drucker_p对应数据所示.
本文还比较了Hill48屈服函数与各向异性Drucker屈服函数对DP590塑性各向异性的预测精度. Hill48屈服函数基本形式如式(3)所示:
$ F\left(\sigma_{y y}-\sigma_{z z}\right)^{2}+G\left(\sigma_{z z}-\sigma_{x x}\right)^{2}+H\left(\sigma_{x x}-\sigma_{y y}\right)^{2}+\\2 L \sigma_{y z}^{2}+2 M \sigma_{z x}^{2}+2 N \sigma_{x y}^{2}=2 \bar{\sigma}^{2}. $ | (3) |
使用厚向异性系数r值对其参数进行标定, 结果如表 2所示.
根据表 2中屈服函数的参数, 绘制屈服函数的面内屈服轨迹如图 3所示.由图 3可知, 各向异性Drucker屈服函数预测的屈服轨迹与实验结果十分吻合; 而Hill48屈服准则未能准确预测除轧制方向外的其他两个方向的屈服应力.
试样变形过程中, 断裂应变通常显著高于最大均匀应变.为了表征DP590的硬化行为, 需要建立合适的硬化模型来准确预测大应变下的真实应力.常用的Swift非饱和型硬化模型和Voce饱和型硬化模型均不能很好地预测DP590在大应变下的应力, 因此研究中使用了修正的Voce硬化准则(指数型Voce+Voce硬化准则)表征DP590的硬化行为, 其硬化方程为
$ \sigma(\varepsilon)=\sigma_{\mathrm{Y}}+A_{1}\left(1-\mathrm{e}^{-\beta_{1} \varepsilon}\right)+A_{2}\left(-\mathrm{e}^{-\beta_{2} \varepsilon}\right). $ |
式中:ε是等效塑性应变, σY是初始屈服应力, A1、β1、A2、β2是4个硬化常数.
使用轧制方向的狗骨试样实验数据拟合硬化参数, 用最小二乘法求得:σY=347.53 MPa, A1=313.17 MPa, β1=18.66, A2=377.36 MPa, β2=1.84.
为验证硬化准则的准确性, 将其用于剪切试样的仿真模拟, 并与Swift、Voce硬化准则进行对比, 模拟和实验结果如图 4所示.由图 4可知, Voce+Voce硬化准则可以准确表征DP590高强钢的硬化行为, 而Swift和Voce硬化方程分别高估和低估了DP590在大应变下的应力.
为进一步验证各向异性Drucker屈服函数和Voce+Voce硬化方程对DP590高强钢塑性行为预测的准确性, 将该本构模型通过Fortran语言编译为VUMAT子程序后嵌入ABAQUS有限元软件, 用于模拟不同应力状态下DP590板材的加载过程.模拟选用了R5缺口、圆孔和剪切试样3种试样.
首先, 根据各个试样的几何尺寸建立有限元模型.为提高模拟的计算效率, 利用试样的对称性, 针对缺口试样和圆孔试样均建立了1/8对称模型; 针对剪切试样建立了厚向对称的1/2模型.此外, 为兼顾模拟的精度和计算效率, 划分网格时在塑性变形较大的区域划分较细的网格, 并逐渐增大过渡到弹性变形为主的区域.经过网格敏感性测试后, 各试样大变形区域均划分大小约为0.05 mm的网格.
图 5是各试样3个不同取向(0°, 45°, 90°)的实验和模拟载荷-位移曲线对比, 为便于区分, 圆孔试样和剪切试样曲线的横坐标同时向右平移了2 mm.由图 5可知, 各向异性Drucker屈服函数很好地再现了DP590的3种试样不同取向的载荷-位移曲线, 而且表现出了比Hill48屈服函数更好的预测精度.因此, 各向异性Drucker屈服函数与Voce+Voce硬化准则较好地表征了DP590高强钢的塑性行为, 可以用于标定韧性断裂准则.
在众多韧性断裂准则中, 本研究以DF2012韧性断裂准则[10]为理论基础, 这是一种基于应力三轴度、罗德参数和等效塑性应变构建断裂判据的准则, 其表达式如表 3所示.将3种经典韧性断裂准则Cockcroft-Latham (C-L)[4]、Clift[6]和Rice-Tracey (R-T)[8]与DF2012准则进行对比.由于应力状态和应变强度是控制韧性断裂萌生的最重要因素, 所有模型均表示为应力三轴度、罗德参数和等效塑性应变的形式.
本研究选用R10缺口试样、圆孔试样和剪切试样标定断裂准则的参数.由于实验难以直接获得应力三轴度、罗德参数和断裂应变, 因此, 在单拉实验结果和标定的本构模型基础上, 采用实验-模拟混合法标定韧性断裂准则.
以R10缺口试样为例, 阐述实验-模拟混合法确定断裂参数的过程.首先, 根据图 6中试样断裂位置对应表层单元A的实验-模拟应变演化对比以及试样的力程曲线对比来验证模拟模型的准确性.然后, 使用试样的实验断裂位置处对应单元(即中心层单元B)的等效塑性应变模拟结果预测断裂应变.值得注意的是, 本研究将实验中试样承载能力急剧下降的时刻等效为断裂开始时刻, 将断裂单元在该时刻的等效塑性应变作为断裂应变.使用同样的方法, 可以分别获得各试样的断裂应变, 如表 4所示.
除了断裂应变, 参数标定还需要获得试样断裂单元处的应力三轴度和罗德参数. 图 7是通过模拟获得的3种试样的应力三轴度和罗德参数.
分别计算各试样的应力三轴度和罗德参数均值:
$ \eta_{\text {avg }}=\frac{1}{\bar{\varepsilon}_{f}^{p}} \int_{0}^{\bar{\varepsilon}_{f}^{P}} \eta\left(\bar{\varepsilon}^{p}\right) \mathrm{d} \bar{\varepsilon}^{p}, $ | (8) |
$ L_{\text {avg }}=\frac{1}{\bar{\varepsilon}_{f}^{p}} \int_{0}^{\bar{\varepsilon}_{f}^{p}} L\left(\bar{\varepsilon}^{p}\right) \mathrm{d} \bar{\varepsilon}^{p}. $ | (9) |
计算结果如表 4所示.
根据表 4, 用优化方法求解韧性断裂准则参数, 结果见表 5.使用DF2012构建DP590高强钢在应力三轴度和等效塑性应变(η, ε)二维空间下的断裂轨迹, 如图 8所示, 这是DP590的断裂判据. 图 8还比较了C-L、Clift和R-T等3种准则构建的断裂轨迹.由图 8中断裂轨迹与实验数据的对比可知, DF2012韧性断裂准则具有足够的灵活性, 为DP590构建了合理的断裂轨迹.
为了验证构建的断裂轨迹的准确性, 使用Fortran语言将韧性断裂准则编译后嵌入ABAQUS仿真软件, 并作为断裂判据用于DP590高强钢的拉伸仿真模拟.除了狗骨试样, 本研究中其余5种试样均用于断裂模拟中, 其中圆孔、R10缺口和剪切试样参与了断裂准则的参数标定, R5和R20缺口试样并未参与过参数标定.以圆孔和R20试样为例, 对比了实验和模拟的试样断裂形貌, 结果如图 9所示.
由图 9可知, 有限元模拟较好地还原了试样的断裂形式.为进一步验证韧性断裂准则的准确性, 同样以圆孔和R20试样为例, 绘制实验和模拟载荷-位移对比曲线如图 10所示.将所有试样的实验和模拟断裂位移进行汇总和比较, 如图 11所示.
由图 11可知, DF2012韧性断裂准则较为准确地预测了各试样断裂的发生, 比C-L、Clift和R-T准则精度更高. DF2012韧性断裂准则对剪切试样断裂位移的预测精度明显高于其他准则.模拟结果验证了DF2012韧性断裂准则的灵活性和对DP590断裂预测的准确性, 至少在应力三轴度值为(0~0.66)的应力状态范围内得以验证.这是因为DF2012韧性断裂准则同时考虑了应力三轴度和罗德参数在断裂发生中的作用, 并合理构建了应力三轴度、罗德参数和断裂应变之间的关系.
4 结论基于DP590高强钢6种不同试样的单向加载实验, 本文首先研究了DP590的塑性各向异性并确定本构模型, 然后使用实验-模拟混合法标定了DF2012韧性断裂准则, 并用于预测DP590高强钢板材的断裂, 主要得出以下结论:
1) 对比试样沿轧制方向的0°、45°和90°等3个取向的实验和模拟载荷-位移曲线可知, 与Hill48屈服函数相比, 各向异性Drucker屈服函数可以准确地表征DP590的各向异性塑性行为.
2) Voce+Voce、Swift和Voce等3种硬化准则对DP590在大应变下硬化行为的预测结果差异较大.对比剪切试样的实验和模拟载荷-位移曲线可知, Voce+Voce硬化方程相对其他两者, 能正确表征DP590高强钢的硬化行为.
3) 采用实验-模拟混合法可以方便、准确地标定韧性断裂准则参数.由理论断裂轨迹与实验值对比可知, DF2012韧性断裂准则具有足够的灵活性, 可以很好地构建DP590的断裂判据.
4) 将DF2012韧性断裂准则用于DP590的拉伸断裂模拟, 可以准确地预测不同应力状态试样断裂的发生, 验证了DF2012韧性断裂准则对DP590断裂预测具有足够的精度.
[1] |
SUN Wenlong, CHEN Xiaokai, WANG Lu. Analysis of energy saving and emission reduction of vehicles using light weight materials[J]. Energy Procedia, 2016, 88: 891. DOI:10.1016/j.egypro.2016.06.106 |
[2] |
GURSON A L. Continuum theory of ductile rupture by void nucleation and growth: part I: yield criteria and flow rules for porous ductile media[J]. Journal of Engineering Materials and Technology, 1977, 99(1): 4. DOI:10.1115/1.3443401 |
[3] |
TVERGAARD V, NEEDLEMAN A. Analysis of the cup-cone fracture in a round tensile bar[J]. Acta Metallurgica, 1984, 32(1): 159. DOI:10.1016/0001-6160(84)90213-X |
[4] |
COCKROFT M G, LATHAM D J. Ductility and the workability of metals[J]. Journal Institute of Metals, 1968, 96: 33. |
[5] |
OH S I, CHEN C C, KOBAYASHI S. Ductile fracture in axisymmetric extrusion and drawing: part 2: workability in extrusion and drawing[J]. Journal of Engineering for Industry, 1977, 101(1): 36. DOI:10.1115/1.3439471 |
[6] |
CLIFT S E, HARTLEY P, STURGESS C E, et al. Fracture prediction in plastic deformation processes[J]. International Journal of Mechanical Sciences, 1990, 32(1): 3. DOI:10.1016/0020-7403(90)90148-C |
[7] |
KO Y K, LEE J, HUH H, et al. Prediction of fracture in hub-hole expanding process using a new ductile fracture criterion[J]. Journal of Materials Processing Technology, 2007, 187/188: 361. DOI:10.1016/j.jmatprotec.2006.11.071 |
[8] |
RICE J R, TRACEY D M. On the ductile enlargement of voids in triaxial stress fields[J]. Journal of the Mechanics and Physics of Solids, 1969, 17(3): 202. DOI:10.1016/0022-5096(69)90033-7 |
[9] |
OYANE M, SATO T, OKIMOTO K, et al. Criteria for ductile fracture and their applications[J]. Journal of Mechanical Working Technology, 1980, 4(1): 66. DOI:10.1016/0378-3804(80)90006-6 |
[10] |
LOU Yanshan, HUH H, LIM S J, et al. New ductile fracture criterion for prediction of fracture forming limit diagrams of sheet metals[J]. International Journal of Solids and Structures, 2012, 49(25): 3607. DOI:10.1016/j.ijsolstr.2012.02.016 |
[11] |
CAO T S, GACHET J M, MONTMITONNET P, et al. A Lode-dependent enhanced Lemaitre model for ductile fracture prediction at low stress triaxiality[J]. Engineering Fracture Mechanics, 2014, 124/125: 83. DOI:10.1016/j.engfracmech.2014.03.021 |
[12] |
MOHR D, MARCADET S J. Micromechanically-motivated phenomenological Hosford-Coulomb model for predicting ductile fracture initiation at low stress triaxialities[J]. International Journal of Solids and Structures, 2015, 67/68: 43. DOI:10.1016/j.ijsolstr.2015.02.024 |
[13] |
LEE J Y, LEE M G, BARLAT F, et al. Piecewise linear approximation of nonlinear unloading-reloading behaviors using a multi-surface approach[J]. International Journal of Plasticity, 2017, 93: 114. DOI:10.1016/j.ijplas.2017.02.004 |
[14] |
LOU Yanshan, YOON J W. Anisotropic ductile fracture criterion based on linear transformation[J]. International Journal of Plasticity, 2017, 93: 10. DOI:10.1016/j.ijplas.2017.04.008 |
[15] |
CAO Jun, LI Fuguo, MA Xinkai, et al. A modified elliptical fracture criterion to predict fracture forming limit diagrams for sheet metals[J]. Journal of Materials Processing Technology, 2018, 252: 117. DOI:10.1016/j.jmatprotec.2017.09.018 |
[16] |
LOU Yanshan, HUH H. Prediction of ductile fracture for advanced high strength steel with a new criterion: Experiments and simulation[J]. Journal of Materials Processing Technology, 2013, 213(8): 1298. DOI:10.1016/j.jmatprotec.2013.03.001 |
[17] |
WANG Changsheng, CHEN Jun, XIA C, et al. A new method to calculate threshold values of ductile fracture criteria for advanced high-strength sheet blanking[J]. Journal of Materials Engineering and Performance, 2014, 23(4): 1304. DOI:10.1007/s11665-013-0861-z |
[18] |
陈劼实, 周贤宾. 成形极限预测韧性断裂准则及屈服准则的影响[J]. 北京航空航天大学学报, 2006, 32(8): 972. CHEN Jieshi, ZHOU Xianbin. The influence of ductile fracture criterion and yield criterion on forming limit prediction[J]. Journal of Beijing University of Aeronautics and Astronautics, 2006, 32(8): 972. DOI:10.3969/j.issn.1001-5965.2006.08.023 |
[19] |
王在林, 韩飞, 刘继英, 等. 韧性断裂准则在超高强钢辊弯成形工艺中的应用[J]. 塑性工程学报, 2012, 19(4): 19. WANG Zailin, HAN Fei, LIU Jiying, et al. Application of ductile fracture criterion in roll forming process of ultra-high strength steel[J]. Journal of Plastic Engineering, 2012, 19(4): 19. DOI:10.3969/j.issn.1007-2012.2012.04.004 |
[20] |
余心宏, 翟妮芝, 翟江波. 应用韧性断裂准则预测板料的成形极限图[J]. 锻压技术, 2007, 32(5): 46. YU Xinhong, ZHAI Nizhi, ZHAI Jiangbo. The forming limit diagram of sheet metal is predicted by ductile fracture criterion[J]. Forging & Stamping Technology, 2007, 32(5): 46. DOI:10.3969/j.issn.1000-3940.2007.05.012 |
[21] |
BLABER J A, ADAIR B, ANTONIOU A, et al. Ncorr: Open-source 2D Digital Image Correlation Matlab software[J]. Experimental Mec-hanics, 2015, 55(6): 1106. DOI:10.1007/s11340-015-0009-1 |
[22] |
LOU Yanshan, YOON J W. Anisotropic yield function based on stress invariants for BCC and FCC metals and its extension to ductile fracture criterion[J]. International Journal of Plasticity, 2018, 101: 133. DOI:10.1016/j.ijplas.2017.10.012 |