Parton-Shower Simulations of R-parity Violating Supersymmetric Models

A Facial Aging Simulation Method Using flaccidity deformation criteriaAlexandre Cruz Berg Lutheran University of Brazil.Dept Computer ScienceRua Miguel Tostes, 101. 92420-280 Canoas, RS, Brazil berg@ulbra.tche.br Francisco José Perales LopezUniversitat les Illes Balears.Dept Mathmatics InformaticsCtra Valldemossa, km 7,5E-07071 Palma MallorcaSpainpaco.perales@uib.esManuel GonzálezUniversitat les Illes Balears.Dept Mathmatics InformaticsCtra Valldemossa, km 7,5E-07071 Palma MallorcaSpainmanuel.gonzales@uib.esAbstractDue to the fact that the aging human face encompasses skull bones, facial muscles, and tissues, we render it using the effects of flaccidity through the observation of family groups categorized by sex, race and age. Considering that patterns of aging are consistent, facial ptosis becomes manifest toward the end of the fourth decade. In order to simulate facial aging according to these patterns, we used surfaces with control points so that it was possible to represent the effect of aging through flaccidity. The main use of these surfaces is to simulate flaccidity and aging consequently.1.IntroductionThe synthesis of realistic virtual views remains one of the central research topics in computer graphics. The range of applications encompasses many fields, including: visual interfaces for communications, integrated environments of virtual reality, as well as visual effects commonly used in film production.The ultimate goal of the research on realistic rendering is to display a scene on a screen so that it appears as if the object exists behind the screen. This description, however, is somewhat ambiguous and doesn't provide a quality measure for synthesized images. Certain areas, such as plastic surgery, need this quality evaluation on synthesized faces to make sure how the patient look like and more often how the patient will look like in the future. Instead, in computer graphics and computer vision communities, considerable effort has been put forthto synthesize the virtual view of real or imaginary scenes so that they look like the real scenes.Much work that plastic surgeons put in this fieldis to retard aging process but aging is an inevitable process. Age changes cause major variations in the appearance of human faces [1]. Some aspects of aging are uncontrollable and are based on hereditary factors; others are somewhat controllable, resulting from many social factors including lifestyle, among others [2].1.1.Related WorkMany works about aging human faces have been done. We can list some related work in the simulation of facial skin deformation [3].One approach is based on geometric models, physically based models and biomechanical models using either a particle system or a continuous system.Many geometrical models have been developed, such as parametric model [4] and geometric operators [5]. The finite element method is also employed for more accurate calculation of skin deformation, especially for potential medical applications such as plastic surgery [6]. Overall, those works simulate wrinkles but none of them have used flaccidity as causing creases and aging consequently.In this work is presented this effort in aging virtual human faces, by addressing the synthesis of new facial images of subjects for a given target age.We present a scheme that uses aging function to perform this synthesis thru flaccidity. This scheme enforces perceptually realistic images by preserving the identity of the subject. The main difference between our model and the previous ones is that we simulate increase of fat and muscular mass diminish causing flaccidity as one responsible element for the sprouting of lines and aging human face.In the next section will plan to present the methodology. Also in section 3, we introduce the measurements procedure, defining structural alterations of the face. In section 4, we present a visual facial model. We describe age simulation thrua deformation approach in section 5. In the last section we conclude the main results and future work.2.MethodologyA methodology to model the aging of human face allows us to recover the face aging process. This methodology consists of: 1) defining the variations of certain face regions, where the aging process is perceptible; 2) measuring the variations of those regions for a period of time in a group of people and finally 3) making up a model through the measurements based on personal features.That could be used as a standard to a whole group in order to design aging curves to the facial regions defined.¦njjjpVM2.1Mathematical Background and AnalysisHuman society values beauty and youth. It is well known that the aging process is influenced by several parameters such: feeding, weight, stress level, race, religious factors, genetics, etc. Finding a standard set of characteristics that could possibly emulate and represent the aging process is a difficult proposition.This standard set was obtained through a mathematical analysis of some face measurements in a specific group of people, whose photographs in different ages were available [7]. To each person in the group, there were, at least, four digitized photographs. The oldest of them was taken as a standard to the most recent one. Hence, some face alterations were attained through the passing of time for the same person.The diversity of the generated data has led to the designing of a mathematical model, which enabled the acquiring of a behavior pattern to all persons of the same group, as the form of a curve defined over the domain [0,1] in general, in order to define over any interval [0,Į] for an individual face. The unknown points Įi are found using the blossoming principle [8] to form the control polygon of that face.The first step consisted in the selection of the group to be studied. Proposing the assessment of the face aging characteristics it will be necessary to have a photographic follow-up along time for a group of people, in which their face alterations were measurable.The database used in this work consisted of files of patients who were submitted to plastic surgery at Medical Center Praia do Guaíba, located in Porto Alegre, Brazil.3.MeasurementsAccording to anatomic principles [9] the vectors of aging can be described aswhich alter the position and appearance of key anatomic structures of the face as can be shown in figure 1 which compares a Caucasian mother age 66 (left side) with her Caucasian daughters, ages 37 (right above) and 33 (right below) respectively.Figure 1 - Observation of family groupsTherefore, basic anatomic and surgical principles must be applied when planning rejuvenative facial surgery and treating specific problems concomitantwith the aging process.4.Visual Facial ModelThe fact that human face has an especially irregular format and interior components (bones, muscles and fabrics) to possess a complex structure and deformations of different face characteristics of person to person, becomes the modeling of the face a difficult task. The modeling carried through in the present work was based on the model, where the mesh of polygons corresponds to an elastic mesh, simulating the dermis of the face. The deformations in this mesh, necessary to simulate the aging curves, are obtained through the displacement of the vertexes, considering x(t) as a planar curve, which is located within the (u,v ) unit square. So, we can cover the square with a regular grid of points b i,j =[i/m,j/n]T ; i=0,...,m; j=0,...,n. leading to every point (u,v ) asfrom the linear precision property of Bernstein polynomials. Using comparisons with parents we can distort the grid of b i,j into a grid b'i,j , the point (u,v )will be mapped to a point (u',v') asIn order to construct our 3D mesh we introduce the patch byAs the displacements of the vertexes conform to the certain measures gotten through curves of aging and no type of movement in the face is carried through, the parameters of this modeling had been based on the conformation parameter.4.1Textures mappingIn most cases the result gotten in the modeling of the face becomes a little artificial. Using textures mapping can solve this problem. This technique allows an extraordinary increase in the realism of the shaped images and consists of applying on the shaped object, existing textures of the real images of the object.In this case, to do the mapping of an extracted texture of a real image, it is necessary that the textureaccurately correspond to the model 3D of that is made use [9].The detected feature points are used for automatic texture mapping. The main idea of texture mapping is that we get an image by combining two orthogonal pictures in a proper way and then give correct texture coordinates of every point on a head.To give a proper coordinate on a combined image for every point on a head, we first project an individualized 3D head onto three planes, the front (x, y), the left (y, z) and the right (y, z) planes. With the information of feature lines, which are used for image merging, we decide on which plane a 3D-head point on is projected.The projected points on one of three planes arethen transferred to one of feature points spaces suchas the front and the side in 2D. Then they are transferred to the image space and finally to the combined image space.The result of the texture mapping (figure 2) is excellent when it is desired to simulate some alteration of the face that does not involve a type of expression, as neutral. The picture pose must be the same that the 3D scanned data.¦¦¦ mi nj lk n j m i lk k j i w B v B u B b w v u 000,,)()()(')',','(¦¦ m i nj n jmij i v B u B b v u 00,)()(),(¦¦ m i nj n j m i j i v B u B b v u 00,)()(')','(¦¦¦ mi nj lk n j m i lk k j i w B v B u B b w v u 000,,)()()(')',','(Figure 2 - Image shaped with texturemapping5.Age SimulationThis method involves the deformation of a face starting with control segments that define the edges of the faces, as¦¦¦ mi nj lk n j m i lk k j i w B v B u B b w v u 000,,)()()(')',','(Those segments are defined in the original face and their positions are changed to a target face. From those new positions the new position of each vertex in the face is determined.The definition of edges in the face is a fundamental step, since in that phase the applied aging curves are selected. Hence, the face is divided in influencing regions according to their principal edges and characteristics.Considering the face morphology and the modeling of the face aging developed [10], the face was divided in six basic regions (figure 3).The frontal region (1) is limited by the eyelids and the forehead control lines. The distance between these limits enlarges with forward aging.The orbitary region (2) is one of the most important aging parameters because a great number of wrinkles appears and the palpebral pouch increases [11]. In nasal region (3) is observed an enlargement of its contour.The orolabial region (4) is defined by 2 horizontal control segments bounding the upper and lower lips and other 2 segments that define the nasogenian fold. Figure 3 - Regions considering the agingparametersThe lips become thinner and the nasogenian fold deeper and larger. The mental region (5) have 8 control segments that define the low limit of the face and descend with aging. In ear curve (6) is observed an enlargement of its size. The choice of feature lines was based in the characteristic age points in figure 6.The target face is obtained from the aging curves applied to the source face, i.e., with the new control segment position, each vertex of the new image has its position defined by the corresponding vertex in the target face. This final face corresponds to the face in the new age, which was obtained through the application of the numerical modeling of the frontal face aging.The definition of the straight-line segment will control the aging process, leading to a series of tests until the visual result was adequate to the results obtained from the aging curves. The extremes of the segments are interpolated according to the previously defined curves, obtained by piecewise bilinear interpolation [12].Horizontal and vertical orienting auxiliary lines were defined to characterize the extreme points of the control segments (figure 4). Some points, that delimit the control segments, are marked from the intersection of the auxiliary lines with the contour of the face, eyebrow, superior part of the head and the eyes. Others are directly defined without the use of auxiliary lines, such as: eyelid hollow, eyebrow edges, subnasion, mouth, nasolabial wrinkle andnose sides.Figure 4 - Points of the control segmentsOnce the control segments characterize the target image, the following step of the aging process can be undertaken, corresponding to the transformations of the original points to the new positions in the target image. The transformations applied to the segments are given by the aging curves, presented in section 4.In the present work the target segments are calculated by polynomial interpolations, based on parametric curves [12].5.1Deformation approachThe common goal of deformation models is to regulate deformations of a geometric model by providing smoothness constraints. In our age simulation approach, a mesh-independent deformation model is proposed. First, connected piece-wise 3D parametric volumes are generated automatically from a given face mesh according to facial feature points.These volumes cover most regions of a face that can be deformed. Then, by moving the control pointsof each volume, face mesh is deformed. By using non-parallel volumes [13], irregular 3D manifolds are formed. As a result, smaller number of deformvolumes are necessary and the number of freedom incontrol points are reduced. Moreover, based on facialfeature points, this model is mesh independent,which means that it can be easily adopted to deformany face model.After this mesh is constructed, for each vertex on the mesh, it needs to be determined which particularparametric volume it belongs to and what valueparameters are. Then, moving control points ofparametric volumes in 3D will cause smooth facialdeformations, generating facial aging throughflaccidity, automatically through the use of the agingparameters. This deformation is written in matricesas , where V is the nodal displacements offace mesh, B is the mapping matrix composed ofBernstein polynomials, and E is the displacementvector of parametric volume control nodes.BE V Given a quadrilateral mesh of points m i,j ,, we define acontinuous aged surface via a parametricinterpolation of the discretely sampled similaritiespoints. The aged position is defined via abicubic polynomial interpolation of the form with d m,n chosen to satisfy the known normal and continuity conditions at the sample points x i,j .>@>M N j i ,...,1,...,1),(u @@>@>1,,1,),,( j j v i i u v u x ¦3,,),(n m n m n m v u d v u x An interactive tool is programmed to manipulate control points E to achieve aged expressions making possible to simulate aging through age ranges. Basic aged expression units are orbicularis oculi, cheek, eyebrow, eyelid, region of chin, and neck [14]. In general, for each segment, there is an associated transformation, whose behavior can be observed by curves. The only segments that do not suffer any transformation are the contour of the eyes and the superior side of the head.5.2Deformation approachThe developed program also performs shape transformations according to the created aging curves, not including any quantification over the alterations made in texture and skin and hair color. Firstly, in the input model the subjects are required to perform different ages, as previouslymentioned, the first frame needs to be approximately frontal view and with no expression.Secondly, in the facial model initialization, from the first frame, facial features points are extracted manually. The 3D fitting algorithm [15] is then applied to warp the generic model for the person whose face is used. The warping process and from facial feature points and their norms, parametric volumes are automatically generated.Finally, aging field works to relieve the drifting problem in template matching algorithm, templates from the previous frame and templates from the initial frame are applied in order to combine the aging sequence. Our experiments show that this approach is very effective. Despite interest has been put in presenting a friendly user interface, we have to keep in mind that the software system is research oriented. In this kind of applications an important point is the flexibility to add and remove test facilities. 6.Results The presented results in the following figuresrefer to the emulations made on the frontalphotographs, principal focus of this paper, with theobjective to apply the developed program to otherpersons outside the analyzed group. The comparisonswith other photographs of the tested persons dependon their quality and on the position in which theywere taken. An assessment was made of the new positions, of the control segments. It consisted in: after aging a face, from the first age to the second one, through the use of polynomial interpolation of the control segments in the models in the young age, the new positions are then compared with the ones in the model of a relative of older age (figure 5). The processed faces were qualitatively compared with theperson’s photograph at the same age. Figure 5 - Synthetic young age model,region-marked model and aged modelAlso the eyelid hollow, very subtle falling of the eyebrow, thinning of the lips with the enlarging of the nasion and the superior part of the lip, enlargingof the front and changing in the nasolabial wrinkle.7.ConclusionsModelling biological phenomena is a great deal of work, especially when the biggest part of the information about the subject involves only qualitative data. Thus, this research developed had has a challenge in the designing of a model to represent the face aging from qualitative data.Due to its multi-disciplinary character, the developed methodology to model and emulate the face aging involved the study of several other related fields, such as medicine, computing, statistics and mathematics.The possibilities opened by the presented method and some further research on this field can lead to new proposals of enhancing the current techniques of plastic face surgery. It is possible to suggest the ideal age to perform face lifting. Once the most affected aging regions are known and how this process occurs over time. Also missing persons can be recognized based on old photographs using this technique. AcknowledgementsThe project TIN2004-07926 of Spanish Government have subsidized this work.8. References[1] Burt, D. M. et al., Perc. age in adult Caucasianmale faces, in Proc. R. Soc., 259, pp 137-143,1995.[2] Berg, A C. Aging of Orbicularis Muscle inVirtual Human Faces. IEEE 7th InternationalConference on Information Visualization, London, UK, 2003a.[3] Beier , T., S. Neely, Feature-based imagemetamorphosis, In Computer Graphics (Proc.SIGGRAPH), pp. 35-42, 1992.[4] Parke, F. I. P arametrized Models for FacialAnimation, IEEE Computer & Graphics Applications, Nov. 1982.[5] Waters, K.; A Muscle Model for Animating ThreeDimensional Facial Expression. Proc SIGGRAPH'87,Computer Graphics, Vol. 21, Nº4, United States, 1987. [6] Koch, R.M. et alia.. Simulation Facial SurgeryUsing Finite Element Models, Proceedings of SIGGRAPH'96, Computer Graphics, 1996.[7] Kurihara, Tsuneya; Kiyoshi Arai, ATransformation Method for Modeling and Animation of the Human Face from Photographs, Computer Animatio n, Springer-Verlag Tokyo, pp.45-58, 1991.[8] Kent, J., W. Carlson , R. Parent, ShapeTransformation for Polygon Objects, In Computer Graphics (Proc. SIGGRAPH), pp. 47-54, 1992. [9] Sorensen, P., Morphing Magic, in ComputerGraphics World, January 1992.[10]Pitanguy, I., Quintaes, G. de A., Cavalcanti, M.A., Leite, L. A. de S., Anatomia doEnvelhecimento da Face, in Revista Brasileira deCirurgia, Vol 67, 1977.[11]Pitanguy, I., F. R. Leta, D. Pamplona, H. I.Weber, Defining and measuring ageing parameters, in Applied Mathematics and Computation , 1996.[12]Fisher, J.; Lowther, J.; Ching-Kuang S. Curveand Surface Interpolation and Approximation: Knowledge Unit and Software Tool. ITiCSE’04,Leeds, UK June 28–30, 2004.[13]Lerios, A. et al., Feature-Based VolumeMetamorphosis, in SIGGRAPH 95 - Proceedings,pp 449-456, ACM Press, N.Y, 1995.[14]Berg, A C. Facial Aging in a VirtualEnvironment. Memória de Investigación, UIB, Spain, 2003b.[15]Hall, V., Morphing in 2-D and 3-D, in Dr.Dobb's Journal, July 1993.。

装备环境工程第20卷第8期·90·EQUIPMENT ENVIRONMENTAL ENGINEERING2023年8月重大工程装备气流清扫的超音速喷管气动设计及其性能的对比分析赵宏星1，卢耀辉1，王北昆1，唐波1，罗银生2，陈德君2，毛荣生2（1.西南交通大学 机械工程学院，成都 610031；2.唐山百川智能机器股份有限公司，河北 唐山 063000）摘要：目的提出使用拉瓦尔喷管产生高速气流清扫固体表面附着的水膜。

方法设计中心轴对称锥形喷管（Taper-A）、中心轴对称Sivell法喷管（Sivell-A）、中心轴对称短化喷管（MLN-A）和二维锥形型线喷管（Taper-2D），建立包括外流场的LES数值仿真模型，并进行仿真，分析研究外流场结构，并基于韦伯数判据，分析超音速喷管的清扫性能。

结果喷管短化设计方法可以将喷管长度缩短50%。

外流场速度呈波动衰减趋势，特征线法喷管的气流膨胀更充分。

缩短喷管长度会减小内流场的附面层厚度，因此MLN-A在速度波动中的能量耗散较少；喷管过长也会降低清扫性能，MLN-A和Taper-2D在x L>2区域的最大等效水膜厚度小于0.2 μm，但MLN-A有效清扫面积比Taper-A的喷管提高15%以上，清扫性能最优。

结论喷管短化设计方法可以有效缩短喷管。

喷管结构对清扫性能影响较大。

MLN-A喷管的清扫性能最优。

关键词：超音速喷管；气动设计；大涡模拟；外流场；清扫性能中图分类号：U270.1+1 文献标识码：A 文章编号：1672-9242(2023)08-0090-08DOI：10.7643/ issn.1672-9242.2023.08.012Aerodynamic Design and Performance Comparison of SupersonicNozzle Using Airflow SweepingZHAO Hong-xing1, LU Yao-hui1, WANG Bei-kun1, TANG Bo1, LUO Yin-sheng2, CHEN De-jun2, MAO Rong-sheng2(1. School of Mechanical Engineering, Southwest Jiaotong University, Chengdu 610031, China;2. Tangshan Baichuan Intelligent Machine Co., Ltd., Hebei Tangshan 063000, Chin)ABSTRACT: It is proposed to use a Laval nozzle to generate high-speed airflow to clean the water film attached to the solid surface. A central axisymmetric conical nozzle (Taper-A), a central axisymmetric Sivell nozzle (Sivell-A), a central axisymmet-ric minimum length nozzle (MLN-A), and a two-dimensional conical nozzle (Taper-2D) were designed. An LES numerical收稿日期：2023-03-27；修订日期：2023-05-09Received：2023-03-27；Revised：2023-05-09基金项目：四川省科技计划项目（2022YFG0251）Fund：Sichuan Science and Technology Programme (2022YFG0251)作者简介：赵宏星（1998—），男，硕士研究生，主要研究方向为空气动力学。

专利名称：Volumetric hair simulation发明人：Lena Petrovic,John R. Anderson申请号：US11072049申请日：20050303公开号：US07468730B2公开日：20081223专利内容由知识产权出版社提供专利附图：摘要：A volumetric representation of a hair simulation model determines collective hair attributes. To determine inter-hair collisions, vertices include average velocities of the adjacent portions of the model. The average velocities determine target velocities.Forces for the model are determined from the target velocity values. To direct hair to adesired pose, vertices include target and current density values representing the density of adjacent portions of the model in the desired pose and current position, respectively. The differences in density values determine pressure forces applied to the model. To determine the illumination of the hair, vertices include density values representing the density of adjacent portions of the model. The density values define a hair surface, and signed distance values relative to the surface are determined for the vertices. Normal vectors are determined from the gradients of the signed distance values at locations corresponding the positions of the hairs.申请人：Lena Petrovic,John R. Anderson地址：Oakland CA US,San Anselmo CA US国籍：US,US代理机构：Law Office of Jonathan Hollander PC更多信息请下载全文后查看。

132机械设计与制造Machinery Design&Manufacture第9期2018年9月滚动摩擦系数的测定及EDEM仿真分析刘万锋1，徐武彬$’2,李冰$’2,李玉凤1(1.广丽雅大学机械工程学院，广西柳州'4'〇〇〇;2广酿校临海机械装备设计9造及控制重点实验室培育基地，广西钦州'3'000)摘要：为方便且准确的测定物料滚动摩擦系数，通过对平抛运动装置的改进设定了一种测定物料滚动摩擦系数的方法.将所需测定物料打磨成球状，然后让其在半径为！的圆e铁轨上成！角处无初速度自由滚落，测定其自由停止时水平运 动的距离及水平抛出点的高度，通过计算推导求得滚动摩擦系数。

根据此参数在ED EM中建立安息角仿真模型，通过安 息角试验测定与仿真结果对比验证仿真模型的准确性。

结果表明：此方法所测定的滚动摩擦系数可以用来作为EDEM仿 真参数的标定，此方法所测定的滚动摩擦系可行的。

这为以后物料参数的标定提供了一定的理论基础。

关键词:滚动摩擦系数;安息角;ED EM仿真;S P C法中图分类号:TH16;TH243.1文献标识码:A文章编号：1001-3997(2018)09-0132-04Measurement and Simulation of Rolling Friction CoefficientLIU Wan-feng1，XU W u-bin1，2, LI Bing1，2, LI Yu-Feng1(1.School of Mechanical Engineering，Guangxi University of Science and Technology，Guangxi Liuzhou545000, China;2.Guangxi Colleges and Universities Key Laboratory Breeding Base of Costal Mechanical Equipment Design，Manufacturing and Control，Guangxi Qinzhou535000，China)Abstract：/# order to d etermine the rolling friction coefficient of the materials conveniently and accurately，a method for measuring the rolling friction coefficient of materials is set up with an improved flat polishing motion device.Polishing the materials which needed to be measured into a b a ll，then putting it on at an angle o f! in an arc track with the radius as the arc of R and let it free roll without initial velocity，measuring the distance of the horizontal movement and the height of the initial position.Set up a Repose angle simulation model in EDEM，comparing with the result of repose angle from experiment and simulation to verify the accuracy of the simulation model.The results shows that the error ratio between experimental dates and simulation results less than 3]，it is also a valuable method to test the rolling f riction coefficient using the rolling f riction coefficient from the experiment as the calibration of simulation parameters in EDEM.It provide a theoretical basis for the calibration of material parameters in the future.K e y W o r d s：Rolling Friction Coefficient；Repose A ngle；E D E M Simulation；S P C M e t h o d1引言颗粒介质是人类生产生活中常见的一种物质组态，其离散现出一些不同与一般体体体质<1-3=，但是人们对颗粒物料的认识还远远没有到达成熟的地步。

第二章室内外计算参数第一节一般术语第2.1.1条计算参数 design conditions特指设计计算过程中所采用的表征空气状态或变化过程及太阳辐射的物理量。

常用的计算参数有干球温度、湿球温度、含湿量、比焓、风速和压力等。

第2.1.2条室内外计算参数 indoor and outdoor design conditions设计计算过程中所采用的室内空气计算参数、室外空气计算参数和太阳辐射照度等参数的统称。

第2.1.3条空气温度 air temperature暴露于空气中但又不受直接辐射的温度表所指示的温度。

一般指干球温度。

第2.1.4条干球温度 dry-bulb temperature干球温度表所指示的温度。

第2.1.5条湿球温度 wet-bulb temperature湿球温度表所指示的温度。

第2.1.6条黑球温度 black globe temperature黑球温度表所指示的温度。

第2.1.7条露点温度 dew-point temperature在大气压力一定、某含湿量下的未饱和空气因冷却达到饱和状态时的温度。

第2.1.8条空气湿度 air humidity表征空气中水蒸汽含量多少或干湿程度的物理量。

第2.1.9条绝对湿度 absolute humidity单位体积的湿空气中所含水蒸汽的质量。

第2.1.10条相对湿度 relative humidity空气实际的水蒸汽分压力与同温度下饱和状态空气的水蒸汽分压力之比，用百分率表示。

第2.1.11条历年值 annual(value)逐年值。

特指整编气象资料时，所给出的以往一段连续年份中每一年的某一时段的平均值或极值。

第2.1.12条累年值 normals多年值。

特指整编气象资料时，所给出的以往一段连续年份的某一时段的累计平均值或极值。

第2.1.13条历年最冷月 annual coldest month每年逐月平均气温最低的月份。

第2.1.14条历年最热月 annual hottest month每年逐月平均气温最高的月份。

第42卷第2期2008年2月浙 江 大 学 学 报(工学版)Journal o f Zhejiang U niv ersity (Engineer ing Science)Vol.42No.2Feb.2008收稿日期:2006 09 17.浙江大学学报(工学版)网址:w w w.journals.z /eng基金项目:国家自然科学基金资助项目(10572125).作者简介:马利(1975-),女,浙江临海人,博士后,从事计算力学研究.E mail:mali@通讯联系人:陶伟明,男,教授,博导.E mail:taowm @z SPH 耦合有限元方法的水射流弹塑性碰撞模拟马 利1,陶伟明2,郭乙木2,郑津洋1(1.浙江大学化工机械研究所,浙江杭州310027; 2.浙江大学工程力学系,浙江杭州310027)摘 要:针对有限元分析(FEA )处理流固耦合及超大变形问题时所存在的困难,提出了光滑粒子流体动力学(SPH )耦合有限元的方法来模拟水射流与刚性表面和塑性表面的碰撞过程.在SP H 计算中,采用计算几何中的Vo ro no i 图给每个SP H 粒子赋予质量,并采用可变光滑长度考虑大变形中粒子相互间距的变化.分析了水射流与刚性表面的弹性碰撞,给出了射流变形与压力传播过程.对于水射流与塑性表面的碰撞,将SP H 方法与有限元方法通过接触算法加以耦合,其中水射流以SPH 粒子建模,而目标体以有限元方法建模.通过该耦合模型研究了在不同射流速度下切割深度与成坑孔径的变化,为理解切割机理和优化工作参数提供了参考.关键词:水射流切割;光滑粒子流体动力学(SP H);有限元分析(FEA );碰撞中图分类号:O 39 文献标识码:A 文章编号:1008 973X(2008)02 0259 05Elastic /plastic impact simulation of water jet using smoothed particlehydrodynamics and finite element methodM A Li 1,T AO Wei ming 2,GU O Yi mu 2,ZH ENG Jin yang 1(1.I ns titute of M echanical M achiner y and P r ocess Eq uip ment,Zhej iang Univ er sity ,H angz hou 310027,China;2.D ep ar tment o f M echanics ,Zhej iang Univ er sity ,H angz hou 310027,China)Abstract:T o overco me the difficulty of fluid so lid interaction and extra larg e defo rmatio n pro blem o f finite element analy sis (FEA),this wo rk presented a coupled smoothed particle hydrody namics (SPH )and finite element method to sim ulate the process of w ater jet impacting w ith the rig id and plastic tar get material.In the SPH com putation,the Vo ronoi diagr am in com putational geometry w as adopted to assig n the initial particle m ass and the variable sm oothed leng th w as used to co nsider the large distance chang e amo ng the particles induced by the large defor mation.According to the elastic im pact com putational result,the pro cess of w ater jet colliding with the rigid target and then becoming into the state of flow w as analyzed,also the defor mation pro cess of the w ater jet w ith the pr essure propagatio n w as given.And in the plastic impact co mputation,the SPH w as combined w ith FEA by contact alg orithm,w here the w ater jet w as molded by SPH particles and the target mater ial w as molded by finite elem ents.The relations of cut depth and crater size w ith different w ater jet velocity w ere studied,w hich w as helpful to understand the w ater jet cutting mechanism and optimize the o perating parameter s.Key words:w ater jet cutting;sm oothed particle hydrodynam ics (SPH );finite element analy sis (FEA);impact pr ocess高速液体与固体表面的碰撞广泛存在于自然界和工业生产活动中,例如雨滴打在高速飞行的飞机蒙皮上,有可能导致材料的破坏;在材料加工技术中,水射流切割技术就是利用高压水从精细小孔射出,从而产生强大的切割力来切割金属或复合材料,这些都是典型的液体高速冲击固体的现象.水射流切割模拟的困难在于找到一种有效的方法来处理这类流固耦合的问题,不仅包括界面接触问题,而且涉及流体的大运动/变形、固体的断裂等.虽然在流体与固体领域已经各自发展了成熟的数值计算方法,如有限差分法、有限元法等,但是应用有限元法处理上述液体高速冲击固体问题,由于涉及极度大变形,如液体的飞溅和固体的断裂等,在计算过程中可能会引起网格的严重扭曲或失效,从而影响计算精度.即使基于任意拉格朗日 欧拉描述(arbitrary Lagrange Euler ian form ulation),对一个复杂的三维问题不断进行网格的重构以及接触界面的识别和处理,在计算上仍然是很困难的.光滑粒子流体动力学(smoothed particle hy drodynam ics,SPH )方法是近20多年来逐步发展起来的一种无网格方法,基本思想是利用相互作用的质点组来描述连续流体(或固体),各个物质点上承载各种物理量,包括质量、速度等,通过求解质点组的动力学方程和跟踪每个质点的运动轨道,求得整个系统的力学行为.最初的SPH 方法是Lucy [1]、Gingo ld 等人[2]于1977年分别提出的.Jo hnso n 等人[3]提出了归一化的光滑函数算法,该算法提高了SPH 的计算精度,并能够通过分片试验.他们将SPH 法应用于圆杆的撞击以及弹体侵彻的数值计算,得到了很好的结果[4 5].在流体方面,Mo nag han[6]应用SPH 方法计算了水坝突然坍塌后水体的运动.本文通过综合SPH 方法和有限元方法来模拟水射流与刚性表面和塑性表面的碰撞过程,分析了水射流与刚性表面的碰撞,给出了射流变形与压力传播过程.在塑性碰撞计算中,分别用SPH 粒子和有限元来模拟水射流与目标体金属,通过 no des to surfaces 的接触算法将SPH 与FEA 进行耦合,计算了在不同射流速度下的切割效果,包括切割深度与成坑孔径,这将有助于理解水射流切割机理,进行生产上的参数优化.1 SPH 基本理论SPH 方法的核心思想是利用一个作用在近程范围内的核函数描述空间的粒子及其邻粒子的相互关系.设F(x)为一场函数,则有F h(x I )=w (x I-x J,h)F (x J)d J.(1)式中:w (x I -x J ,h)为核函数,h 为核函数的影响半径.关于核函数应满足的条件可参见文献[2].将式(1)化为离散化求和形式:F (x I )=NJ =1w IJm J J F (x J ).(2)式中: 、m 分别为射流的密度和质量.得到F (x I )=NJ =1m J Jw IJ F (x J ).(3)对密度场应用式(2),有(x I )=NJ =1wJm J .(4)在流体动力学中,可用如下方程组来描述流体的运动和状态:运动方程 = d vd t,连续方程 dd t + v =0.能量方程 d E d t = d d t.(5)式中: 、v 、E 、 分别为射流的应力、速度、内能和应变.对于式(5)中的运动方程,应用SPH 方法并加以离散化得d v I d t =1 I NJ =1m JJ J w IJ .(6)式(6)通过一定形式的变形得到其对称形式:d v I d t=NJ =1m JI 2I + J2Jw IJ .(7)根据式(7),可由t n时刻的应力求得t n时刻的加速度,进而由如下显式中心差分格式:vn +1/2I =vn -1/2I+12( t n +1/2+ t n -1/2)d v nI d t,x n +1I=x n I+ t n +1/2v n +1/2I,(8)求出tn +1时刻的位移与速度,从而得到t n +1/2时刻的应变率为d n +1/2Id t=12[ v +( v )T ]n +1/2I .(9)对于式(5)中的连续方程,类似以上推导,有 n +1I= n I + tn +1/2NJ =1m J(v I-v J ) w IJn +1/2.(10)再结合流体的本构方程和能量方程,可得t n +1时刻的应力 n +1I 和内能E n +1I .SPH 常用的核函数有B 样条函数和二次光滑函数,根据文献[4]的讨论,选用性能较优的Quad ratic 光滑函数: w IJ =1 h 3IJ38r 2IJ -32r IJ +32;0 r 2.0.(11)式中:h 为核函数作用范围的尺寸度量,也称为光滑260浙 江 大 学 学 报(工学版) 第42卷长度;r 为粒子间距与光滑长度的比值.w IJ 的一阶导数为d w IJ d r IJ =1h 3IJ34r IJ -32;0 r 2.0.(12)2 可变光滑长度对于复杂流体大变形过程,粒子的相互距离在每一时间步都发生很大的变化,因此计算过程中必须考虑随粒子间距的变化而重新计算合适的光滑长度.光滑长度与粒子的分布有密切的关系,因此光滑长度对时间的导数应遵循连续性方程变化规律:d (h(t))d t=h(t) v .随着时间的推移,每个粒子的光滑长度开始变得不同.当粒子相互远离时,其光滑长度增加;当粒子相互靠近时,其光滑长度也相应减少.一般保持每个粒子的邻域内所包含的邻粒子个数大致相当.3 初始质量分配在无网格计算中,Organ [7 8]等人利用计算几何中的Voro noi 多边形作为积分网格.对于SPH 方法,较难确定的是每个粒子所占的空间体积,此处借用Voronoi 图定义初始时刻每个粒子的质量.当然也可以采用规则排列的粒子,特别对于立方体型的模型.但是对于较为复杂的区域或形状边界,Voronoi 图具有更好的适应性.图1 Voronoi 图示意F ig.1 Schematic of V o ronoi cell过相邻粒子的连线作垂直平分线,这些垂直平分线相交而构成的围绕粒子的多边形即为Voro noi 图,如图1所示.根据Vor ono i 多边形的面积,确定每个粒子所负载的质量.计算中,每个粒子承载的质量为m I = S I H .式中: S I 为Vo ronoi 多边形的面积, H 为每层粒子沿轴向所占的高度.对于最外层的粒子,无法按照以上规则构筑其Vor onoi 多边形,因此,当实际建模时可在外围多布置一层粒子,然后截取所需要的粒子及体元结构.当初始建模时,应该事先估算粒子的个数与分布,尽量使每个粒子的排列间距均匀.在初始状况下粒子排列严重不均匀会对计算的稳定性和可靠性产生不利影响.4 数值模拟4.1 水射流与刚性表面的碰撞建立了水射流冲击刚性平面的SPH 数值模拟模型.假设射流的形状是圆柱形,半径为0.24cm,高为0.75cm,以初始速度100m/s 垂直撞击刚性平面.粒子模型如图4(a )所示,沿长度方向分为16层,每层圆截面上布置92个粒子,共计1472个粒子.计算了水射流与平面接触开始后90 s 内的变形过程,图2给出了不同时刻水射流变形状况,以及伴随这一过程的压力变化.从不同时刻的变形图像可以看到整个射流撞击刚性平面并在平面上发生流动的过程.在碰撞开始时,底面压力最先增加到p 4.57 107Pa ,在随后的时间内压力逐渐向圆柱顶部传播,并在传播过程中逐渐降低.对底部液体来说,在碰撞瞬间整个圆柱底部压力升到p ,但此时只有边缘部分的液体可以自由横向流动,靠近轴线部分仍保持压缩状态,直到圆柱周边的稀疏波到达为止.稀疏波到达以后,刚性平面上的压力降到稳定流动状态.图3给出了底面上沿径向分布的粒子的速度v R 的变化,由内到外其编号分别为1、2、8、21、39、63,粒子1是位于底面轴线上的节点,粒子63是位于底面周线上的节点.可以看到,边缘上的液体在碰撞瞬间即开始向外流动,越靠近轴线的液体获得流动需要的时间越滞后于周线液体,并且获得的流动速度逐渐降低.从第40 s 左右开始,相邻粒子的速度线由互相平行到互相纽结,说明低速流体被相邻的高速流体加速,反之高速流体被减速,这种现象可以认为是流体进入了类似湍流状态.但是这种类似于湍流的行为究竟是反映了流动的本质还是由于SPH 方法数值计算存在的误差而造成的,尚需进一步研究.4.2 水射流与塑性表面的碰撞(水射流切割)水射流切割实际上就是射流与塑性表面的碰撞过程.对这一过程进行数值模拟,将SPH 方法和有261第2期马利,等:SPH 耦合有限元方法的水射流弹塑性碰撞模拟图2 不同时刻射流变形和压力传播过程F ig.2 Defor mation process w ith pressur e pr opagatio nat differenttime图3 不同节点粒子径向速度变化曲线Fig.3 Cur ves of r adial velocit y of par ticles at differ ent nodes限元方法通过 no des to sur faces 接触算法进行耦合[9 12].通过几何条件判断SPH 粒子与有限元单元表面是否接触,若满足接触条件则由罚函数方法确定接触力,这样就将SPH 粒子作用力施加到了有限元单元表面.根据水射流切割的实际工作参数,建立了高为3mm ,半径为0.62mm 的水射流SPH 粒子模型,目标体金属以1010钢为例,尺寸为10mm 10m m 5m m ,相关的材料属性如下:材料密度 t =7845kg/m 3、弹性模量E =207GPa 、泊松比 =0.29、屈服应力 s =128.5M Pa ,失效应变 f =0.2.利用耦合SPH 与FEA 方法分别对不同射流速度下的水射流切割过程进行了数值模拟,图4为当射流速度v =1100m /s 时的数值模拟结果.在数值模拟过程中,目标材料的失效由塑性应变阈值控制,当单元的塑性应变达到失效应变值 f 时,可以认为单元失效而在计算中将该单元删除.根据数值模拟的结果,可以这样描述模拟切割机理:首先SPH 粒子与目标体金属表面发生碰撞,在碰撞中心区域形成初始孔径,随着碰撞过程的延续,一部分粒子继续加深孔穴,同时一部分粒子与孔穴的壁面碰撞而扩大成坑.作为有限长度射流,射流的头部在碰撞过程中形成 蘑菇状 ,当应力波到达射流尾部时,由于界面的反射,压缩波转变为拉伸波,导致射流尾部液体向上飞溅.图4 水射流切割模拟F ig.4 Simulatio n o f water jet penet ratio n图5给出了射流速度与最大切割深度的关系.可见,最大切割深度随射流速度的增加而增大,两者的关系近似为线性关系.但是当射流速度超过1000m/s 时,切割深度与成坑孔径均显著增加.图6是射流速度与成坑孔径的关系.可见,随着射流速度的增加,成坑孔径也基本上线性增加,这与实验现象一致.大的成坑孔径意味着切缝质量不高,因此,应该在保证切割深度的前提下,尽量采用低的射流速度.262浙 江 大 学 学 报(工学版) 第42卷图5 切割深度与射流速度的关系Fig.5 Relationship betw een w ater jet velocity and cutdepth图6 成坑孔径与射流速度的关系F ig.6 R elatio nship betw een w ater jet velo city andcr ater diameter5 结 语本文应用SPH 方法与FEA 方法研究了在水射流切割过程中射流与刚性表面和塑性表面的碰撞过程,分别分析了与刚性表面碰撞过程中水射流的变形与压力变化,并在塑性碰撞过程中计算了射流速度与切割深度及成坑孔径的关系,为了解水射流切割机理和优化工作参数提供了参考数据.数值模拟的结果显示了这一方法在处理流体(或固体)大变形问题时的有效性.利用SPH 粒子建立水射流模型,并与有限元模型进行耦合,建立的耦合模型适用于处理类似的流固耦合问题.整个计算都是基于拉格朗日系下的描述,较之欧拉描述或ALE 描述更为方便有效.参考文献(References):[1]L U CY L B.A numerical approach to the testing of fusion process [J].The Astronomical Journal,1977,88:1013 1024.[2]GI NG OL D R A ,M ON A GH A N J J.Smoo thed particlehy dr odynamics:t heo ry and application t o no n spherical star s [J].Monthly Notice Royal Astronomy Society,1977,181(4):375 389.[3]JO HN SO N G R,BSISSEL S R.N or ma lized smoothingfunctions fo r SPH impact co mputatio ns [J].Interna tional Journal of Numerical and Mathematical Engineer ing,1996,39(16):2725 2741.[4]JO HN SO N G R,ST RYK R A ,BSISSE S R.SP H forhigh velocity impact co mputations [J].Computer Meth ods in Applied Mechanics and Engineering,1996,139(1 4):347 373.[5]JO HN SO N G R,BSISSEL S R ,ST RY K R A.A generalized particle alg or ithm for high v elocit y impact co mpu tations [J].Computational Mechanics,2000,25(2/3):245 256.[6]M ON A GH A N J J.Simulating free sur face flow s w ithSPH [J].Journal of C omputational Physics,1994,110(2):399 406.[7]OR GA N D,FL EM IN G M ,T ERR Y T .Continuousmeshless appr ox imations for nonco nv ex bo dies by diffractio n and t ranspar ency [J].C om putational Mechan ics,1996,18(3):225 235.[8]ZH OU J X,W EN J B,ZH A N G H Y,et al.A no dalinteg ratio n and post processing technique based on V or onoi diagr am fo r Ga lerkin meshless methods [J].Computer Methods in Applied Mechanics and Engineer ing,2003,192(35/36):3831 3843.[9]CA M P BEL L J,V IG N JEV IC R,L IBERSK Y L.A contact algo rithm fo r smo othed part icle hydrodynamics [J].C omputer Methods in Applied Mechanics and Engineer ing,2000,184(1):49 65.[10]JOH N SO N G R.L inking of L ag rangian par ticle methods to standard finite element methods fo r high velocity impact computatio ns [J].Nuclear Engineering and D e sign,1994,150:265 274.[11]AT T A WA Y S W ,HEIN ST EIN M W ,SW EGL E JW.Co upling of smo oth par ticle hydro dy namics w ith the f inite element method [J].Nuclear Engineering and Design,1994,150:199 205.[12]BO NET J,KU L A SEGA RA M S,ROD RIGU EZ P A ZM X,et al.Var iational for mulatio n fo r the smoo th par ticle hydrodynamics (SP H)simulatio n of fluid andso lid pro blems [J].C omputer Methods in Applied M e chanics and Engineering,2004,193(12 14):1245 1256.263第2期马利,等:SPH 耦合有限元方法的水射流弹塑性碰撞模拟。

requirement of comfort for outdoor activities. In order to improve the outdoor micro-climate, a canopy that covered the project site was proposed. The simulation results showed that optimizing the shading coefficient of the canopy while designing the openings forms can improve the overly heating conditions during summer.Key words outdoor microclimate; thermal sensational index; Ladybug & Honeybee; computational fluid dynamics (C FD)近年来，城市空间布局中室外或半室外的活动场所逐渐增加，鼓励市民增加室外活动、亲近自然，为健康生活创造了良好的物理空间条件。

因此，热舒适度成为室外或半室外活动场所的评价要点之一。

极端的热环境可能对室外活动人员造成感官上的不适，更有甚者会对身体造成一定的伤害。

例如：降低器官生理机能、提升心血管疾病的发病率等[1，2]。

在全球气候变暖的大趋势下，更加需要关注城市微气候中人员活动的安全性和舒适性[3]。

影响室外热环境的因素较为复杂。

其中，一类因素是自然环境条件，主要包括风速、太阳辐射等；另一类因素是人为设计，主要包括建筑群的布局安排、下垫面的形式等[4]。

除了前述客观因素以外，有的评价标准也将人员主观行为列为影响因素之一[5]。

针对人员热舒适，国内外已有多种舒适度评价指标，各个评价指标包含的参数各不相同，且评价标准的侧重点各不同。

本文将对常见的几种室外热舒适评价指标进行文献综述对比，并通过结合实际项目案例进行室外热舒适度的模拟分析计算。

Tribology International 33(2000)47–58/locate/tribointSimulation of rough surfaces with FFTJiunn-Jong Wu*Department of Mechanical Engineering,Chinese Cultural University,Taipei,Taiwan,ROC Received 30June 1999;received in revised form 1February 2000;accepted 7February 2000AbstractA numerical procedure for the simulation of three-dimensional surfaces has been developed.The method is based on Fast Fourier Transform (FFT).This method can simulate surfaces with given spectral density or auto-correlation function (ACF).Although this method cannot guarantee that each profile of the generated surface has the correct ACF,the average ACF of profiles for the generated surfaces is very close to the given one.©2000Elsevier Science Ltd.All rights reserved.Keywords:Simulation;FFT;Surface roughness1.IntroductionThe measurement and characterization of rough sur-faces are important in tribology.Numerous methods have been developed.In many applications and analyses,a growing need is developing not only to characterize but also to generate or synthesize surfaces with desired properties,especially spectral density or auto-correlation function (ACF),in more than two directions to meet the engineering requirement.Many synthesizing methods have been developed to meet the requirements,such as ACF or producing a sur-face with a given spectral density.But also,they have their restrictions.Shinozuka and Jan [1],in 1972,proposed a multivari-ate–multidimensional random process with a known spectral density.In 1978,Patir [2]used ‘the linear trans-formation on random matrixes’and proposed a pro-cedure for generating an N ×M matrix of roughness amplitudes z i,j ,having a Gaussian distribution of heights,with a given n ×m auto-correlation matrix.Besides,time series and FFT (Fast Fourier Transform)are popular in generating surfaces.In 1982,Watson and Spedding [3]used ARMA (Auto-Regressive and Moving Average method).Newland [4],in 1984,used FFT.In 1990,Gu and Huang [5]used ARMA,but generated non-Gaussian*Tel.:+886-2-28239548.E-mail address:jjw5277@ (J.-J.Wu).0301-679X/00/$-see front matter ©2000Elsevier Science Ltd.All rights reserved.PII:S 0301-679X (00)00016-5surfaces.In 1991,You and Ehmann [6]used FFT and ARMA.In 1992,Hu and Tonder [7]used finite impulse response (FIR)filters.In 1998,by using Hu and Tond-er’s method,Chilamankuri and Bhushan [8]found that the average ACF of profiles for generated surface is very close to the desired one.In 1995,Ganti and Bhushan [9]used FFT to generate a one-dimensional fractal profile.All these methods can be applied for specific purposes.Among these methods,it is found that current time series methods (ARMA)can consider few orders system,which means only the features of ACF near the origin are simulated.Also it is found that FFT is the fast and convenient tool in generating surfaces.Therefore,in this paper,we focus our attention on the method using FFT to generate surfaces with required properties.Thus,Newland’s and Hu and Tonder’s work will be discussed in Section 2.Finally,a new simulation method is pro-posed to simulate rough surfaces with given spectral density or ACF.2.Generation of a rough surface with IFFT 2.1.Aliasing of FFTSuppose there is an infinite series {z r |r =0,1,2,ϱ}.The FFT can be applied to a finite series {z r |r =0,1,2,%,N Ϫ1}only.FFT is defined as Z˜k ϭ1N ͸N Ϫ1r ϭ0z r exp ͩϪi 2p kr Nͪfor k ϭ0,1,2,%,N Ϫ148J.-J.Wu/Tribology International33(2000)47–58NomenclatureA k,l FFT of random series h r,sb Correlation distanceC r,s Real circular auto-correlation function of afinite surfaceC˜r,s Circular auto-correlation function computed from surface or used to generate surfaces C¯r,s Mean circular auto-correlation function obtained from generated surfacesh r,s FIRfilter in Hu and Tonder’s methodk,l,p,q,r,s Variable integerH k,l Transfer function in Hu and Tonder’s methodM Number of surface points in x directionm Number of ACF values in x direction,used by Hu and Tonder’s methodN Number of surface points in y directionn Number of ACF values in y direction,used by Hu and Tonder’s methodR r,s Real auto-correlation functionR˜r,s Auto-correlation function obtained from generated surfaces,or used to generate surfaces R¯r,s Mean auto-correlation function obtained from generated surfaces⌬R p,q Difference between real ACF and the ACF of surface generated by Hu and Tonder’s method S r,s Real spectral densityS˜r,s Spectral density obtained from generated surfaces,or used to generate surfacesz r Profile heightz r,s Surface heightZ k FFT of infinite series z r(r=0,1,2,%,ϱ).Z˜k FFT offinite series z r(r=0,1,2,%,NϪ1).Z p,q FFT of surface height z r,s(r=0,1,2,%,MϪ1,s=0,1,2,%,NϪ1).h i,j A series of random number⌬Sampling intervalf k Phase of Z k,or a random phasew Wave numbers Standard deviation of surface heightswhich is the Fourier Transform for afinite series {z r|r=0,1,2,%,NϪ1}While doing the Fast Fourier Transform(FFT),there is aliasing happening.That is,the high frequency components contribute to the{z r}series and distort the Fourier coefficients calculated by the FFT below Nyquist frequency(p/⌬,where⌬is the sampling interval).If the Fourier Transform of an infinite series z r is Z k(r,k=0,1,2,ϱ),then Z˜k and Z˜N−k are the sum of contri-butions from the Fourier Transform Z at frequencies k, NϪk.Z˜kϭZ˜N−kϭZ kϩZ N−k(1) Therefore,its spectral density obtained from FFT isS˜kϭS˜N−kϭ|Z˜k|2ϭ|Z˜N−k|2ϭ|Z kϩZ N−k|2ϭ|S1/2k exp(i f k)ϩS1/2N−k exp(i f N−k)|2ϭS kϩS N−kϩ2S1/2k S1/2N−k cos(f kϪf N−k)where f k and f NϪk are the phases of Z k and Z NϪk,and S k is its real spectral density.While doing FFT of afinite series z r(r=0,1,2,%,NϪ1),we can only know Z˜k.It is impossible to know Z k and Z NϪk from Z˜k and Z˜N−k.But we can be sure that(|Z k|Ϫ|Z N−k|)Ͻ|Z˜k|ϭ|Z˜N−k|Ͻ(|Z k|ϩ|Z N−k|)(S1/2kϪS1/2N−k)ϽS˜1/2kϭS˜1/2N−kϽ(S1/2kϩS1/2N−k)Therefore,by using FFT,we cannot tell Z k from Z NϪk and also S k from S NϪk.2.2.ACF and circular ACFSuppose the ACF of the infinite series{z s|s=0,1,2,ϱ} is R r.R rϭE{z s z s+r}sϭ0,1,2,%,ϱWhile doing FFT,the circular auto-correlation func-tion(cir-ACF)C r of thisfinite series can be obtained by[4]C rϭ1N͸NϪ1sϭ0z s z s+r49J.-J.Wu /Tribology International 33(2000)47–58where z s +r =z s +r −N for s +r ՆN .The cir-ACF C r can also be obtained from its true ACF R r .C r ϭC N −r ϭN −r N R r ϩrNR N −rIt is also the IFFT of the spectral density of the profile.C r ϭ͸N Ϫ1r ϭ0S˜k exp ͩi 2p kr N ͪ2.3.FormulaFrom the above properties,we propose that a simu-lated surface can be generated byz p ,q ϭ͸M Ϫ1k ϭ0͸N Ϫ1l ϭ0ͱS˜k ,l exp ͩi 2p ͫf k ,l ϩkp M ϩlq Nͬͪp ϭ0,1,2,%,(M Ϫ1)(2)q ϭ0,1,2,%,(N Ϫ1)where f k ,l is a set of independent random phase anglesuniformly distributed between 0and 2p .And in order to make z to be real,there is a restriction that f k ,l ϭϪf M −k ,N −l k ϭ0,1,2,%,M /2f k ,0ϭϪf M −k ,0l ϭ0,1,2,%,N /2f 0,l ϭϪf 0,N −lf 0,0ϭf 0,N /2ϭf M /2,0ϭf M /2,N /2ϭ0S˜k ,l is the spectral density used to generate the surfaces.For a two-dimensional surface,the spectral density S ˜k ,l could be obtained from its real spectral density by S˜M −k ,N −l ϭS ˜k ,N −l ϭS ˜M −k ,l ϭS ˜k ,l ϭS k ,l k ϭ0,1,2,%,M /2l ϭ0,1,2,%,N /2Then,by using Eq.(2),a surface can be generated.If we want to obtain a surface with given ACF,then R˜N −r ,M −s ϭR ˜r ,M −s ϭR ˜N −r ,s ϭR ˜r ,s ϭR r ,s r ϭ0,1,%,M /2s ϭ0,1,%,N /2S˜k ,l can obtained from R ˜r ,s S ˜k ,l ϭ1MN ͸M Ϫ1r ϭ0͸N Ϫ1s ϭ0R˜r ,s exp ͩϪi 2p ͫkr M ϩls N ͬͪk ϭ0,1,%,(M Ϫ1)(3)l ϭ0,1,%,(N Ϫ1)Then,by using Eq.(2),a rough surface can be simulated with given ACF.In fact,Ganti and Bhushan’s method [9]is very simi-lar to ours.But their method can only be applied to aone-dimensional profile,but can not be extended into a two-dimensional surface.parison with current method2.4.1.Newland’s methodNewland also generated rough surface with Eq.(2).But they have different ideas about spectral density and ACF.In Newland’s method,the circular auto-correlation function (cir-ACF)is used.If ACF R r ,s of a surface isknown,the new ACF R˜r ,s could be obtained by R ˜r ,s ϭC r ,s ϭ(M −r )M (N −s )N R r ,s ϩr M (N −s )N RM −r ,sϩ(M −r )M s N R r ,N −s ϩr M sN R M −r ,N −s (4)r ϭ0,1,%,M /2s ϭ0,1,%,N /2The value of R˜r ,s is equal to cir-ACF C r ,s .Then by using FFT,the spectral density S˜k ,l can be obtained from R˜r ,s by Eq.(3).Then surface can be generated by Eq.(2).If the spectral density is known at the very beginning,Newland believes that there is a circular relation similar to Eq.(4)for spectral density.Hence,the spectral den-sity can be obtained byS˜k ,l ϭ(M −k )M (N −l )N S k ,l ϩk M (N −l )N S M −k ,l ϩ(M −k )M l N S k ,N −l ϩk M l N SM −k ,N −lk ϭ0,1,%,(M Ϫ1)l ϭ0,1,%,(N Ϫ1)where S k ,l is the true spectral density.Then,as before,the profile can be obtained by Eq.(2).Hence,z p ,q shouldhave the spectral density S˜k ,l .2.4.2.Hu and Tonder’s methodHu and Tonder [7]proposed an FIR (finite impulse response filter)method.Supposed {h i ,j }is a series of Gaussian random numbers,and the simulated rough sur-face heights are taken asz p ,q ϭ͸m Ϫ1r ϭ0͸n Ϫ1s ϭ0h r ,s h p +r ,q +sp ϭ0,1,2,%,M Ϫ1(5)q ϭ0,1,2,%,N Ϫ1where h r ,s is called FIR filter and is the coefficient which50J.-J.Wu /Tribology International 33(2000)47–58defines the system.Thus,an M ×N -point surface with a given m ×n ACF is generated.m and n should be less than M and N .Usually,they are chosen to be m =M /2and n =N /2.Assume that h r ,s and H k ,l have the following relation-ship:(6)H k ,l ϭ͸m Ϫ1r ϭ0͸n Ϫ1s ϭ0h r ,s exp ͩi 2pͫkr m ϩls nͬͪh r ,s ϭ1mm ͸m Ϫ1k ϭ0͸n Ϫ1l ϭ0H k ,l exp ͩϪi 2pͫkr m ϩls nͬͪk ,r ϭ0,1,2,%,(m Ϫ1)l ,s ϭ0,1,2,%,(n Ϫ1)Therefore,the FFT of Eq.(5)gives Z k ,l ϭH k ,l A k ,lwhere A k ,l is the FFT of h ,and H k ,l is called the trans-fer function.Because h is a random series,its spectral density is a constant (S hh ).S hh ϭAA Јϭconstantwhere A Јis the conjugate of A .Similar to Eq.(3),by using FFT,S˜k ,l can be obtained from a new ACF R˜r ,s ,which can be obtained by R˜m −r ,n −s ϭR ˜r ,n −s ϭR ˜m −r ,s ϭR ˜r ,s ϭR r ,s r ϭ0,1,%,m /2s ϭ0,1,%,n /2Thus,the transfer function could be obtained from H k ,l ϭͫS˜k ,l S hhͬ1/2Then,h r ,s could be obtained from Eq.(6),so the sur-faces can be generated by Eq.(5).Hu and Tonder also proposed another method to implement the above pro-cedure by FFT.However,there are mathematical mis-takes in their method (see Appendix A.1).3.ExampleWe generate a 50×50-point surface with given spectraldensity S (w x ,w y )=1/(w 2x +w 2y )and sampling interval 0.01.The perspective view of the generated surface is shown in Fig.1.Fig.2shows its height distribution.It is found that the height follows Gaussian distribution.As for the spectral density and ACF of simulated sur-face,we use one-dimensional profiles ortwo-dimen-Fig.1.Simulated surface.sional surfaces as examples.We compare the spectral density and ACF of profiles or surfaces generated by Newland’s,Hu and Tonder’s and our methods.3.1.From spectral densityNewland and our method can start from given spec-tral density.Fig.3shows the spectral densities of profiles that gen-erated by our method and Newland’s method with given spectral density S (w )=1/(w 2+1).Both methods can gen-erate a surface with good spectral density.But from the figure,we find that our method is better than Newland’s.From Fig.4,we find that,in Newland’s method,there is about 30%error at some frequency.But in our new method,there is about 10%error only.Fig.5shows similar spectral densities with given spectral density S (w )=1/(w 3+1).It is also found that our method is better thanNewland’s.Fig.2.Height distribution of Fig.1.51J.-J.Wu /Tribology International 33(2000)47–58Fig.3.Spectral density of profiles generated with S (w )=1/(w 2+1).Fig.4.S (w )generated profile /S (w )true of Fig.3.Fig.5.Spectral density of profiles generated with S (w )=1/(w 3+1).3.2.Profile generated from ACFIn this section,profiles are generated from the given ACF R (X )=exp(ϪX /b ),where b is the correlation dis-tance,which makes ACF to fall to 1/e .We generate pro-files by Hu and Tonder’s,Newland’s and our methods with this given ACF.For each method,we choose the correlation distance b to be 50,10and 2.Therefore,a total of nine cases are investigated.For each case,five profiles are generated.Thus,the variation of the ACF for each case can be observed.We will compare the ACF of generated profiles with the given ACF.The ACF of generated profile is defined by R r ϭ1N −r͸N Ϫ1Ϫrs ϭ0z s z s +rWhile generating profiles,we start from a 512-point ACF.Therefore,for Hu and Tonder’s method,n =512,and N =1024(see Appendix A.1).That is,1024-point profiles are generated.For Newland’s and our methods,N =512(see Appendix A.2).That is,512-point profiles are generated.For all three methods,the sampling inter-val is 1.Figs.6–8shows the ACF of profiles generated by Hu and Tonder’s method with correlation distance 50,10and 2,respectively.In Fig.6,it is found that ACF does not meet the requirement.Fig.7is slightly better.Fig.8is much better,but it is still not good enough.Strictly speaking,it is not good at all.Therefore,we can assert that profiles generated by Hu and Tonder’s method do not have the correct ACF most of the times.Figs.9–11show the ACF of profiles generated by Newland’s method with correlation distance 50,10and 2,respectively.We find that ACF is good for X Ͻ20.Figs.12–14show the ACF of profiles generatedbyFig.6.ACF’s of five 1024-point profiles generated by Hu and Tond-er’s method with R (X )=exp(ϪX /50),sampling interval 1.52J.-J.Wu /Tribology International 33(2000)47–58Fig.7.ACF’s of five 1024-point profiles generated by Hu and Tond-er’s method with R (X )=exp(ϪX /10),sampling interval1.Fig.8.ACF’s of five 1024-point profiles generated by Hu and Tond-er’s method with R (X )=exp(ϪX /2),sampling interval 1.our methods with correlation distance 50,10and 2,respectively.While generating profiles,we use the same random series as those of Newland’s method.The ACF of profiles generated by our method are not much differ-ent from those by Newland’s method.ACF is good for X Ͻ20,too.Generally speaking,Newland’s and our methods are better than Hu and Tonder’s method,but are not good enough.For a 512-point profile,ACF is good for X Ͻ20only.It is acceptable,but not very good.3.3.Surface generated from ACFChilamankuri and Bhushan [8]use Hu and Tonder’s method to generate 256×256-point surfaces.They found that the average ACF are good.Therefore,we usedaFig.9.ACF’s of five 512-point profiles generated by Newland’s method with R (X )=exp(ϪX /50),sampling interval1.Fig.10.ACF’s of five 512-point profiles generated by Newland’s method with R (X )=exp(ϪX /10),sampling interval 1.similar method to generate such surface,and compare the average ACF with the given ACF.In this section,surfaces are generated with given ACF R (X ,Y )=exp(Ϫ√X 2+Y 2/b ),where b is the correlation dis-tance.We choose the correlation distance b to be 2and 10.We will compare the average ACF of generated pro-files with the given ACF.128×128-point ACF is used to generate surfaces.Therefore,for Hu and Tonder’s method,m =n =128,and M =N =256.That is,256×256-point surface is generated.But,for Newland’s and our methods,M =N =128.That is,128×128-point surfaces are generated.Fig.15show the average ACF for the surfaces gener-ated by Hu and Tonder’s method for correlation distance 10.Fig.16shows that for Newland’s method.Fig.17shows that for our method.It is found that our method is the best one.The average ACF of Newland’s method53J.-J.Wu /Tribology International 33(2000)47–58Fig.11.ACF’s of five 512-point profiles generated by Newland’s method with R (X )=exp(ϪX /2),sampling interval1.Fig.12.ACF’s of five 512-point profiles generated by Wu’s method with R (X )=exp(ϪX /50),sampling interval 1.is slightly lower than the desired ACF.Hu and Tonder’s method is the worst one.As for the correlation distance 2,the average ACF of surfaces generated by these three methods is very close to the desired ACF,which is shown in Fig.18.256×256-point surfaces for Newland’s and our methods are also generated.Fig.19shows the average ACF of profiles for the surface generated by Newland’s method.Fig.20shows that by our method.It is found that the average ACF’s are much better than those of 128×128-point surfaces,and our method is slightly better than Newland’s.From these examples,it is found that that,Newland’s and our methods are better than Hu and Tonder’s method.Our method is slightly better than Newland’s.For a very small correlation distance,all three methods areacceptable.Fig.13.ACF’s of five 512-point profiles generated by Wu’s method with R (X )=exp(ϪX /10),sampling interval1.Fig.14.ACF’s of five 512-point profiles generated by Wu’s method with R (X )=exp(ϪX /2),sampling interval 1.4.Discussion4.1.Profile generated from spectral densitySuppose {z r |r =0,1,2,%,ϱ}is the profile that we want to generate,and its FFT is Z k .But we can only generate a profile with finite length {z r |r =0,1,2,…,(N Ϫ1)},and itsFFT is Z˜k .Since,there is aliasing effect,the FFT has the following relationship.Z˜k ϭZ ˜N −k ϭZ k ϩZ N −k The spectral density of generated profile is S˜k ϭS ˜N −k ϭ|Z ˜k |2ϭ|Z ˜N −k |2ϭ|Z k ϩZ N −k |2If the amplitude at high frequency is small,54J.-J.Wu /Tribology International 33(2000)47–58Fig.15.ACF’s of five 256×256-point surfaces generated by Hu and Tonder’s method with R (X )=exp(Ϫ√X 2+Y 2/10)sampling interval1.Fig.16.ACF’s of five 128×128-point surfaces generated by Newl-and’s method with R (X )=exp(Ϫ√X 2+Y 2/10),sampling interval 1.Z k ϩZ N −k ϷZ kThen,the spectral density at low frequency can be used to generate the profile.We take Z˜k ϭZ ˜N −k ϷZ k S˜k ϭS ˜N −k ϭ|Z ˜k |2ϷS k Then,|Z ˜k |can be obtained from S k .The profile can be generated.Newland’s method is wrong in estimating spectral density.In Fig.3,the spectral density is S (w )=1/(w 2+1).At high frequency,S (w )ෂ1/w 2.Since S˜k ϭN −k N S k ϩk N S N −k the errorisFig.17.ACF’s of five 128×128-point surfaces generated by Wu’s method with R (X )=exp(Ϫ√X 2+Y 2/10),sampling interval1.Fig.18.ACF’s of a 256×256-point surface generated by Hu and Ton-der’s method and 128×128-point surface by Newland’s and Wu’s methods with R (X )=exp(Ϫ√X 2+Y 2/10),sampling interval 1.⌬ϭ1ϪS ˜k S k Ϸk N ͫ1Ϫw 2kw 2N −k ͬϭk N ͫ1Ϫͩk N −k ͪ2ͬfor k ϭ0,1,2,%,ͩN2Ϫ1ͪSo,the maximum occurs at k Ϸ0.33N .The error is about 25%.There is nearly 30%error in Fig.4.If the amplitude at high frequency is not small enough,|Z k |Ϸ|Z N −k |then the phases are important.Then,S ˜k ϭ|Z ˜k |2ϭ|Z k ϩZ N −k |2ϭS k ϩS N −k ϩ2S 1/2k S 1/2N −k cos(f k Ϫf N −k )SinceE {cos(f k Ϫf N −k )}ϭ055J.-J.Wu /Tribology International 33(2000)47–58Fig.19.ACF’s of five 256×256-point surfaces generated by Newl-and’s method with R (X )=exp(Ϫ√X 2+Y 2/10),sampling interval1.Fig.20.ACF’s of five 256×256-point surfaces generated by Wu’smethod with R (X )=exp(Ϫ√X 2+Y 2/10),sampling interval 1.It is reasonable to assume thatS˜k ϭS k ϩS N −k Therefore,if we consider the spectral density at high frequency,by using Eq.(2),the profile can be generated.At this case,if we consider the spectral density at high frequency,the error of Newland’s method for spectral density is S (w )=1/(w 2+1)⌬ϭ1ϪS ˜k S k +S N −k Ϸk w 2k +(N −k )w 2kN (w 2k +w 2N −k )for k ϭ0,1,2,%,ͩN2Ϫ1ͪSo,the maximum occurs at k Ϸ0.5N .The error is about 50%.From the above examples,we know how our method is correct and how Newland’s method is wrong.4.2.ACF and cir-ACFSuppose {z r |r =0,1,2,%,ϱ}is the profile that we want to generate.But we can only generate a profile with finite length {z r |r =0,1,2,…,(N Ϫ1)}.If the ACF of the {z r |r =0,1,2,%,ϱ}is R r ,the cir-ACF C r of the generated profile {z r |r =0,1,2,…,(N Ϫ1)}should be C r ϭN −r N R r ϩr NR N −rSo,C r R r as long as R r R N Ϫr .So,cir-ACF and ACF are different most of the times.However,it is quite different for a profile generated by Newland’s or our method.The ACF and the cir-ACF of the generated profile areR˜r ϭ1N −r͸N Ϫr Ϫ1k ϭ0z k z k ϩrC˜r ϭ1N ͸N Ϫ1k ϭ0z k z k ϩr where z k ϩr ϭz k ϩr −N for k ϩr ՆN .The mean ACF and mean cir-ACF are defined as R¯r ϭE [R ˜r ]C¯r ϭE [C ˜r ]The mean ACF R¯r and the mean cir-ACF C ¯r of the gen-erated profile should be the same (see Appendix A.2).R¯r ϭR ¯N −r ϭC ¯r ϭC ¯N −r Therefore,it is impossible to generate a profile withgiven ACF R r and cir-ACF C r ,simultaneously.If our method is used,the ACF R˜r used to generate the profile isR˜N −r ϭR ˜r ϭR r for r ϽN /2Then the mean ACF and cir-ACF of the generated profile will beR¯r ϭR ¯N −r ϭC ¯r ϭC ¯N −r ϭR r The generated profile has the correct mean ACF R¯r for r ϽN /2,but not for N /2Յr ϽN .The generated profile does not have the correct cir-ACF,either.If Newland’s method is used,the ACF R ˜r used to gen-erate the profile isR ˜N −r ϭR ˜r ϭC r ϭC N −r ϭN −r N R r ϩrNR N −rfor r ϽN /2The mean ACF and the mean cir-ACF of the generated profile will beR ¯r ϭR ¯N −r ϭC ¯r ϭC ¯N −r ϭN −r N R r ϩrNR N −rThe generated profile has the correct mean cir-ACF,butdoes not have the correct ACF.56J.-J.Wu/Tribology International33(2000)47–58So,no method meets the requirements of both R r and C r.It seems that our method is better;because,by using our method,the mean ACF is correct for rϽN/2.But, by using Newland’s method,no ACF is correct.From Figs.9–14,it is found that,in most cases,both Newland’s and our methods cannot generate profiles with correct ACF.It is because correct mean ACF does not guarantee that ACF of each profile is correct.There-fore,for one-dimensional profiles,our method is theor-etically correct,but not practical.From Figs.15–20,it is found that,the average ACF of all profiles for the surfaces generated by our method is very close to the real ACF.Newland’s method is slightly worse.Hu and Tonder’s method is the worst.4.3.Hu and Tonder’s methodIn Hu and Tonder’s method,a surface can be gener-ated by the following equation.z p,qϭ͸mϪ1rϭ0͸nϪ1sϭ0h r,s h p+r,q+spϭ0,1,2,%,MϪ1qϭ0,1,2,%,NϪ1whereE{h ij h kl}ϭͭ1if i=k,j=l0otherwiseIn fact,this equation is similar to Patir’s method[2], which is the MA(moving average)method.This is why Hu and Tonder said that their method is similar to MA method.In Patir’s method,h r,s can be obtained by solvingR p,qϭ͸mϪp rϭ0͸nϪq sϭ0h r,s h p+r,q+spϭ0,1,2,%,mϪ1(7) qϭ0,1,2,%,nϪ1This equation can be solved by an interactive technique, such as the Newton Method.However,Hu and Tonder used the technique of FFT to solve them.For an FFT,they should beR p,qϭ͸mϪ1rϭ0͸nϪ1sϭ0h r,s h p+r,q+spϭ0,1,2,%,mϪ1qϭ0,1,2,%,nϪ1and alsoh k+r,l+sϭh k+r−m,l+s−n for kϩrՆm,lϩsՆnThe h r,s solved by Hu and Tonder is not the solution of Eq.(7).We define⌬R p,q as the difference between the desired ACF and the ACF of surface generated by Hu and Tonder’s method.Then⌬R p,q is⌬R p,qϭ͸mϪ1rϭm͸nϪ1sϭn h r,s h p+r,q+spϭ0,1,2,%,NϪ1qϭ0,1,2,%,MϪ1If⌬R p,q is small enough,Hu and Tonder’s method could be a good approximation.Wefind that,if correlation distance is large or the number of points for the generated surface is small,⌬R is not negligible.We alsofind that,for a real random series{h r|r=1,%,n},E{h r h s}is not exactly0for r s, especially when n is not large enough.This may also cause the error.Therefore,only if correlation distance is very small and the number of points for the generated surface is very large,is Hu and Tonder’s method a good approximation for Eq.(7).Chilamankuri and Bhushan[8],Yu and Bhushan[10] and Poon and Bhushan[11]generated surfaces with vari-ous correlation distances.In their work,they found that the generated surfaces have the desired ACF.In fact,in most of their cases,their correlation distance is less than 0.22µm and their sampling interval is0.078µm.It is similar to a correlation distance less than3and sampling interval1of our cases.These cases are acceptable.But, some of their cases are not as good as they thought.4.4.SummaryNewland did not clarify cir-ACF and cir-spectral den-sity.Hu and Tonder did not clarify the periodic property of FFT.Therefore,both methods are incorrect.By using our method,a surface with given spectral density or ACF can be simulated as long as two-dimen-sional ACF or spectral density is known.If only one-dimensional ACF is known,then,two-dimensional ACF can be easily obtained byR r,sϭR k where kϭͱr2+s2Then,two-dimensional spectral density can be obtained from the one-dimensional ACF.If only one-dimensional spectral density S p(w)is known,then,Nayak’s[12]for-mula can be applied.S s(w)ϭ1p͵ϱw(w21Ϫw2)1/2ͭSЉp(w1)w1ϪSЈp(w1)w21ͮd w1This equation can be solved numerically.So,the two-dimensional spectral density can be obtained.Then,the surfaces can be obtained by Eq.(2).57J.-J.Wu /Tribology International 33(2000)47–58However,there is still limitation in our method.Although the mean ACF of surfaces generated by our method is correct,it does not mean that ACF of each profile is correct.Our method can generate surfaces with given spectral density,and can generate surfaces having correct average ACF.5.ConclusionFFT is a convenient tool in generating surfaces.How-ever,FFT cannot evaluate R r from R N Ϫr .Therefore,FFT has its restrictions.FFT cannot generate profiles with both correct ACF and correct cir-ACF.In this paper,a simulation method is proposed.In this method,FFT is used to generate surfaces with given spectral density or ACF.But,since FFT cannot differen-tiate ACF from cir-ACF,the mean ACF of surfaces gen-erated by our method is correct,but the mean cir-ACF is incorrect.Although this method cannot guarantee that the generated profiles having the correct ACF,it can guarantee that the average ACF of all profiles for the generated surfaces is correct.Appendix AA.1.Relation between FFT and convolutionWe use one-dimensional profile as an example.Hu and Tonder used the following equation to generate a profile.z r ϭ͸n Ϫ1s ϭ0h s h s +rͭr =0,1,2,%,(N −1)s =0,1,2,%,(n −1)(8)It is obvious that h should be {h r |r ϭ0,1,2,%,(N ϩn Ϫ1)}Thus,the length of {z }is N ,the length of {h }is n and the length of {h }is (N +n Ϫ1).All three variables have different length.Assume that h and H k have the following relationship.H k ϭ͸n Ϫ1r ϭ0h r exp ͩϪi 2p krn ͪh r ϭ1N ͸n Ϫ1k ϭ0H k exp ͩi 2p kr nͪ(9)k ,r ϭ0,1,2,%,(n Ϫ1)Therefore,Z k ϭH k A kwhere A k is the FFT of h .Since h ’s are random series,their spectral density S h is a constant K .Cir-ACF and spectral density of the generated profile should have the following properties.C r ϭ1N ͸n Ϫ1s ϭ0z s z s +r ϭ͸n Ϫ1k ϭ0h k h k +r s 2h where E {h i h j }ϭͭs 2h i =j 0i j S˜k ϭZ ∗Z ϭ|H k |2A ∗A ϭ|H k |2S h ϭK |H k |2where h k +r =h k +r ϪN for k +r Նn .In Hu and Tonder’s method,a new ACF is gener-ated by R˜n −r ϭR ˜r ϭR r r ϭ0,1,%,n /2Then,S ˜k is obtained by R ˜r ,then H k is obtained from S˜k ,S ˜k ϭ1n ͸n /2Ϫ1k ϭϪn /2R˜r exp ͩϪi 2p kr n ͪH k ϭͩS ˜k Kͪ1/2Then,by Eqs.(9)and (8),the generated profile will have ACF R r .R r ϭC r ϭ1n ͸n Ϫ1s ϭ0z s z s +r ϭ͸n Ϫ1k ϭ0h k h k +r s 2h(10)But,if we compute the ACF directly from the profile generated by Eq.(8),we find that R˜r ϭ1n −r͸n Ϫr Ϫ1k ϭ0z k z k ϩr ϭn −rn͸n Ϫr Ϫ1k ϭ0h k h k ϩr s 2h ϽR rC˜r ϭ1n ͸n Ϫ1k ϭ0z k z k ϩr ϭn −r n͸n Ϫr Ϫ1k ϭ0h k h k ϩr s 2hϩrn ͸n Ϫ1k ϭN Ϫrh k h k ϩr s 2h ϽC rwhich do not meet our requirement of ACF.For example,for a non-periodic eight-point profile,by using Eq.(8)with {h }of n ,we find that z 1ϭh 1h 1ϩh 2h 2ϩh 3h 3ϩh 4h 4z 2ϭh 1h 2ϩh 2h 3ϩh 3h 4ϩh 4h 5z 3ϭh 1h 3ϩh 2h 4ϩh 3h 5ϩh 4h 6z 4ϭh 1h 4ϩh 2h 5ϩh 3h 6ϩh 4h 7z 5ϭh 1h 5ϩh 2h 6ϩh 3h 7ϩh 4h 8z 6ϭh 1h 6ϩh 2h 7ϩh 3h 8ϩh 4h 9z 7ϭh 1h 7ϩh 2h 8ϩh 3h 9ϩh 4h 10z 8ϭh 1h 8ϩh 2h 9ϩh 3h 10ϩh 4h 11From Eq.(10),its ACF R 1and cir-ACF C 1should be R 1ϭC 1ϭ(h 1h 2ϩh 2h 3ϩh 3h 4ϩh 4h 1)s 2h。