07-velocity and Jacobian - 360文档中心

英文版大学物理-第七章

translational kinetic energy of the molecules of a
Chapter 7 The Kinetic Theory of Gases
7-1 Temperature and Thermal Equilibrium 7-2 Thermodynamic Variables and the
Equation of State 7-3 Pressure and Molecular Motion 7-4 The Microscopic Interpretation of
R = 8.31 J/molK. (the universal gas constant)
An alternative form of ideal gas law:
pV m RT N RT
M
NA
The Boltzmann’s constant
k R 1.38 1023 J / K
NA
pV NkT ,
We get 2 ( 1 mv2 ) kT ,
32
The average translational
K 1 mv2 3 kT kinetic energy of molecules in
2
2
an ideal gas directly depends
only on the temperature, not on
A sample of gas consists of many identical molecules. The molecules are very far apart in comparison to
their size; The direction of motion of any molecule is random; The molecules are treated as if they were hard

扩展线位移法快速实现机器人末端线速度规划

扩展线位移法快速实现机器人末端线速度规划陈琳;倪崇琦;戴骏;潘海鸿【摘要】将用于数控机床的扩展线位移方法引入到六自由度串联机器人末端运动线速度规划算法中,以提高机器人末端运动线速度控制效率.该算法借鉴数控机床中将刀具实际运动线速度与编程指定速度相互关联的扩展线位移方法思想,将机器人末端运动线速度转换为机器人运动矢量合角速度最终实现机器人末端运动线速度规划.通过搭建的机器人实验平台进行实验,实验结果表明提出的扩展线位移方法在机器人末端沿直线轨迹和圆弧轨迹运动时的实际合角速度与理论推算合角速度基本一致,实际关节角速度与理论关节角速度趋于一致.表明该方法能够有效实现机器人末端线速度轨迹规划.【期刊名称】《机械设计与制造》【年(卷),期】2016(000)009【总页数】4页(P176-179)【关键词】六自由度串联机器人;末端运动线速度规划;扩展线位移法;矢量合角速度【作者】陈琳;倪崇琦;戴骏;潘海鸿【作者单位】广西大学机械工程学院,广西南宁530004;广西制造系统与先进制造技术重点实验室,广西南宁530004;广西大学机械工程学院,广西南宁530004;广西大学机械工程学院,广西南宁530004;广西大学机械工程学院,广西南宁530004;广西制造系统与先进制造技术重点实验室,广西南宁530004【正文语种】中文【中图分类】TH16;TP242.2目前工业机器人已经广泛应用在许多工业领域中，在诸如焊接、喷涂等作业场合，不仅对机器人末端位置和姿态有一定要求，对工作时速度、加速度都有着严格要求，以满足实际生产需要。

针对机器人速度分析规划方法主要包括CAD变量几何法及雅可比矩阵方法。

文献［1］利用CAD变量几何方法求解平面多自由度机器人臂速度和加速度，但求解精度不高。

雅可比矩阵方法在机器人运动速度规划中具有十分重要的作用，文献［2］针对机器人任意点速度推导出其雅可比矩阵计算公式，并编制出相应程序；文献［3］基于螺旋理论建立少自由度操作空间与关节空间的速度一一映射关系，根据雅可比矩阵分析少自由度机器人的奇异性；文献［4］采用雅可比矩阵方法研究6R串联机器人关节空间速度和操作空间运动速度之间的映射关系；文献［5］采用速度雅可比矩阵在已知各关节位置及速度前提下对机器人进行正、逆速度分析；文献［6］利用全局雅可比矩阵对四足变结构机器人进行速度分解控制。

雅可比矩阵

v d 1 d x lim D lim xDt w Dt 0 Dt Dt 0 D lim J (q)qDt J (q)dq
Dt 0
J11 J 21 v J 31 w J 41 J 51 J 61
教材例题2.1：逆雅可比矩阵的示例：例2.1 如图2.2所示的二自由度机械手，手部沿固定坐标系X0轴正向以1.0 m/s的速度移动，杆长l1=l2=0.5 m。设在某瞬时θ1=30°， θ2=60°，求相应瞬时的关节速度。
解由式(2.6)知，二自由度机械手速度雅可比为
因此，逆雅可比为
2.1.3 机器人雅可比讨论机器人的奇异形位分为两类： (1) 边界奇异形位：当机器人臂全部伸展开或全部折回时，使手部处于机器人工作空间的边界上或边界附近，出现逆雅可比奇异，机器人运动受到物理结构的约束。这时相应的机器人形位叫做边界奇异形位。 (2) 内部奇异形位：两个或两个以上关节轴线重合时，机器人各关节运动相互抵消，不产生操作运动。这时相应的机器人形位叫做内部奇异形位。
对力雅可比矩阵的补充说明：
虚功方程力雅可比分析：
2.2.3 机器人静力计算
机器人操作臂静力计算可分为两类问题： (1) 已知外界环境对机器人手部的作用力F，(即手部端点力 F-F′)，利用式(2.20)求相应的满足静力平衡条件的关节驱动力矩τ。 (2) 已知关节驱动力矩τ，确定机器人手部对外界环境的作用力或负载的质量。第二类问题是第一类问题的逆解。逆解的关系式为
或写成
根据虚位移原理，机器人处于平衡状态的充分必要条件是对任意符合几何约束的虚位移有δW=0，并注意到虚位移δq和δX之间符合杆件的几何约束条件。利用式δX=Jδq，将式(2.18)写成

The University of Florida Sparse Matrix Collection

The University of Florida Sparse Matrix CollectionTimothy A.Davis1andYifan Hu2The University of Florida Sparse Matrix Collection is a large,widely available,and actively grow-ing set of sparse matrices that arise in real applications.Its matrices cover a wide spectrum of do-mains,include those arising from problems with underlying2D or3D geometry(such as structural engineering,computationalﬂuid dynamics,model reduction,electromagnetics,semiconductor de-vices,thermodynamics,materials,acoustics,computer graphics/vision,robotics/kinematics,and other discretizations)and those that typically do not have such geometry(such as optimization, circuit simulation,economic andﬁnancial modeling,theoretical and quantum chemistry,chemical process simulation,mathematics and statistics,power networks,and other networks and graphs). The collection is widely used by the sparse matrix algorithms community for the development and performance evaluation of sparse matrix algorithms.The collection includes software for accessing and managing the collection,from MATLAB,Fortran,and C,as well as online search capabil-ity.Graph visualization of the matrices is provided,and a new multilevel coarsening scheme is proposed to facilitate this task.Categories and Subject Descriptors:G.1.3[Numerical Analysis]:Numerical Linear Algebra—linear systems(direct methods),sparse and very large systems;G.4[Mathematics of Com-puting]:Mathematical Software—algorithm analysis,eﬃciency;G.2[Discrete Mathematics]: Graph TheoryGeneral Terms:Algorithms,Experimentation,PerformanceAdditional Key Words and Phrases:sparse matrices,performance evaluation,graph drawing, multilevel algorithms1.INTRODUCTIONAlthough James Wilkinson’s foundational work in numerical analysis touched only lightly upon sparse matrix computations,his deﬁnition of a sparse matrix is widely used([Gilbert et al.1992],for example).Wilkinson deﬁned a matrix as“sparse”if it has enough zeros that it pays to take advantage of them.He actually stated his deﬁnition in the negation:1Dept.of Computer and Information Science and Engineering,Univ.of Florida,Gainesville,FL, USA.davis@cise.uﬂ.edu.http://www.cise.uﬂ.edu/∼davis.Portions of this work were supported by the National Science Foundation,under grants9111263,9223088,9504974,9803599,0203720, and0620286.2AT&T Labs Research,Florham Park,NJ,USA.yifanhu@.Permission to make digital/hard copy of all or part of this material without fee for personal or classroom use provided that the copies are not made or distributed for proﬁt or commercial advantage,the ACM copyright/server notice,the title of the publication,and its date appear,and notice is given that copying is by permission of the ACM,Inc.To copy otherwise,to republish, to post on servers,or to redistribute to lists requires prior speciﬁc permission and/or a fee.c 20YY ACM0098-3500/20YY/1200-0001$5.00ACM Transactions on Mathematical Software,Vol.V,No.N,M20YY,Pages1–0??.2·T.A.Davis and Y.HuThe matrix may be sparse,either with the non-zero elements concen-trated on a narrow band centered on the diagonal or alternatively theymay be distributed in a less systematic manner.We shall refer to amatrix as dense if the percentage of zero elements or its distributionis such as to make it uneconomic to take advantage of their presence.[Wilkinson1971]In other words,if you can save time or memory(usually both)by exploiting the zeros,then the matrix is sparse.An interesting aspect of this deﬁnition is that it is dependent not only on the matrix,but on the algorithm used on the matrix.For example,a matrix may be“sparse”for an iterative method for solving linear systems or for a graph theoretic algorithm,but not for a sparse factorization method.This article describes the University of Florida Sparse Matrix Collection(here-after referred to as the UF Collection),which contains sparse matrices arising in a wide range of applications.Section2gives the motivation for collecting sparse matrices from real applications and making them widely available.Section3de-scribes the current state of the collection and the breadth of problems the matrices represent.Section4describes the algorithm used,and new techniques developed, for visualizing the matrices.Section5describes the three data formats used for storing the matrices,and the kinds of auxiliary information available for each ma-trix.Section6describes the four methods for searching and downloading matrices of interest:a MATLAB interface(UFget),a Java interface(UFgui),the matrix web pages(via a standard browser and a web-based search tool),and Amazon Web Services TM.Examples of how the UF Collection can be used for performance eval-uation are given in Section7.The future of the UF Collection depends critically upon the submission of new matrices,as discussed in Section8.Finally,Section9 summarizes the contributions of the UF Collection and its associated software.In this paper,a graph or mesh is said to have2D or3D geometry if its vertices have a position on an xy or xyz plane that naturally arises from the problem being solved.A sparse matrix is said to have2D or3D geometry if its nonzero pattern is the adjacency matrix of such a graph.The notation|A|refers to the number of nonzeros in a matrix.2.MOTIVATIONThe role of sparse matrices from real applications in the development,testing, and performance evaluation of sparse matrix algorithms has long been recognized. Theﬁrst established collection was started by Duﬀand Reid(1970to1979,[Duﬀand Reid1979]),and then compiled into the Harwell-Boeing collection by Duﬀ, Grimes,and Lewis[Duﬀet al.1989].This collection provided the starting point of University of Florida Sparse Matrix Collection.Since the start of our collection in 1992,additional matrices have been added over the years,and other collections have been made,many of which have also been incorporated into the UF Collection(such as[Bai et al.1996;2008;Batagelj and Mrvar2008;Boisvert et al.2008;Boisvert et al.1997;Dumas2008;Gay2008;Gould et al.2008;Koster2008;M´e sz´a ros2008; Mittelmann2008;Resende et al.1995;Rudnyi2008;Rudnyi et al.2006;Saad2008; Schenk2008]).ACM Transactions on Mathematical Software,Vol.V,No.N,M20YY.University of Florida sparse matrix collection·3 The Matrix Market[Boisvert et al.1997]is the most similar collection to the UF Collection.Both collections include a search tool,and both categorize the matrices by application domain and problem source.Both provide matrices in similarﬁle formats.Both provide a web page for each matrix,with basic statistics andﬁgures. They diﬀer in size,with the UF Collection containing much larger matrices and 4.5times as many matrices.The latest matrix added to the Matrix Market was in 2000,whereas UF collection is constantly being updated with new matrices.The largest matrix in the Matrix Market has dimension90,449with2.5million nonzeros, whereas the largest matrix in the UF Collection has a dimension of28million with 760million nonzeros.Nearly every matrix in the Matrix Market is also included in the UF Collection.However,the Matrix Market does include matrix generators; the UF Collection has no matrix generators.Nearly all research articles that include a section on the performance analysis of a sparse matrix algorithm include results on matrices from real applications or parametrized matrices that mimic those that could arise in practice.Since maintaining a large collection of matrices from real applications is not trivial,an alternative isﬁrst considered,namely,parametrized and random matrix generators.2.1Parameterized and random matricesRandomized or repeatable parametrized sequences of matrices can be easily con-structed via simple matrix generators.Examples are listed below,as a comparison and contrast with the UF Collection,which does not include any matrix generators.(1)Simple discretizations of the Laplace operator on square2D and3D meshes.(2)The L-shaped meshes of[George and Liu1981].(3)Least-squares problems from square2Dﬁnite-element meshes[George et al.1983].(4)The LAPACK test matrix generators[Anderson et al.1999],which can generatebanded matrices and sparse matrices with random nonzero pattern.These also appear in the Matrix Market.(5)The Non-Hermitian Eigenvalue Problem(NEP)matrix generators[Bai et al.1996;2008].These include regular2D meshes and random patterns;additional matrices are made available only as operators(y=Ax,where A is not explicitly represented).(6)Zlatev’s matrix generators,which create a sparse Toeplitz structure(whereselected diagonals are all nonzero)[Zlatev1991].(7)Higham’s Matrix Computational Toolbox[Higham2002],part of which appearsas the gallery function in MATLAB.The gallery contains three functions that generate parametrized sparse matrices from regular2D meshes(neumann, poisson,and wathen)and two that generate banded matrices(toeppen and tridiag).No randomization is used,except in the values(but not pattern) of the wathen matrices.Many of these matrix generators also appear in the Matrix Market.(8)Matrices with purely random nonzero patterns[Erd¨o s and R´e nyi1959;Gilbert1959].They can be generated by sprand,sprand,and sprandsym in MATLAB [Gilbert et al.1992].ACM Transactions on Mathematical Software,Vol.V,No.N,M20YY.4·T.A.Davis and Y.Hu(9)The YM11subroutine in HSL can generate random sparse matrices,with op-tions for ensuring structural nonsingularity and bandedness(a matrix is struc-turally nonsingular if there exists a permutation so that the matrix has a zero-free diagonal)[Duﬀ2001].(10)In between the purely-random and purely-parametrized classes of matricesare partially randomized matrices with speciﬁc structure,as exempliﬁed by the CONTEST toolbox[Taylor and Higham2009].They provide a set of graph generators,some of which use an underlying2D or3D geometry and some that do not.For example,they include an implementation of the small-world graphs of[Kleinberg2000],which are2D meshes with additional randomized edges to distant vertices in the work generators in CONTEST that do not have geometry include the scale-free graph models of[Barab´a si and Albert1999].The prime advantage of random and parametrized matrices is that they are very easy to generate in any size desired.Matrices from real applications are very diﬃcult to generate and it is hard to vary them in size.Another key advantage of random and parametrized matrices is that asymptotic results can sometimes be derived.For example,the nested dissection ordering applied to a2D s-by-s mesh leads to an asymptotically optimal ordering for sparse Cholesky factorization,with31(n log2n)/8+O(n)nonzeros in the Cholesky factor L,and requiring829(n3/2)/84+O(n log n)ﬂoating point operations to compute, where n=s2is the dimension of the matrix[George and Liu1981].Purely random matrices may be useful for testing some sparse matrix applica-tions.For example,they can be useful for testing the convergence of iterative methods,assuming the generator has some control over the eigenvalue spectra(as is the case for sprandsym,for example).However,under Wilkinson’s deﬁnition they are not truly sparse when factorized by direct methods.With modest assumptions, purely random n-by-n matrices with O(n)nonzero entries require O(n3)time and O(n2)memory to factorize,because of catastrophicﬁll-in[Duﬀ1974].Catastrophic ﬁll-in very rarely occurs in matrices arising in real applications,and even when it does it is an indication that direct factorization methods are not applicable to that problem.Thus,performance obtained on purely random matrices will not indicate how well a sparse matrix factorization method will work on matrices from real ap-plications.Purely random networks also do not accurately model the characteristics of real networks[Watts and Strogatz1998].In between the two extremes of highly structured matrices(banded,or square 2D and3D meshes,for example)and purely randomized matrices is another class of matrix generators that can create matrices with some regular structure and a carefully selected random structure.The small-world graphs from the CONTEST toolbox are one such example.These are constructed to capture essential properties of large networks,such as the small-world property,scale-free degree distribution, motifs,and graphlet frequency[Taylor and Higham2009].These models are very useful for network algorithms,but they have not yet been shown to mimic the per-formance of sparse matrix factorization methods on matrices from real applications. ACM Transactions on Mathematical Software,Vol.V,No.N,M20YY.University of Florida sparse matrix collection·5Fig.1.Sparse matrix pattern of a RAH-66Comanche helicopter2.2Sparse matrices from real applicationsWhile useful,random and parametrized matrices have their limitations.This mo-tivates the development of the UF Collection,which focuses on matrices from real applications1.One of the matrices in the UF Collection is a complete representation of a Boe-ing/Sikorsky RAH-66Comanche helicopter,with3D geometry(obtained from Alex Pothen).Figure1shows the sparsity pattern of the matrix.Nested dissection is still useful for this problem,but asymptotic results cannot be derived for a complex structure such as this one.A simple square3D mesh would not be a good approxi-mation to this graph,since it would not capture the properties of the long helicopter rotors,or the topology of the hole where the rear rotor is located.This matrix is one of the few in the collection where the3D coordinates of the vertices are given;a picture of the graph drawing when these coordinates are used is shown in Figure2, on the left.For comparison,a force-directed graph drawing[Hu2005],which recre-ates these coordinates based only on the connectivity of the graph,is shown on the right.The method for creating thisﬁgure is discussed in Section4.The structure is warped,with the tail rotor twisted out of its housing(to the right)and itsﬁve main rotors shrunken(on the top),but it is still clear that the graph represents some kind of3D problem.Edge colors in the force-directed graphs represent the distance between two vertices.Figure3is the nonzero pattern of matrix arising from circuit simulation(left)and its force-directed graph drawing(right).The graphs of the Comanche helicopter and this electronic circuit indicate that one is a3D problem and the other is a network with no2D or3D geometry.This diﬀerence is not clear by merely comparing the nonzero patterns of their matrices.Parameterized matrices can capture many key properties of matrices from real applications,such as the regularity of a3D mesh,or the small-world properties of a network,but they cannot capture the rich and complex structure of matrices such1The UF Collection does include a handful of random matrices,which remain in the collection for historical reasons.ACM Transactions on Mathematical Software,Vol.V,No.N,M20YY.6·T.A.Davis and Y.HuFig.2.Graph of a RAH-66Comanche helicopter,using given3D coordinates(left)and its force-directed graph drawing(right)Fig.3.Sparse matrix pattern of an electronic circuit from Steve Hamm,Motorola(left)and its force-directed graph(right)as those presented above.2.3Collecting the matricesThe ideal matrix collection would be informally representative of matrices that arise in practice.It should cast a broad net so as to capture matrices from every applica-tion domain that relies on sparse matrix methods.For example,matrices arising in circuit simulation(and other network-related domains)diﬀer greatly from matrices arising from the discretization of2D and3D physical domains;this can be clearly seen in the precedingﬁputationalﬂuid dynamics matrices diﬀer from structural engineering matrices,and both are vastly diﬀerent from matrices arising in linear programming orﬁnancial portfolio optimization.The collection should be kept up to date,since matrices of interest grow in size each year as computer mem-ories get larger.New application domains also appear,such as eigenvalue problems arising in web connectivity matrices[Kamvar2008;Kamvar et al.2004;Page et al. 1998],which have existed only since the mid1990’s.ACM Transactions on Mathematical Software,Vol.V,No.N,M20YY.University of Florida sparse matrix collection·7 Sparse matrix algorithm developers use the matrices in the UF Collection to develop their methods,since theoretical asymptotic time/memory complexity anal-ysis only goes so far.If there are no matrices available to developers from a given application domain,it is quite possible that when their methods are used in that domain,the performance results will be disappointing.This provides a strong mo-tivation for computational scientists to submit their matrices to a widely available collection such as this one,so that gaps can be avoided.Thus,new application areas are always being added to the collection.Our strategy for adding matrices to the collection is simple,although admittedly ad hoc.Theﬁrst author maintains a collection of codes for the direct solution of sparse linear systems.End-users of this software are uniformly requested to submit matrices to the collection.Additional matrices are requested when they are found cited in articles and conference presentations,which includes a wider range of matrices(such as graphs arising in network analysis).Any matrices received are included,unless they are clear duplications of matrices already in the collection. Small matrices are not included,unless they are subsets of a larger collection(such as the Pajek data set).This strategy does introduce an unavoidable source of bias,but we have at-tempted to avoid this bias by relying on matrices collected by others.The UF Collection includes many sets of matrices collected in this way,such as those col-lected by Saad,who develops iterative methods for sparse linear systems[Saad 2003].Not all matrices in the collection arise from a sparse linear system.For example,many linear programming problems have been included in the UF Collec-tion[Gay2008;M´e sz´a ros2008;Mittelmann2008;Resende et al.1995].Network and combinatorial problems are also included,such as a matrix of prime numbers from Problem7of Trefethen’s100-digit challenge[Trefethen2002].Once a matrix is collected,its maintenance is the responsibility of theﬁrst author. This guards against arbitrary modiﬁcations in the matrices,which can occur when the matrix submitters modify their matrix problems.Repeatability of experiments based on the UF Collection is thus ensured.3.DESCRIPTION OF THE COLLECTIONAs of March2010the UF Sparse Matrix Collection consists of2272problems,each of which contains at least one sparse matrix(typically just one).It often represents a sparse linear system of the form Ax=b,where b is sometimes provided as well as the sparse matrix A.Many other problems are eigenvalue problems,and many matrices come from unstructured graphs and networks with no associated linear system.In some cases,a problem consists of a sequence of closely related matrices, such as Jacobian matrices arising in the solution of a nonlinear system.Some problems include two sparse matrices,such as a stiﬀness matrix and a mass matrix from a structural engineering eigenvalue problem.In this case,A is the stiﬀness matrix,and the mass matrix appears as an auxiliary matrix in the problem.3.1Application areasThe collection is divided into157diﬀerent matrix groups,with more groups added as new matrices are submitted to the collection.A complete list of these groups is too long to include here;details are given on the collection’s web site.Each groupACM Transactions on Mathematical Software,Vol.V,No.N,M20YY.8·T.A.Davis and Y.HuTable I.Partial list of sources of matricesNon-Hermitian Eigenvalue Problems[Bai et al.1996]Pajek Networks[Batagelj and Mrvar2008]Multistage stochasticﬁnancial modeling[Berger et al.1995]The Matrix Market(collection)[Boisvert et al.1997]Univ.of Utrecht circuit simulation[Bomhof and van der Vorst2000]MRI reconstruction M.Bydder,UCSDHarwell-Boeing Collection[Duﬀet al.1989]Combinatorial problems[Dumas2008]Frequency domain,nonlinear analog circuits[Feldmann et al.1996]NETLIB Linear Programming Test Problems[Gay2008]Symmetric sparse matrix benchmarks[Gould et al.2008]Linear programming problems[Gupta1996]Stanford/Berkeley Web Matrices[Kamvar et al.2004]Xyce circuit simulation matrices[Keiter et al.2003]Computer vision problems[Kemelmacher2005]PARASOL Matrix Collection[Koster2008]Linear programming test set[M´e sz´a ros2008]2D and3D semiconductor physics[Miller and Wang1991]Linear programming test set[Mittelmann2008]QAPLIB,quadratic assignment[Resende et al.1995]Oberwolfach Model Reduction Benchmarks[Rudnyi et al.2006]SPARSKIT matrix collection[Saad2008]Univ.of Basel Collection[Schenk2008]DNA electrophoresis matrices[van Heukelum et al.2002]Chemical engineering problems[Zitney1992;Zitney et al.1996]typically consists of a set of matrices from a single source.A few of the matrix groups in the UF Collection consist of entire sets of matrices from another sparse matrix collection.In this case,the group may consist of matrices from very diﬀerent application areas.For example,the Harwell-Boeing collection forms a single group [Duﬀand Reid1979;Duﬀet al.1989].Sources(papers and web sites)for some of the matrix groups are listed in Table I,which gives aﬂavor of the range of problems the collection contains.The group of a problem forms part of the full name of the problem.For exam-ple,the full name of the west0479problem in the Harwell-Boeing collection(the HB group in the UF sparse matrix collection)is HB/west0479.In addition,all problems are tagged with a string,kind,which indicates the application domain of the prob-lem.The group and kind of a problem are not the same.Some groups have many diﬀerent kinds of problems,and problems of the same kind can appear in diﬀerent groups.For example,the HB/west0479problem is in the HB group,and its kind is tagged as a“chemical process simulation problem.”Five other groups include chemical process simulation problems,namely,the Bai,Grund,Mallya,VanVelzen, and Zitney groups,all of which are named after the person from whom the matrices were obtained.A complete list of the kind strings for all problems is given in Table II.The table is split into two categories:problems with no underlying geometry and problems with2D or3D geometry.The collection contains68matrices with random nonzero pattern.They appear in the collection only because they already occur as widely-used test problems in ACM Transactions on Mathematical Software,Vol.V,No.N,M20YY.University of Florida sparse matrix collection·9 Table II.Summary of Problem.kind for all2272problems1516problems with no2D/3D geometry70chemical process simulation problem251circuit simulation problem299combinatorial problem11counter-example problem68economic problem4frequency-domain circuit simulation problem23least squares problem342linear programming problem135optimization problem56power network problem10statistical/mathematical problem61theoretical/quantum chemistry problem88directed graph8bipartite graph23undirected graph68random graph756problems with2D/3D geometry13acoustics problem166computationalﬂuid dynamics problem12computer graphics/vision problem44electromagnetics problem28materials problem42model reduction problem3robotics problem35semiconductor device problem288structural problem31thermal problem942D/3D problem(other than those listed above)another collection that was subsequently added to the UF sparse matrix collection. In retrospect,their inclusion in the UF Collection was a mistake,but matrices are never removed once added,since this would cause confusion when replicating results that use the collection.We guarantee statements such as“all matrices with property X in the collection as of date Y”always give a single consistent answer.3.2Matrix statisticsThe2272matrices in the collection come from359diﬀerent authors and50diﬀerent editors.A matrix editor is the one who collected the matrix for inclusion into any established collection(not just the UF Collection2,but also for others such as the Harwell/Boeing collection3,the Matrix Market4,or the GRID/TSLE collection5).A matrix author is the one who created the matrix.Each matrix has an editor and an author(sometimes the same person).Editors contributing at least10matrices are listed in Table III(sorted by the number of matrices collected).Figures4and5plot the matrix size(dimension and number of nonzeros)versus 2UF Sparse Matrix Collection:http://www.cise.uﬂ.edu/research/sparse/matrices3Harwell/Boeing collection:/nag/hb/hb.shtml4Matrix Market:/MatrixMarket5GRID/TLSE collection:http://gridtlse.enseeiht.fr:8080/websolve/ACM Transactions on Mathematical Software,Vol.V,No.N,M20YY.10·T.A.Davis and Y.HuFig.4.Matrix dimension(the largest of row/column dimension if rectangular)versus year created. The solid line is the cumulative sum.Fig.5.Number of nonzeros in each matrix versus year created.The solid line is the cumulative sum.ACM Transactions on Mathematical Software,Vol.V,No.N,M20YY.University of Florida sparse matrix collection·11 Table III.Primary editors of the UF Collection#matrices editor720T.Davis293J.-G.Dumas255I.Duﬀ,R.Grimes,J.Lewis166 C.Meszaros103O.Schenk78Z.Bai,D.Day,J.Demmel,J.Dongarra72V.Batagelj67Y.Ye61 A.Baggag,Y.Saad48 D.Gay43R.Fourer39N.Gould,Y.Hu,J.Scott38 E.Rudnyi35G.Kumfert,A.Pothen34 A.Curtis,I.Duﬀ,J.Reid28J.Chinneck27N.Gould25H.Mittelmann23R.Boisvert,R.Pozo,K.Remington,ler,R.Lipman,R.Barrett,J.Dongarra23 F.Grund22J.Koster13I.Lustig12H.Simon11M.ResendeFig.6.Overall histograms of matrix dimensions and nonzerosACM Transactions on Mathematical Software,Vol.V,No.N,M20YY.12·T.A.Davis and Y.Huthe year in which the matrices were created.The solid line in theﬁgures is the cumulative sum of the data plotted.Bothﬁgures show an exponential growth in problem sizes,similar to how computer memory sizes have grown since1970. Note that small problems are still being added to the collection.This is because a matrix group often includes a range of related problems,from small test cases to the largest problems of interest.The outlier matrix in1971is from the Edinburgh Associative Thesaurus,located at and obtained from the Pajek data set[Batagelj and Mrvar 2008].It is a graph with23,219vertices and325,592edges that wasﬁrst constructed in1971[Kiss et al.1973].Figure6plots two histograms of the overall distribution of matrices in the col-lection.4.VISUALIZING THE COLLECTIONMany basic facts,including symmetry,structural rank,ordering statistics,as well as a plot of the sparsity pattern,are given for each matrix in the collection.However these facts alone do not always give suﬃcient information about the matrix.For example,does this matrix come from an application involving2D or3D mesh?Or from a small-world network?Are there other structures in the applications that are not discernible from a plot of the sparsity pattern?To help answer these questions,a visualization of each matrix in the form of graph drawing is provided.If the matrix is structurally symmetric,it is taken as the adjacency matrix of an undirected graph,where two vertices i and j are connected if the(i,j)-th entry of the matrix is nonzero.Rectangular or structurally unsymmetric matrices are treated as bipartite graphs.More speciﬁcally,the augmented matrix0A A T0is used as the adjacency matrix of an undirected bipartite graph whose vertices are composed of the rows and columns of the matrix.We provide two graphs for square matrices with unsymmetric structure:the graph of A+A T and the augmented matrix above.The graph drawings depend only on the nonzero pattern,not the values.The basic algorithm used for drawing the graphs is a multilevel force-directed algorithm[Hu2005].However this algorithm fails to produce informative drawings for some matrices.We propose here a new coarsening scheme to deal with these cases.In the following we brieﬂy describe the basic algorithm,followed by the new coarsening scheme.4.1The graph drawing algorithmThe basic algorithm[Hu2005]employs a force-directed model[Fruchterman and Reingold1991].Vertices are treated as particles with electrical charges that push each other away.This repulsive force F r(i,j)between any two vertices i and j is proportional to the inverse of their distance,F r(i,j)=K2x i−x j ,i=j.ACM Transactions on Mathematical Software,Vol.V,No.N,M20YY.。

Motion control of robot manipulators

Figure 1: Puma Robot Manipulator Robot manipulators are basically multi{degree{of{freedom positioning devices. The robot, as the \plant to be controlled", is a multi{input/multi{output, highly coupled, nonlinear mechatronic system. The main challenges in the motion control problem are the complexity of the dynamics, and uncertainties, both parametric and dynamic. Parametric uncertainties arise from imprecise 1
Motion Control of Robot Manipulators
Mark W. Spong The Coordinated Science Laboratory, University of Illinois at Urbana{Champaign, 1308 W. Main St., Urbana, Ill. 61801 USA.
1.2 Kinematics
Kinematics refers to the geometric relationship between the motion of the robot in Joint Space and the motion of the end{e ector in Task Space without consideration of the forces that produce the motion. The Forward Kinematics Problem is to determine the mapping

利用拉格朗日方法对三连杆机械臂动力学建模与并进行控制仿真

Dynamic Modelling and Control Simulation of aThree-link Robotic Manipulator1.Problem to solveA three-link manipulator is shown as following figure.The geometric parameters of the robot are specified in following figure.L3m3m2m1The links of the robot are slender rods with uniformed masses which are m1,m2and m3repectively. Based the abovementioned information,please:(1)Find the solution to the forward kinematics and inverse kinematic for this robot and derive the Jacobian transformation which transfers the joint rates to the tip velocity.(2)Derive the dynamic model for this robot applying Lagrange formulation and verify the correctness of the dynamic model.(3)Design a controller for this robot according to model-based patitioning method or other approaches based on the established dynamic model to enable the robot have the ability to follow given trajectories.(4)Verify the correctness of the above derivation applying the code in robotics toolbox.(5)Conduct the control simulation for this robot to following a given tip trojectory in Cartesian space. (optional)目录1．求解该机器人正运动学和逆运动学解，并推导雅可比变换矩阵： (3)1.1正运动学解 (3)1.2逆运动学求解 (3)1.3关节角速率传递到末端效应器速度的雅可比矩阵变换 (4)2．应用拉格朗日公式推导该机器人的动力学模型： (6)3、机器人轨迹跟踪控制器设置 (10)4、前三题验证 (14)4.1、正运动学验证 (14)4.2、逆运动学验证 (14)4.3、雅可比变换矩阵验证 (14)4.4、动力学方程验证 (14)5、笛卡尔空间轨迹跟踪 (15)sdh_kinematics.m (19)sdh_jacobian_mat.m (19)sdh_inverse_kinematics.m (19)sdh_rvctool_inverse_kinematics.m (20)mdh_dynamics.m (21)M_theta.m (22)V_G_theta.m (23)forward_dynamics.m (25)Car_M_theta.m (28)Car_V_G_theta.m (30)trans_J.m (33)Kine_theta.m (34)J_theta.m (36)1．求解该机器人正运动学和逆运动学解，并推导雅可比变换矩阵：1.1正运动学解DH 参数表i 1i -α1a i -1d i -1i θ-1001θ2o901L 02θ302L 03θ43L 0连杆坐标变换矩阵：111101000000100001c s s c T θθθθ-⎛⎫⎪ ⎪= ⎪⎪⎝⎭2211222000-10000001c s L T s c θθθθ-⎛⎫ ⎪ ⎪= ⎪⎪⎝⎭33233230000010001c s L s c T θθθθ-⎛⎫⎪ ⎪= ⎪⎪⎝⎭334100010000100001L T ⎛⎫ ⎪ ⎪= ⎪ ⎪⎝⎭从操作臂末端到基座的坐标变换矩阵为：0012341234T T T T T ==12312313123212111231231312321211232332322---001c c c s s L c c L c c L c s c s s c L s c L s c L s s c L s L s ++⎛⎫⎪++ ⎪⎪+ ⎪⎝⎭1.2逆运动学求解假设：1231112131421222324012312343132333441424344r r r r r r r r T T T T r r r r r r r r θθθ⎛⎫ ⎪ ⎪= ⎪ ⎪⎝⎭（）（）（）两侧同时左乘()11T -：111100-0000100001c s s c ⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭11121314212223243132333441424344r r r r r r r r r r r r r r r r ⎛⎫ ⎪ ⎪ ⎪⎪⎝⎭=123234T T T =2323122323232322323000100001c s L L c L c s c L s L s -++⎛⎫⎪- ⎪ ⎪+ ⎪⎝⎭使左右两边（2,4）相等：110x y s p c P +=令11cos x P ρϕ=，11sin y P ρϕ=，代入上式：11111sin()0ρϕθθϕ-=⇒=或者1πϕ+对（1）式左乘()112T -：221222120--00-10001c s L c s c L s ⎛⎫⎪ ⎪ ⎪ ⎪⎝⎭111100-0000100001c s s c ⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭11121314212223243132333441424344r r r r r r r r r r r r r r r r ⎛⎫ ⎪⎪ ⎪ ⎪⎝⎭=2334T T ⇒1212212121221211----0001c c s c s L c c s s s c L s s c ⎛⎫⎪ ⎪ ⎪ ⎪⎝⎭11121314212223243132333441424344r r r r r r r r r r r r r r r r ⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭=332333333-000010001c s L L c sc L s +⎛⎫⎪⎪⎪ ⎪⎝⎭使左右两边（1，4）和（2，4）分别相等得：()()12122122331212212331,42,4-x y z x y z c c P s c P s P L c L L c c s P s c P c P L s L s ++-=+⎧⎪⎨-++=⎪⎩：：两式平方求和可得：22222222111111111232332222x y z x y x y c p s p p L p p c s L p c L p s L L L L c ++++--=++2222222212311111111312--2222x y z x y x y L L L c p s p p p p c s L p c L p s c KL L ++++--∴==解得：033arccos -180180o K θθ⎡⎤=±∈⎣⎦（，）再根据（2，4）：22233sin()L s ρθϕ+=33222arcsinL s θϕρ⇒=-解得：或者3322-arcsinL s πϕρ-1.3关节角速率传递到末端效应器速度的雅可比矩阵变换各连杆的速度和角速度在相应的坐标系下表示为：0000ω⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦0000v ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦11011001110ˆ0R Z ωωθθ⎡⎤⎢⎥=+=⎢⎥⎢⎥⎣⎦ ()11000100010v R v p ω=+=212212*********ˆs R Z c θωωθθθ⎡⎤⎢⎥=+=⎢⎥⎢⎥⎣⎦ ()22111211121100-v R v p L ωθ⎡⎤⎢⎥=+=⎢⎥⎢⎥⎣⎦23133233223323123ˆs R Z c θωωθθθθ⎡⎤⎢⎥=+=⎢⎥⎢⎥+⎣⎦ ()()23233222322232321221-L s v R v p L c L L c θωθθ⎡⎤⎢⎥=+=⎢⎥⎢⎥+⎣⎦ 231434323123s c θωωθθθ⎡⎤⎢⎥==⎢⎥⎢⎥+⎣⎦()()()23244333433342332331222231-L s v R v p L c L L L L c L c θωθθθ⎡⎤⎢⎥=+=++⎢⎥⎢⎥++⎣⎦ 为了得到这些速度相对于固定基坐标系的表达，用旋转矩阵04R 对他们作旋转变换，即：1231231001234123412312312323---0c c c s s R R R R R s c s s c s c ⎛⎫ ⎪== ⎪ ⎪⎝⎭通过这个变换可以得到：()()()()123123123200444412312312332332323122223111212312312123323312323123311212---0-- c c c s s L s v R v s c s s c L c L L s c L L c L c L s L s c L s c L c c s L c L c s L c s L c L c c L θθθθθθθ⎡⎤⎛⎫⎢⎥ ⎪==++⎢⎥⎪ ⎪⎢⎥++⎝⎭⎣⎦+++-+-⎡⎤⎣⎦=++ ()()()312312123323312323123312232332233323232330c c L s c s L c L s s L s s L s c s L c c L c L c θθθθθθ⎛⎫ ⎪+-+-⎡⎤ ⎪⎣⎦ ⎪ ⎪++++⎝⎭可以写出关节角速率传递到基坐标系下末端效应器速度的雅可比矩阵为：()()()11212312321233233123312301121231232123323312331232232332233323323-0L s L s c L s c L c c s L c L c s L c s J L c L c c L c c L s c s L c L s s L s s L s c s L c c L c L c θ++-+-⎛⎫⎪=++-+- ⎪ ⎪++⎝⎭2．应用拉格朗日公式推导该机器人的动力学模型：三自由度均质连杆的惯性张量为：121112110001*********c I m L m L ⎛⎫ ⎪⎪ ⎪= ⎪ ⎪ ⎪⎝⎭222222220001*********c I m L m L ⎛⎫ ⎪⎪ ⎪= ⎪ ⎪ ⎪⎝⎭323332330001*********c I m L m L ⎛⎫ ⎪⎪ ⎪= ⎪ ⎪ ⎪⎝⎭第i 个连杆的动能i k 可以表示为：1122i i ii iic TT i i i c c c c i ik m v v v T ωω=+则连杆1的动能为：222222111111111111111()222126k m L m L m L θθθ=+⋅⋅= 连杆2的动能为：222222221221222222222222122221221222111122221111 2626c T k m L L c L I m L m L c m L L c m L θθωωθθ⎡⎤⎛⎫⎛⎫=+++⋅⋅⎢⎥ ⎪ ⎪⎝⎭⎝⎭⎢⎥⎣⎦⎛⎫=+++ ⎪⎝⎭连杆3的动能为：33233333331122c T c k m v I ωω=+⋅⋅2223333cz cy cx c v v v v ++=33332232322232323223232223232311()2211()22c c c c x L c L c x L s L s y L s L s y l c L c θθθθθθ⎧⎧=+=-⋅-⋅+⎪⎪⎪⎪⇒⎨⎨⎪⎪=+=⋅+⋅+⎪⎪⎩⎩ 312232311()2c zL L c L c θ=++ 222222331332332331223132332322312222223233323323333332332311111()22622111111 ()()262632k m L m L c m l c m L L c m L L c m L L c c m L m L m L L c m L m l m l l c θθθθθ=+++++++++++ 总动能为：()222222132332333233232222222122323122313232322311114421 24k L L L L c L L L L c L L c L c L L c L L c L L c c θθθθθθθ⎛⎫⎛⎫=++++++ ⎪⎪⎝⎭⎝⎭⎛⎫+++++ ⎪⎝⎭，连杆1的势能：10u =连杆2的势能为：22222211sin 22u m gL m L gθ=+连杆3的势能为：332232332311()22u m g L s L s m g L L ⎛⎫=+++ ⎪⎝⎭总势能为：22323323223231111u()L s L gL L L 2222m m g m g s m m g θ⎛⎫⎛⎫=+++++ ⎪ ⎪⎝⎭⎝⎭对总动能和总势能分别求偏导：2222111212222122112222313223323312232323312223111331 + 23L m L m L m L c m L L c m L m L c m L c m L L c m L L c m L L c c θθθθ∂⎛⎫=+++ ⎪∂⎝⎭⎛⎫+++++ ⎪⎝⎭ 2222222323332332333233321111 3332L m L m L m L m L L c m L m L L c θθθθ∂⎛⎫⎛⎫=+++++ ⎪ ⎪∂⎝⎭⎝⎭2233333323323111332L m L m L m L L c θθθ∂⎛⎫=++ ⎪∂⎝⎭ 11122222222222121221222213323231222222231221313231323223132322310111323111 222L k uL k u m L s c m L L s m L s c L s c m L L s m L L s m L L s c m L L c c θθθθθθθθθθθθθθ∂∂∂=-=∂∂∂∂∂∂=-=----∂∂∂---- 2223223223222223323231213231323223132332333323323332311221111322211 22m gL c m gL c m gL c L k u m L s c m L L s L L c s m L L s m L L s m gL c θθθθθθθθθ---∂∂∂=-=----∂∂∂--1222222111211221121222121221212221311222222322221322133232323133231312133321 2()233L d m L m L m L s c m L c m L L c m L L s m L dt m L s c m L c m L c s m L c m θθθθθθθθθθθθθθθθθ⎛⎫∂ ⎪∂⎝⎭=+-++-+-+-++-1222131221313232313132313232322132322323132322312()() L L s m L L c m L L s m L L c m L L c s m L L c s m L L c c θθθθθθθθθθθθθ+-++--++ 22222222322332323323323323332323333233311133311 22L d m L m L m L m L L s m L L c m L dt m L L s m L L c θθθθθθθθθθ⎛⎫∂ ⎪∂⎝⎭=++-++-+3223333323233233233211113322L d m L m L m L L s L L c dt θθθθθθ⎛⎫∂ ⎪∂⎝⎭=+-+操作臂的运动方程为：d k k udt τθθθ∂∂∂-+=∂∂∂ 可求得：1222222111121122112122212122121222131112222322221322133232321213332 2(3L d L m L m L m L s c L c m L L c m L L s m L dt m L s c m L c m L c s θτθθθθθθθθθθθθθθ⎛⎫∂ ⎪∂∂⎝⎭=-=+-++-+∂-+-+ 22313323131222131221313232313132313232322132322323131)232() ()m L c m L L s m L L c m L L s m L L c m L L c s m L L c s m L θθθθθθθθθθθθθθθ+-+-++--++ 232231L c c θ 22222222232233232332333323332222222222222212122132221332323131221313231111133331111 3232 L d L m L m L m L m L L c L m L L c dt m L s c m L L s m L s c L s c m L L s m L L s θτθθθθθθθθθθθθθ⎛⎫∂ ⎪∂∂⎝⎭=-=+++++∂++++++ 2223232231323223132333323323222322332311111 22222m L L s c m L L c s L L s m L L s m L c g m L c g m L c g θθθθθ++--+++322222333232332333332323131323132232322313233233231111132332111 222L d L m L m L L c m L L s c m L L s dt m L L c s m L L s m L c g θτθθθθθθθθ⎛⎫∂ ⎪∂∂⎝⎭=-=++++∂+++ 由上式可以看出操作臂的动力学方程为：=()(,)()M V G τθθθθθ++ 其中：()2222112122221222222231322332331223132332322322222232333232333232223332323311331032111103332=1110323m L m L m L c m L L c m L m L c m L c m L L c m L L c m L L c c m L m L m L m L L c m L m L L c M m L m L L c m L θ⎛+++++++++ ++++++⎝⎫⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪ ⎪⎭()222221121212221322221332323231312221313232313232322132322323122222222121221322222()233()()1132m L s c m L L s m L s c m L c s m L L s m L L s m L L c s m L L c s m L s c m L L s m L s V θθθθθθθθθθθθθθθθθθθθθθθ----+--+--+++= ，22222213323231312213132312223232231323223132333323323222223323231313231323223132332113211122211113222c m L s c m L L s m L L s m L L s c m L L c s m L L s m L L s m L s c m L L s m L L c s m L L s θθθθθθθθθθθθθ⎡⎢⎢⎢⎢⎢⎢+++⎢⎢⎢++--⎢⎢⎢⎢+++⎣ ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎦()222322332333230112212G m L c g m L c g m L c g m L c g θ⎡⎤⎢⎥⎢⎥⎢⎥=++⎢⎥⎢⎥⎢⎥⎣⎦3、机器人轨迹跟踪控制器设置由于该空间三自由度机器人模型为非线性系统，并且需要跟踪轨迹，是一个多输入多输出控制系统，故将控制规律分解成为基于模型的控制部分和伺服控制部分，利用α-β分解法，将运动方程解耦，得到线性化解耦控制规律，选择适合的K p 和K v ，利用PD 控制完成轨迹跟踪。

黄永刚单晶塑性有限元umat子程序

SUBROUTINE UMAT(stress,statev,ddsdde,sse,spd,scd,1 rpl, ddsddt, drplde, drpldt,2 stran,dstran,time,dtime,temp,dtemp,predef,dpred,cmname,3 ndi,nshr,ntens,nstatv,props,nprops,coords,drot,pnewdt,4 celent,dfgrd0,dfgrd1,noel,npt,layer,kspt,kstep,kinc)c WRITE (6,*) 'c NOTE: MODIFICATIONS TO *UMAT FOR ABAQUS VERSION 5.3 (14 APR '94) cc (1) The list of variables above defining the *UMAT subroutine,c and the first (standard) block of variables dimensioned below,c have variable names added compared to earlier ABAQUS versions.cc (2) The statement: include 'aba_param.inc' must be added as below.cc (3) As of version 5.3, ABAQUS files use double precision only.c The file aba_param.inc has a line "implicit real*8" and, sincec it is included in the main subroutine, it will define the variablesc there as double precision. But other subroutines still need thec definition "implicit real*8" since there may be variables that arec not passed to them through the list or common block.cc (4) This is current as of version 5.6 of ABAQUS.cc (5) Note added by J. W. Kysar (4 November 1997). This UMAT has beenc modified to keep track of the cumulative shear strain in eachc individual slip system. This information is needed to correct anc error in the implementation of the Bassani and Wu hardening law.c Any line of code which has been added or modified is precededc immediately by a line beginning CFIXA and succeeded by a linec beginning CFIXB. Any comment line added or modified will beginc with CFIX.cc The hardening law by Bassani and Wu was implemented incorrectly.c This law is a function of both hyperbolic secant squared and hyperbolicc tangent. However, the arguments of sech and tanh are related to the *total*c slip on individual slip systems. Formerly, the UMAT implemented thisc hardening law by using the *current* slip on each slip system. Thereinc lay the problem. The UMAT did not restrict the current slip to be ac positive value. So when a slip with a negative sign was encountered, thec term containing tanh led to a negative hardening rate (since tanh is anc odd function).c The UMA T has been fixed by adding state variables to keep track of thec *total* slip on each slip system by integrating up the absolute valuec of slip rates for each individual slip system. These "solution dependentc variables" are available for postprocessing. The only required changec in the input file is that the DEPV AR command must be changed.cC----- Use single precision on Cray byC (1) deleting the statement "IMPLICIT*8 (A-H,O-Z)";C (2) changing "REAL*8 FUNCTION" to "FUNCTION";C (3) changing double precision functions DSIGN to SIGN.CC----- Subroutines:CC ROTATION -- forming rotation matrix, i.e. the directionC cosines of cubic crystal [100], [010] and [001]C directions in global system at the initialC stateCC SLIPSYS -- calculating number of slip systems, unitC vectors in slip directions and unit normals toC slip planes in a cubic crystal at the initialC stateCC GSLPINIT -- calculating initial value of current strengthsC at initial stateCC STRAINRA TE -- based on current values of resolved shearC stresses and current strength, calculatingC shear strain-rates in slip systemsCC LATENTHARDEN -- forming self- and latent-hardening matrixCC ITERATION -- generating arrays for the Newton-RhapsonC iterationCC LUDCMP -- LU decompositionCC LUBKSB -- linear equation solver based on LUC decomposition method (must call LUDCMP first) C----- Function subprogram:C F -- shear strain-rates in slip systemsC----- Variables:CC STRESS -- stresses (INPUT & OUTPUT)C Cauchy stresses for finite deformationC STATEV -- solution dependent state variables (INPUT & OUTPUT) C DDSDDE -- Jacobian matrix (OUTPUT)C----- Variables passed in for information:CC STRAN -- strainsC logarithmic strain for finite deformationC (actually, integral of the symmetric part of velocityC gradient with respect to time)C DSTRAN -- increments of strainsC CMNAME -- name given in the *MATERIAL optionC NDI -- number of direct stress componentsC NSHR -- number of engineering shear stress componentsC NTENS -- NDI+NSHRC NSTATV -- number of solution dependent state variables (asC defined in the *DEPVAR option)C PROPS -- material constants entered in the *USER MA TERIAL C optionC NPROPS -- number of material constantsCC----- This subroutine provides the plastic constitutive relation ofC single crystals for finite element code ABAQUS. The plastic slipC of single crystal obeys the Schmid law. The program gives theC choice of small deformation theory and theory of finite rotationC and finite strain.C The strain increment is composed of elastic part and plasticC part. The elastic strain increment corresponds to latticeC stretching, the plastic part is the sum over all slip systems ofC plastic slip. The shear strain increment for each slip system isC assumed a function of the ratio of corresponding resolved shearC stress over current strength, and of the time step. The resolvedC shear stress is the double product of stress tensor with the slipC deformation tensor (Schmid factor), and the increment of currentC strength is related to shear strain increments over all slipC systems through self- and latent-hardening functions.C----- The implicit integration method proposed by Peirce, Shih andC Needleman (1984) is used here. The subroutine provides an option C of iteration to solve stresses and solution dependent stateC variables within each increment.C----- The present program is for a single CUBIC crystal. However,C this code can be generalized for other crystals (e.g. HCP,C Tetragonal, Orthotropic, etc.). Only subroutines ROTATION andC SLIPSYS need to be modified to include the effect of crystalC aspect ratio.CC----- Important notice:CC (1) The number of state variables NSTATV must be larger than (or CFIX equal to) TEN (10) times the total number of slip systems inC all sets, NSLPTL, plus FIVE (5)CFIX NSTATV >= 10 * NSLPTL + 5C Denote s as a slip direction and m as normal to a slip plane.C Here (s,-m), (-s,m) and (-s,-m) are NOT consideredC independent of (s,m). The number of slip systems in each setC could be either 6, 12, 24 or 48 for a cubic crystal, e.g. 12C for {110}<111>.CC Users who need more parameters to characterize theC constitutive law of single crystal, e.g. the frameworkC proposed by Zarka, should make NSTATV larger than (or equal C to) the number of those parameters NPARMT plus nine timesC the total number of slip systems, NSLPTL, plus fiveCFIX NSTATV >= NPARMT + 10 * NSLPTL + 5CC (2) The tangent stiffness matrix in general is not symmetric ifC latent hardening is considered. Users must declare "UNSYMM"C in the input file, at the *USER MATERIAL card.CPARAMETER (ND=150)C----- The parameter ND determines the dimensions of the arrays inC this subroutine. The current choice 150 is a upper bound for aC cubic crystal with up to three sets of slip systems activated.C Users may reduce the parameter ND to any number as long as larger C than or equal to the total number of slip systems in all sets.C For example, if {110}<111> is the only set of slip systemC potentially activated, ND could be taken as twelve (12).cinclude 'aba_param.inc'cCHARACTER*8 CMNAMEEXTERNAL Fdimension stress(ntens),statev(nstatv),1 ddsdde(ntens,ntens),ddsddt(ntens),drplde(ntens),2 stran(ntens),dstran(ntens),time(2),predef(1),dpred(1),3 props(nprops),coords(3),drot(3,3),dfgrd0(3,3),dfgrd1(3,3)DIMENSION ISPDIR(3), ISPNOR(3), NSLIP(3),2 SLPDIR(3,ND), SLPNOR(3,ND), SLPDEF(6,ND),3 SLPSPN(3,ND), DSPDIR(3,ND), DSPNOR(3,ND),4 DLOCAL(6,6), D(6,6), ROTD(6,6), ROTATE(3,3),5 FSLIP(ND), DFDXSP(ND), DDEMSD(6,ND),6 H(ND,ND), DDGDDE(ND,6),7 DSTRES(6), DELATS(6), DSPIN(3), DVGRAD(3,3),8 DGAMMA(ND), DTAUSP(ND), DGSLIP(ND),9 WORKST(ND,ND), INDX(ND), TERM(3,3), TRM0(3,3), ITRM(3)DIMENSION FSLIP1(ND), STRES1(6), GAMMA1(ND), TAUSP1(ND),2 GSLP1(ND), SPNOR1(3,ND), SPDIR1(3,ND), DDSDE1(6,6),3 DSOLD(6), DGAMOD(ND), DTAUOD(ND), DGSPOD(ND),4 DSPNRO(3,ND), DSPDRO(3,ND),5 DHDGDG(ND,ND)C----- NSLIP -- number of slip systems in each setC----- SLPDIR -- slip directions (unit vectors in the initial state)C----- SLPNOR -- normals to slip planes (unit normals in the initialC state)C----- SLPDEF -- slip deformation tensors (Schmid factors)C SLPDEF(1,i) -- SLPDIR(1,i)*SLPNOR(1,i)C SLPDEF(2,i) -- SLPDIR(2,i)*SLPNOR(2,i)C SLPDEF(3,i) -- SLPDIR(3,i)*SLPNOR(3,i)C SLPDEF(4,i) -- SLPDIR(1,i)*SLPNOR(2,i)+C SLPDIR(2,i)*SLPNOR(1,i)C SLPDEF(5,i) -- SLPDIR(1,i)*SLPNOR(3,i)+C SLPDIR(3,i)*SLPNOR(1,i)C SLPDEF(6,i) -- SLPDIR(2,i)*SLPNOR(3,i)+C SLPDIR(3,i)*SLPNOR(2,i)C where index i corresponds to the ith slip systemC----- SLPSPN -- slip spin tensors (only needed for finite rotation)C SLPSPN(1,i) -- [SLPDIR(1,i)*SLPNOR(2,i)-C SLPDIR(2,i)*SLPNOR(1,i)]/2C SLPSPN(2,i) -- [SLPDIR(3,i)*SLPNOR(1,i)-C SLPDIR(1,i)*SLPNOR(3,i)]/2C SLPSPN(3,i) -- [SLPDIR(2,i)*SLPNOR(3,i)-C SLPDIR(3,i)*SLPNOR(2,i)]/2C where index i corresponds to the ith slip systemC----- DSPDIR -- increments of slip directionsC----- DSPNOR -- increments of normals to slip planesCC----- DLOCAL -- elastic matrix in local cubic crystal systemC----- D -- elastic matrix in global systemC----- ROTD -- rotation matrix transforming DLOCAL to DCC----- ROTATE -- rotation matrix, direction cosines of [100], [010]C and [001] of cubic crystal in global systemCC----- FSLIP -- shear strain-rates in slip systemsC----- DFDXSP -- derivatives of FSLIP w.r.t x=TAUSLP/GSLIP, whereC TAUSLP is the resolved shear stress and GSLIP is the C current strengthCC----- DDEMSD -- double dot product of the elastic moduli tensor withC the slip deformation tensor plus, only for finiteC rotation, the dot product of slip spin tensor withC the stressCC----- H -- self- and latent-hardening matrixC H(i,i) -- self hardening modulus of the ith slipC system (no sum over i)C H(i,j) -- latent hardening molulus of the ith slipC system due to a slip in the jth slip system C (i not equal j)CC----- DDGDDE -- derivatice of the shear strain increments in slipC systems w.r.t. the increment of strainsCC----- DSTRES -- Jaumann increments of stresses, i.e. corotationalC stress-increments formed on axes spinning with theC materialC----- DELATS -- strain-increments associated with lattice stretchingC DELATS(1) - DELATS(3) -- normal strain increments C DELATS(4) - DELATS(6) -- engineering shear strain C incrementsC----- DSPIN -- spin-increments associated with the material elementC DSPIN(1) -- component 12 of the spin tensorC DSPIN(2) -- component 31 of the spin tensorC DSPIN(3) -- component 23 of the spin tensorCC----- DVGRAD -- increments of deformation gradient in the currentC state, i.e. velocity gradient times the increment ofC timeCC----- DGAMMA -- increment of shear strains in slip systemsC----- DTAUSP -- increment of resolved shear stresses in slip systemsC----- DGSLIP -- increment of current strengths in slip systemsCCC----- Arrays for iteration:CC FSLIP1, STRES1, GAMMA1, TAUSP1, GSLP1 , SPNOR1, SPDIR1,C DDSDE1, DSOLD , DGAMOD, DTAUOD, DGSPOD, DSPNRO, DSPDRO, C DHDGDGCCC----- Solution dependent state variable STATEV:C Denote the number of total slip systems by NSLPTL, whichC will be calculated in this code.CC Array STATEV:C 1 - NSLPTL : current strength in slip systemsC NSLPTL+1 - 2*NSLPTL : shear strain in slip systemsC 2*NSLPTL+1 - 3*NSLPTL : resolved shear stress in slip systemsCC 3*NSLPTL+1 - 6*NSLPTL : current components of normals to slipC slip planesC 6*NSLPTL+1 - 9*NSLPTL : current components of slip directionsCCFIX 9*NSLPTL+1 - 10*NSLPTL : total cumulative shear strain on eachCFIX slip system (sum of the absoluteCFIX values of shear strains in each slipCFIX system individually)CFIXCFIX 10*NSLPTL+1 : total cumulative shear strain on allC slip systems (sum of the absoluteC values of shear strains in all slipC systems)CCFIX 10*NSLPTL+2 - NSTA TV-4 : additional parameters users may needC to characterize the constitutive lawC of a single crystal (if there areC any).CC NSTATV-3 : number of slip systems in the 1st setC NSTATV-2 : number of slip systems in the 2nd setC NSTATV-1 : number of slip systems in the 3rd setC NSTATV : total number of slip systems in allC setsCCC----- Material constants PROPS:CC PROPS(1) - PROPS(21) -- elastic constants for a general elasticC anisotropic materialCC isotropic : PROPS(i)=0 for i>2C PROPS(1) -- Young's modulusC PROPS(2) -- Poisson's ratioCC cubic : PROPS(i)=0 for i>3C PROPS(1) -- c11C PROPS(2) -- c12C PROPS(3) -- c44CC orthotropic : PORPS(i)=0 for i>9C PROPS(1) - PROPS(9) are D1111, D1122, D2222, C D1133, D2233, D3333, D1212, D1313, D2323,C respectively, which has the same definitionC as ABAQUS for orthotropic materialsC (see *ELASTIC card)CC anisotropic : PROPS(1) - PROPS(21) are D1111, D1122,C D2222, D1133, D2233, D3333, D1112, D2212,C D3312, D1212, D1113, D2213, D3313, D1213,C D1313, D1123, D2223, D3323, D1223, D1323,C D2323, respectively, which has the sameC definition as ABAQUS for anisotropicC materials (see *ELASTIC card)CCC PROPS(25) - PROPS(56) -- parameters characterizing all slipC systems to be activated in a cubicC crystalCC PROPS(25) -- number of sets of slip systems (maximum 3),C e.g. (110)[1-11] and (101)[11-1] are in theC same set of slip systems, (110)[1-11] andC (121)[1-11] belong to different sets of slipC systemsC (It must be a real number, e.g. 3., not 3 !)CC PROPS(33) - PROPS(35) -- normal to a typical slip plane inC the first set of slip systems,C e.g. (1 1 0)C (They must be real numbers, e.g.C 1. 1. 0., not 1 1 0 !)C PROPS(36) - PROPS(38) -- a typical slip direction in theC first set of slip systems, e.g.C [1 1 1]C (They must be real numbers, e.g.C 1. 1. 1., not 1 1 1 !)CC PROPS(41) - PROPS(43) -- normal to a typical slip plane inC the second set of slip systemsC (real numbers)C PROPS(44) - PROPS(46) -- a typical slip direction in theC second set of slip systemsC (real numbers)CC PROPS(49) - PROPS(51) -- normal to a typical slip plane inC the third set of slip systemsC (real numbers)C PROPS(52) - PROPS(54) -- a typical slip direction in theC third set of slip systemsC (real numbers)CCC PROPS(57) - PROPS(72) -- parameters characterizing the initialC orientation of a single crystal inC global systemC The directions in global system and directions in localC cubic crystal system of two nonparallel vectors are neededC to determine the crystal orientation.CC PROPS(57) - PROPS(59) -- [p1 p2 p3], direction of firstC vector in local cubic crystalC system, e.g. [1 1 0]C (They must be real numbers, e.g.C 1. 1. 0., not 1 1 0 !)C PROPS(60) - PROPS(62) -- [P1 P2 P3], direction of firstC vector in global system, e.g.C [2. 1. 0.]C (It does not have to be a unitC vector)CC PROPS(65) - PROPS(67) -- direction of second vector inC local cubic crystal system (real C numbers)C PROPS(68) - PROPS(70) -- direction of second vector inC global systemCCC PROPS(73) - PROPS(96) -- parameters characterizing the visco-C plastic constitutive law (shearC strain-rate vs. resolved shearC stress), e.g. a power-law relationCC PROPS(73) - PROPS(80) -- parameters for the first set ofC slip systemsC PROPS(81) - PROPS(88) -- parameters for the second set of C slip systemsC PROPS(89) - PROPS(96) -- parameters for the third set ofC slip systemsCCC PROPS(97) - PROPS(144)-- parameters characterizing the self-C and latent-hardening laws of slipC systemsCC PROPS(97) - PROPS(104)-- self-hardening parameters for the C first set of slip systemsC PROPS(105)- PROPS(112)-- latent-hardening parameters for C the first set of slip systems and C interaction with other sets ofC slip systemsCC PROPS(113)- PROPS(120)-- self-hardening parameters for the C second set of slip systemsC PROPS(121)- PROPS(128)-- latent-hardening parameters for C the second set of slip systems C and interaction with other sets C of slip systemsCC PROPS(129)- PROPS(136)-- self-hardening parameters for the C third set of slip systemsC PROPS(137)- PROPS(144)-- latent-hardening parameters for C the third set of slip systems and C interaction with other sets ofC slip systemsCCC PROPS(145)- PROPS(152)-- parameters characterizing forward time C integration scheme and finiteC deformationCC PROPS(145) -- parameter theta controlling the implicitC integration, which is between 0 and 1C 0. : explicit integrationC 0.5 : recommended valueC 1. : fully implicit integrationCC PROPS(146) -- parameter NLGEOM controlling whether the C effect of finite rotation and finite strainC of crystal is considered,C 0. : small deformation theoryC otherwise : theory of finite rotation andC finite strainCCC PROPS(153)- PROPS(160)-- parameters characterizing iterationC methodCC PROPS(153) -- parameter ITRATN controlling whether theC iteration method is used,C 0. : no iterationC otherwise : iterationCC PROPS(154) -- maximum number of iteration ITRMAXCC PROPS(155) -- absolute error of shear strains in slipC systems GAMERRCC----- Elastic matrix in local cubic crystal system: DLOCALDO J=1,6DO I=1,6DLOCAL(I,J)=0.END DOEND DOCHECK=0.DO J=10,21CHECK=CHECK+ABS(PROPS(J))END DOIF (CHECK.EQ.0.) THENDO J=4,9CHECK=CHECK+ABS(PROPS(J))END DOIF (CHECK.EQ.0.) THENIF (PROPS(3).EQ.0.) THENC----- Isotropic materialGSHEAR=PROPS(1)/2./(1.+PROPS(2))E11=2.*GSHEAR*(1.-PROPS(2))/(1.-2.*PROPS(2))E12=2.*GSHEAR*PROPS(2)/(1.-2.*PROPS(2))DO J=1,3DLOCAL(J,J)=E11DO I=1,3IF (I.NE.J) DLOCAL(I,J)=E12END DODLOCAL(J+3,J+3)=GSHEAREND DOELSEC----- Cubic materialDO J=1,3DLOCAL(J,J)=PROPS(1)DO I=1,3IF (I.NE.J) DLOCAL(I,J)=PROPS(2)END DODLOCAL(J+3,J+3)=PROPS(3)END DOEND IFELSEC----- Orthotropic metarialDLOCAL(1,1)=PROPS(1)DLOCAL(1,2)=PROPS(2)DLOCAL(2,1)=PROPS(2)DLOCAL(2,2)=PROPS(3)DLOCAL(1,3)=PROPS(4)DLOCAL(3,1)=PROPS(4)DLOCAL(2,3)=PROPS(5)DLOCAL(3,2)=PROPS(5)DLOCAL(3,3)=PROPS(6)DLOCAL(4,4)=PROPS(7)DLOCAL(5,5)=PROPS(8)DLOCAL(6,6)=PROPS(9)END IFELSEC----- General anisotropic materialID=0DO J=1,6DO I=1,JID=ID+1DLOCAL(I,J)=PROPS(ID)DLOCAL(J,I)=DLOCAL(I,J)END DOEND DOEND IFC----- Rotation matrix: ROTA TE, i.e. direction cosines of [100], [010]C and [001] of a cubic crystal in global systemCCALL ROTATION (PROPS(57), ROTA TE)C----- Rotation matrix: ROTD to transform local elastic matrix DLOCAL C to global elastic matrix DCDO J=1,3J1=1+J/3J2=2+J/2DO I=1,3I1=1+I/3I2=2+I/2ROTD(I,J)=ROTATE(I,J)**2ROTD(I,J+3)=2.*ROTATE(I,J1)*ROTA TE(I,J2)ROTD(I+3,J)=ROTA TE(I1,J)*ROTATE(I2,J)ROTD(I+3,J+3)=ROTA TE(I1,J1)*ROTA TE(I2,J2)+2 ROTA TE(I1,J2)*ROTA TE(I2,J1)END DOEND DOC----- Elastic matrix in global system: DC {D} = {ROTD} * {DLOCAL} * {ROTD}transposeCDO J=1,6DO I=1,6D(I,J)=0.END DOEND DODO J=1,6DO I=1,JDO K=1,6DO L=1,6D(I,J)=D(I,J)+DLOCAL(K,L)*ROTD(I,K)*ROTD(J,L) END DOEND DOD(J,I)=D(I,J)END DOEND DOC----- Total number of sets of slip systems: NSETNSET=NINT(PROPS(25))IF (NSET.LT.1) THENWRITE (6,*) '***ERROR - zero sets of slip systems'STOPELSE IF (NSET.GT.3) THENWRITE (6,*)2 '***ERROR - more than three sets of slip systems'STOPEND IFC----- Implicit integration parameter: THETATHETA=PROPS(145)C----- Finite deformation ?C----- NLGEOM = 0, small deformation theoryC otherwise, theory of finite rotation and finite strain, Users C must declare "NLGEOM" in the input file, at the *STEP card CIF (PROPS(146).EQ.0.) THENNLGEOM=0ELSENLGEOM=1END IFC----- Iteration?C----- ITRATN = 0, no iterationC otherwise, iteration (solving increments of stresses andC solution dependent state variables)CIF (PROPS(153).EQ.0.) THENITRATN=0ELSEITRATN=1END IFITRMAX=NINT(PROPS(154))GAMERR=PROPS(155)NITRTN=-1DO I=1,NTENSDSOLD(I)=0.END DODO J=1,NDDGAMOD(J)=0.DTAUOD(J)=0.DGSPOD(J)=0.DO I=1,3DSPNRO(I,J)=0.DSPDRO(I,J)=0.END DOEND DOC----- Increment of spin associated with the material element: DSPIN C (only needed for finite rotation)CIF (NLGEOM.NE.0) THENDO J=1,3DO I=1,3TERM(I,J)=DROT(J,I)TRM0(I,J)=DROT(J,I)END DOTERM(J,J)=TERM(J,J)+1.D0TRM0(J,J)=TRM0(J,J)-1.D0END DOCALL LUDCMP (TERM, 3, 3, ITRM, DDCMP)DO J=1,3CALL LUBKSB (TERM, 3, 3, ITRM, TRM0(1,J)) END DODSPIN(1)=TRM0(2,1)-TRM0(1,2)DSPIN(2)=TRM0(1,3)-TRM0(3,1)DSPIN(3)=TRM0(3,2)-TRM0(2,3)END IFC----- Increment of dilatational strain: DEVDEV=0.D0DO I=1,NDIDEV=DEV+DSTRAN(I)END DOC----- Iteration starts (only when iteration method is used)1000 CONTINUEC----- Parameter NITRTN: number of iterationsC NITRTN = 0 --- no-iteration solutionCNITRTN=NITRTN+1C----- Check whether the current stress state is the initial stateIF (STATEV(1).EQ.0.) THENC----- Initial stateCC----- Generating the following parameters and variables at initialC state:C Total number of slip systems in all the sets NSLPTLC Number of slip systems in each set NSLIPC Unit vectors in initial slip directions SLPDIRC Unit normals to initial slip planes SLPNORCNSLPTL=0DO I=1,NSETISPNOR(1)=NINT(PROPS(25+8*I))ISPNOR(2)=NINT(PROPS(26+8*I))ISPNOR(3)=NINT(PROPS(27+8*I))ISPDIR(1)=NINT(PROPS(28+8*I))ISPDIR(2)=NINT(PROPS(29+8*I))ISPDIR(3)=NINT(PROPS(30+8*I))CALL SLIPSYS (ISPDIR, ISPNOR, NSLIP(I), SLPDIR(1,NSLPTL+1),2 SLPNOR(1,NSLPTL+1), ROTATE)NSLPTL=NSLPTL+NSLIP(I)END DOIF (ND.LT.NSLPTL) THENWRITE (6,*)2 '***ERROR - parameter ND chosen by the present user is less than3 the total number of slip systems NSLPTL'STOPEND IFC----- Slip deformation tensor: SLPDEF (Schmid factors)DO J=1,NSLPTLSLPDEF(1,J)=SLPDIR(1,J)*SLPNOR(1,J)SLPDEF(2,J)=SLPDIR(2,J)*SLPNOR(2,J)SLPDEF(3,J)=SLPDIR(3,J)*SLPNOR(3,J)SLPDEF(4,J)=SLPDIR(1,J)*SLPNOR(2,J)+SLPDIR(2,J)*SLPNOR(1,J)SLPDEF(5,J)=SLPDIR(1,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(1,J)SLPDEF(6,J)=SLPDIR(2,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(2,J) END DOC----- Initial value of state variables: unit normal to a slip planeC and unit vector in a slip directionCSTATEV(NSTATV)=FLOAT(NSLPTL)DO I=1,NSETSTA TEV(NSTA TV-4+I)=FLOAT(NSLIP(I))END DOIDNOR=3*NSLPTLIDDIR=6*NSLPTLDO J=1,NSLPTLDO I=1,3IDNOR=IDNOR+1STA TEV(IDNOR)=SLPNOR(I,J)IDDIR=IDDIR+1STA TEV(IDDIR)=SLPDIR(I,J)END DOEND DOC----- Initial value of the current strength for all slip systemsCCALL GSLPINIT (STATEV(1), NSLIP, NSLPTL, NSET, PROPS(97))C----- Initial value of shear strain in slip systemsCFIX-- Initial value of cumulative shear strain in each slip systemsDO I=1,NSLPTLSTA TEV(NSLPTL+I)=0.CFIXASTA TEV(9*NSLPTL+I)=0.CFIXBEND DOCFIXASTATEV(10*NSLPTL+1)=0.CFIXBC----- Initial value of the resolved shear stress in slip systemsDO I=1,NSLPTLTERM1=0.DO J=1,NTENSIF (J.LE.NDI) THENTERM1=TERM1+SLPDEF(J,I)*STRESS(J)ELSETERM1=TERM1+SLPDEF(J-NDI+3,I)*STRESS(J) END IFEND DOSTA TEV(2*NSLPTL+I)=TERM1END DOELSEC----- Current stress stateCC----- Copying from the array of state variables STA TVE the following C parameters and variables at current stress state:C Total number of slip systems in all the sets NSLPTLC Number of slip systems in each set NSLIPC Current slip directions SLPDIRC Normals to current slip planes SLPNORCNSLPTL=NINT(STA TEV(NSTATV))DO I=1,NSETNSLIP(I)=NINT(STATEV(NSTATV-4+I))END DOIDNOR=3*NSLPTLIDDIR=6*NSLPTLDO J=1,NSLPTLDO I=1,3IDNOR=IDNOR+1SLPNOR(I,J)=STATEV(IDNOR)IDDIR=IDDIR+1SLPDIR(I,J)=STA TEV(IDDIR)END DOEND DOC----- Slip deformation tensor: SLPDEF (Schmid factors)DO J=1,NSLPTLSLPDEF(1,J)=SLPDIR(1,J)*SLPNOR(1,J)SLPDEF(2,J)=SLPDIR(2,J)*SLPNOR(2,J)SLPDEF(3,J)=SLPDIR(3,J)*SLPNOR(3,J)SLPDEF(4,J)=SLPDIR(1,J)*SLPNOR(2,J)+SLPDIR(2,J)*SLPNOR(1,J)SLPDEF(5,J)=SLPDIR(1,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(1,J)SLPDEF(6,J)=SLPDIR(2,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(2,J) END DOEND IFC----- Slip spin tensor: SLPSPN (only needed for finite rotation)IF (NLGEOM.NE.0) THENDO J=1,NSLPTLSLPSPN(1,J)=0.5*(SLPDIR(1,J)*SLPNOR(2,J)-2 SLPDIR(2,J)*SLPNOR(1,J))SLPSPN(2,J)=0.5*(SLPDIR(3,J)*SLPNOR(1,J)-2 SLPDIR(1,J)*SLPNOR(3,J))SLPSPN(3,J)=0.5*(SLPDIR(2,J)*SLPNOR(3,J)-2 SLPDIR(3,J)*SLPNOR(2,J))END DOEND IFC----- Double dot product of elastic moduli tensor with the slipC deformation tensor (Schmid factors) plus, only for finiteC rotation, the dot product of slip spin tensor with the stress:C DDEMSDCDO J=1,NSLPTLDO I=1,6DDEMSD(I,J)=0.DO K=1,6DDEMSD(I,J)=DDEMSD(I,J)+D(K,I)*SLPDEF(K,J)END DOEND DOEND DOIF (NLGEOM.NE.0) THENDO J=1,NSLPTLDDEMSD(4,J)=DDEMSD(4,J)-SLPSPN(1,J)*STRESS(1)DDEMSD(5,J)=DDEMSD(5,J)+SLPSPN(2,J)*STRESS(1)IF (NDI.GT.1) THENDDEMSD(4,J)=DDEMSD(4,J)+SLPSPN(1,J)*STRESS(2)DDEMSD(6,J)=DDEMSD(6,J)-SLPSPN(3,J)*STRESS(2)END IFIF (NDI.GT.2) THENDDEMSD(5,J)=DDEMSD(5,J)-SLPSPN(2,J)*STRESS(3)DDEMSD(6,J)=DDEMSD(6,J)+SLPSPN(3,J)*STRESS(3)END IFIF (NSHR.GE.1) THENDDEMSD(1,J)=DDEMSD(1,J)+SLPSPN(1,J)*STRESS(NDI+1)DDEMSD(2,J)=DDEMSD(2,J)-SLPSPN(1,J)*STRESS(NDI+1)DDEMSD(5,J)=DDEMSD(5,J)-SLPSPN(3,J)*STRESS(NDI+1)DDEMSD(6,J)=DDEMSD(6,J)+SLPSPN(2,J)*STRESS(NDI+1) END IFIF (NSHR.GE.2) THENDDEMSD(1,J)=DDEMSD(1,J)-SLPSPN(2,J)*STRESS(NDI+2)DDEMSD(3,J)=DDEMSD(3,J)+SLPSPN(2,J)*STRESS(NDI+2)DDEMSD(4,J)=DDEMSD(4,J)+SLPSPN(3,J)*STRESS(NDI+2)DDEMSD(6,J)=DDEMSD(6,J)-SLPSPN(1,J)*STRESS(NDI+2) END IFIF (NSHR.EQ.3) THENDDEMSD(2,J)=DDEMSD(2,J)+SLPSPN(3,J)*STRESS(NDI+3)DDEMSD(3,J)=DDEMSD(3,J)-SLPSPN(3,J)*STRESS(NDI+3)DDEMSD(4,J)=DDEMSD(4,J)-SLPSPN(2,J)*STRESS(NDI+3)DDEMSD(5,J)=DDEMSD(5,J)+SLPSPN(1,J)*STRESS(NDI+3) END IFEND DOEND IFC----- Shear strain-rate in a slip system at the start of increment:C FSLIP, and its derivative: DFDXSPCID=1DO I=1,NSETIF (I.GT.1) ID=ID+NSLIP(I-1)CALL STRAINRATE (STATEV(NSLPTL+ID), STATEV(2*NSLPTL+ID),2 STA TEV(ID), NSLIP(I), FSLIP(ID), DFDXSP(ID),3 PROPS(65+8*I))END DOC----- Self- and latent-hardening laws。

黄永刚单晶塑性有限元umat子程序

SUBROUTINE UMAT(stress,statev,ddsdde,sse,spd,scd,1 rpl, ddsddt, drplde, drpldt,2 stran,dstran,time,dtime,temp,dtemp,predef,dpred,cmname,3 ndi,nshr,ntens,nstatv,props,nprops,coords,drot,pnewdt,4 celent,dfgrd0,dfgrd1,noel,npt,layer,kspt,kstep,kinc)c WRITE (6,*) 'c NOTE: MODIFICATIONS TO *UMAT FOR ABAQUS VERSION 5.3 (14 APR '94) cc (1) The list of variables above defining the *UMAT subroutine,c and the first (standard) block of variables dimensioned below,c have variable names added compared to earlier ABAQUS versions.cc (2) The statement: include 'aba_param.inc' must be added as below.cc (3) As of version 5.3, ABAQUS files use double precision only.c The file aba_param.inc has a line "implicit real*8" and, sincec it is included in the main subroutine, it will define the variablesc there as double precision. But other subroutines still need thec definition "implicit real*8" since there may be variables that arec not passed to them through the list or common block.cc (4) This is current as of version 5.6 of ABAQUS.cc (5) Note added by J. W. Kysar (4 November 1997). This UMAT has beenc modified to keep track of the cumulative shear strain in eachc individual slip system. This information is needed to correct anc error in the implementation of the Bassani and Wu hardening law.c Any line of code which has been added or modified is precededc immediately by a line beginning CFIXA and succeeded by a linec beginning CFIXB. Any comment line added or modified will beginc with CFIX.cc The hardening law by Bassani and Wu was implemented incorrectly.c This law is a function of both hyperbolic secant squared and hyperbolicc tangent. However, the arguments of sech and tanh are related to the *total*c slip on individual slip systems. Formerly, the UMAT implemented thisc hardening law by using the *current* slip on each slip system. Thereinc lay the problem. The UMAT did not restrict the current slip to be ac positive value. So when a slip with a negative sign was encountered, thec term containing tanh led to a negative hardening rate (since tanh is anc odd function).c The UMA T has been fixed by adding state variables to keep track of thec *total* slip on each slip system by integrating up the absolute valuec of slip rates for each individual slip system. These "solution dependentc variables" are available for postprocessing. The only required changec in the input file is that the DEPV AR command must be changed.cC----- Use single precision on Cray byC (1) deleting the statement "IMPLICIT*8 (A-H,O-Z)";C (2) changing "REAL*8 FUNCTION" to "FUNCTION";C (3) changing double precision functions DSIGN to SIGN.CC----- Subroutines:CC ROTATION -- forming rotation matrix, i.e. the directionC cosines of cubic crystal [100], [010] and [001]C directions in global system at the initialC stateCC SLIPSYS -- calculating number of slip systems, unitC vectors in slip directions and unit normals toC slip planes in a cubic crystal at the initialC stateCC GSLPINIT -- calculating initial value of current strengthsC at initial stateCC STRAINRA TE -- based on current values of resolved shearC stresses and current strength, calculatingC shear strain-rates in slip systemsCC LATENTHARDEN -- forming self- and latent-hardening matrixCC ITERATION -- generating arrays for the Newton-RhapsonC iterationCC LUDCMP -- LU decompositionCC LUBKSB -- linear equation solver based on LUC decomposition method (must call LUDCMP first) C----- Function subprogram:C F -- shear strain-rates in slip systemsC----- Variables:CC STRESS -- stresses (INPUT & OUTPUT)C Cauchy stresses for finite deformationC STATEV -- solution dependent state variables (INPUT & OUTPUT) C DDSDDE -- Jacobian matrix (OUTPUT)C----- Variables passed in for information:CC STRAN -- strainsC logarithmic strain for finite deformationC (actually, integral of the symmetric part of velocityC gradient with respect to time)C DSTRAN -- increments of strainsC CMNAME -- name given in the *MATERIAL optionC NDI -- number of direct stress componentsC NSHR -- number of engineering shear stress componentsC NTENS -- NDI+NSHRC NSTATV -- number of solution dependent state variables (asC defined in the *DEPVAR option)C PROPS -- material constants entered in the *USER MA TERIAL C optionC NPROPS -- number of material constantsCC----- This subroutine provides the plastic constitutive relation ofC single crystals for finite element code ABAQUS. The plastic slipC of single crystal obeys the Schmid law. The program gives theC choice of small deformation theory and theory of finite rotationC and finite strain.C The strain increment is composed of elastic part and plasticC part. The elastic strain increment corresponds to latticeC stretching, the plastic part is the sum over all slip systems ofC plastic slip. The shear strain increment for each slip system isC assumed a function of the ratio of corresponding resolved shearC stress over current strength, and of the time step. The resolvedC shear stress is the double product of stress tensor with the slipC deformation tensor (Schmid factor), and the increment of currentC strength is related to shear strain increments over all slipC systems through self- and latent-hardening functions.C----- The implicit integration method proposed by Peirce, Shih andC Needleman (1984) is used here. The subroutine provides an option C of iteration to solve stresses and solution dependent stateC variables within each increment.C----- The present program is for a single CUBIC crystal. However,C this code can be generalized for other crystals (e.g. HCP,C Tetragonal, Orthotropic, etc.). Only subroutines ROTATION andC SLIPSYS need to be modified to include the effect of crystalC aspect ratio.CC----- Important notice:CC (1) The number of state variables NSTATV must be larger than (or CFIX equal to) TEN (10) times the total number of slip systems inC all sets, NSLPTL, plus FIVE (5)CFIX NSTATV >= 10 * NSLPTL + 5C Denote s as a slip direction and m as normal to a slip plane.C Here (s,-m), (-s,m) and (-s,-m) are NOT consideredC independent of (s,m). The number of slip systems in each setC could be either 6, 12, 24 or 48 for a cubic crystal, e.g. 12C for {110}<111>.CC Users who need more parameters to characterize theC constitutive law of single crystal, e.g. the frameworkC proposed by Zarka, should make NSTATV larger than (or equal C to) the number of those parameters NPARMT plus nine timesC the total number of slip systems, NSLPTL, plus fiveCFIX NSTATV >= NPARMT + 10 * NSLPTL + 5CC (2) The tangent stiffness matrix in general is not symmetric ifC latent hardening is considered. Users must declare "UNSYMM"C in the input file, at the *USER MATERIAL card.CPARAMETER (ND=150)C----- The parameter ND determines the dimensions of the arrays inC this subroutine. The current choice 150 is a upper bound for aC cubic crystal with up to three sets of slip systems activated.C Users may reduce the parameter ND to any number as long as larger C than or equal to the total number of slip systems in all sets.C For example, if {110}<111> is the only set of slip systemC potentially activated, ND could be taken as twelve (12).cinclude 'aba_param.inc'cCHARACTER*8 CMNAMEEXTERNAL Fdimension stress(ntens),statev(nstatv),1 ddsdde(ntens,ntens),ddsddt(ntens),drplde(ntens),2 stran(ntens),dstran(ntens),time(2),predef(1),dpred(1),3 props(nprops),coords(3),drot(3,3),dfgrd0(3,3),dfgrd1(3,3)DIMENSION ISPDIR(3), ISPNOR(3), NSLIP(3),2 SLPDIR(3,ND), SLPNOR(3,ND), SLPDEF(6,ND),3 SLPSPN(3,ND), DSPDIR(3,ND), DSPNOR(3,ND),4 DLOCAL(6,6), D(6,6), ROTD(6,6), ROTATE(3,3),5 FSLIP(ND), DFDXSP(ND), DDEMSD(6,ND),6 H(ND,ND), DDGDDE(ND,6),7 DSTRES(6), DELATS(6), DSPIN(3), DVGRAD(3,3),8 DGAMMA(ND), DTAUSP(ND), DGSLIP(ND),9 WORKST(ND,ND), INDX(ND), TERM(3,3), TRM0(3,3), ITRM(3)DIMENSION FSLIP1(ND), STRES1(6), GAMMA1(ND), TAUSP1(ND),2 GSLP1(ND), SPNOR1(3,ND), SPDIR1(3,ND), DDSDE1(6,6),3 DSOLD(6), DGAMOD(ND), DTAUOD(ND), DGSPOD(ND),4 DSPNRO(3,ND), DSPDRO(3,ND),5 DHDGDG(ND,ND)C----- NSLIP -- number of slip systems in each setC----- SLPDIR -- slip directions (unit vectors in the initial state)C----- SLPNOR -- normals to slip planes (unit normals in the initialC state)C----- SLPDEF -- slip deformation tensors (Schmid factors)C SLPDEF(1,i) -- SLPDIR(1,i)*SLPNOR(1,i)C SLPDEF(2,i) -- SLPDIR(2,i)*SLPNOR(2,i)C SLPDEF(3,i) -- SLPDIR(3,i)*SLPNOR(3,i)C SLPDEF(4,i) -- SLPDIR(1,i)*SLPNOR(2,i)+C SLPDIR(2,i)*SLPNOR(1,i)C SLPDEF(5,i) -- SLPDIR(1,i)*SLPNOR(3,i)+C SLPDIR(3,i)*SLPNOR(1,i)C SLPDEF(6,i) -- SLPDIR(2,i)*SLPNOR(3,i)+C SLPDIR(3,i)*SLPNOR(2,i)C where index i corresponds to the ith slip systemC----- SLPSPN -- slip spin tensors (only needed for finite rotation)C SLPSPN(1,i) -- [SLPDIR(1,i)*SLPNOR(2,i)-C SLPDIR(2,i)*SLPNOR(1,i)]/2C SLPSPN(2,i) -- [SLPDIR(3,i)*SLPNOR(1,i)-C SLPDIR(1,i)*SLPNOR(3,i)]/2C SLPSPN(3,i) -- [SLPDIR(2,i)*SLPNOR(3,i)-C SLPDIR(3,i)*SLPNOR(2,i)]/2C where index i corresponds to the ith slip systemC----- DSPDIR -- increments of slip directionsC----- DSPNOR -- increments of normals to slip planesCC----- DLOCAL -- elastic matrix in local cubic crystal systemC----- D -- elastic matrix in global systemC----- ROTD -- rotation matrix transforming DLOCAL to DCC----- ROTATE -- rotation matrix, direction cosines of [100], [010]C and [001] of cubic crystal in global systemCC----- FSLIP -- shear strain-rates in slip systemsC----- DFDXSP -- derivatives of FSLIP w.r.t x=TAUSLP/GSLIP, whereC TAUSLP is the resolved shear stress and GSLIP is the C current strengthCC----- DDEMSD -- double dot product of the elastic moduli tensor withC the slip deformation tensor plus, only for finiteC rotation, the dot product of slip spin tensor withC the stressCC----- H -- self- and latent-hardening matrixC H(i,i) -- self hardening modulus of the ith slipC system (no sum over i)C H(i,j) -- latent hardening molulus of the ith slipC system due to a slip in the jth slip system C (i not equal j)CC----- DDGDDE -- derivatice of the shear strain increments in slipC systems w.r.t. the increment of strainsCC----- DSTRES -- Jaumann increments of stresses, i.e. corotationalC stress-increments formed on axes spinning with theC materialC----- DELATS -- strain-increments associated with lattice stretchingC DELATS(1) - DELATS(3) -- normal strain increments C DELATS(4) - DELATS(6) -- engineering shear strain C incrementsC----- DSPIN -- spin-increments associated with the material elementC DSPIN(1) -- component 12 of the spin tensorC DSPIN(2) -- component 31 of the spin tensorC DSPIN(3) -- component 23 of the spin tensorCC----- DVGRAD -- increments of deformation gradient in the currentC state, i.e. velocity gradient times the increment ofC timeCC----- DGAMMA -- increment of shear strains in slip systemsC----- DTAUSP -- increment of resolved shear stresses in slip systemsC----- DGSLIP -- increment of current strengths in slip systemsCCC----- Arrays for iteration:CC FSLIP1, STRES1, GAMMA1, TAUSP1, GSLP1 , SPNOR1, SPDIR1,C DDSDE1, DSOLD , DGAMOD, DTAUOD, DGSPOD, DSPNRO, DSPDRO, C DHDGDGCCC----- Solution dependent state variable STATEV:C Denote the number of total slip systems by NSLPTL, whichC will be calculated in this code.CC Array STATEV:C 1 - NSLPTL : current strength in slip systemsC NSLPTL+1 - 2*NSLPTL : shear strain in slip systemsC 2*NSLPTL+1 - 3*NSLPTL : resolved shear stress in slip systemsCC 3*NSLPTL+1 - 6*NSLPTL : current components of normals to slipC slip planesC 6*NSLPTL+1 - 9*NSLPTL : current components of slip directionsCCFIX 9*NSLPTL+1 - 10*NSLPTL : total cumulative shear strain on eachCFIX slip system (sum of the absoluteCFIX values of shear strains in each slipCFIX system individually)CFIXCFIX 10*NSLPTL+1 : total cumulative shear strain on allC slip systems (sum of the absoluteC values of shear strains in all slipC systems)CCFIX 10*NSLPTL+2 - NSTA TV-4 : additional parameters users may needC to characterize the constitutive lawC of a single crystal (if there areC any).CC NSTATV-3 : number of slip systems in the 1st setC NSTATV-2 : number of slip systems in the 2nd setC NSTATV-1 : number of slip systems in the 3rd setC NSTATV : total number of slip systems in allC setsCCC----- Material constants PROPS:CC PROPS(1) - PROPS(21) -- elastic constants for a general elasticC anisotropic materialCC isotropic : PROPS(i)=0 for i>2C PROPS(1) -- Young's modulusC PROPS(2) -- Poisson's ratioCC cubic : PROPS(i)=0 for i>3C PROPS(1) -- c11C PROPS(2) -- c12C PROPS(3) -- c44CC orthotropic : PORPS(i)=0 for i>9C PROPS(1) - PROPS(9) are D1111, D1122, D2222, C D1133, D2233, D3333, D1212, D1313, D2323,C respectively, which has the same definitionC as ABAQUS for orthotropic materialsC (see *ELASTIC card)CC anisotropic : PROPS(1) - PROPS(21) are D1111, D1122,C D2222, D1133, D2233, D3333, D1112, D2212,C D3312, D1212, D1113, D2213, D3313, D1213,C D1313, D1123, D2223, D3323, D1223, D1323,C D2323, respectively, which has the sameC definition as ABAQUS for anisotropicC materials (see *ELASTIC card)CCC PROPS(25) - PROPS(56) -- parameters characterizing all slipC systems to be activated in a cubicC crystalCC PROPS(25) -- number of sets of slip systems (maximum 3),C e.g. (110)[1-11] and (101)[11-1] are in theC same set of slip systems, (110)[1-11] andC (121)[1-11] belong to different sets of slipC systemsC (It must be a real number, e.g. 3., not 3 !)CC PROPS(33) - PROPS(35) -- normal to a typical slip plane inC the first set of slip systems,C e.g. (1 1 0)C (They must be real numbers, e.g.C 1. 1. 0., not 1 1 0 !)C PROPS(36) - PROPS(38) -- a typical slip direction in theC first set of slip systems, e.g.C [1 1 1]C (They must be real numbers, e.g.C 1. 1. 1., not 1 1 1 !)CC PROPS(41) - PROPS(43) -- normal to a typical slip plane inC the second set of slip systemsC (real numbers)C PROPS(44) - PROPS(46) -- a typical slip direction in theC second set of slip systemsC (real numbers)CC PROPS(49) - PROPS(51) -- normal to a typical slip plane inC the third set of slip systemsC (real numbers)C PROPS(52) - PROPS(54) -- a typical slip direction in theC third set of slip systemsC (real numbers)CCC PROPS(57) - PROPS(72) -- parameters characterizing the initialC orientation of a single crystal inC global systemC The directions in global system and directions in localC cubic crystal system of two nonparallel vectors are neededC to determine the crystal orientation.CC PROPS(57) - PROPS(59) -- [p1 p2 p3], direction of firstC vector in local cubic crystalC system, e.g. [1 1 0]C (They must be real numbers, e.g.C 1. 1. 0., not 1 1 0 !)C PROPS(60) - PROPS(62) -- [P1 P2 P3], direction of firstC vector in global system, e.g.C [2. 1. 0.]C (It does not have to be a unitC vector)CC PROPS(65) - PROPS(67) -- direction of second vector inC local cubic crystal system (real C numbers)C PROPS(68) - PROPS(70) -- direction of second vector inC global systemCCC PROPS(73) - PROPS(96) -- parameters characterizing the visco-C plastic constitutive law (shearC strain-rate vs. resolved shearC stress), e.g. a power-law relationCC PROPS(73) - PROPS(80) -- parameters for the first set ofC slip systemsC PROPS(81) - PROPS(88) -- parameters for the second set of C slip systemsC PROPS(89) - PROPS(96) -- parameters for the third set ofC slip systemsCCC PROPS(97) - PROPS(144)-- parameters characterizing the self-C and latent-hardening laws of slipC systemsCC PROPS(97) - PROPS(104)-- self-hardening parameters for the C first set of slip systemsC PROPS(105)- PROPS(112)-- latent-hardening parameters for C the first set of slip systems and C interaction with other sets ofC slip systemsCC PROPS(113)- PROPS(120)-- self-hardening parameters for the C second set of slip systemsC PROPS(121)- PROPS(128)-- latent-hardening parameters for C the second set of slip systems C and interaction with other sets C of slip systemsCC PROPS(129)- PROPS(136)-- self-hardening parameters for the C third set of slip systemsC PROPS(137)- PROPS(144)-- latent-hardening parameters for C the third set of slip systems and C interaction with other sets ofC slip systemsCCC PROPS(145)- PROPS(152)-- parameters characterizing forward time C integration scheme and finiteC deformationCC PROPS(145) -- parameter theta controlling the implicitC integration, which is between 0 and 1C 0. : explicit integrationC 0.5 : recommended valueC 1. : fully implicit integrationCC PROPS(146) -- parameter NLGEOM controlling whether the C effect of finite rotation and finite strainC of crystal is considered,C 0. : small deformation theoryC otherwise : theory of finite rotation andC finite strainCCC PROPS(153)- PROPS(160)-- parameters characterizing iterationC methodCC PROPS(153) -- parameter ITRATN controlling whether theC iteration method is used,C 0. : no iterationC otherwise : iterationCC PROPS(154) -- maximum number of iteration ITRMAXCC PROPS(155) -- absolute error of shear strains in slipC systems GAMERRCC----- Elastic matrix in local cubic crystal system: DLOCALDO J=1,6DO I=1,6DLOCAL(I,J)=0.END DOEND DOCHECK=0.DO J=10,21CHECK=CHECK+ABS(PROPS(J))END DOIF (CHECK.EQ.0.) THENDO J=4,9CHECK=CHECK+ABS(PROPS(J))END DOIF (CHECK.EQ.0.) THENIF (PROPS(3).EQ.0.) THENC----- Isotropic materialGSHEAR=PROPS(1)/2./(1.+PROPS(2))E11=2.*GSHEAR*(1.-PROPS(2))/(1.-2.*PROPS(2))E12=2.*GSHEAR*PROPS(2)/(1.-2.*PROPS(2))DO J=1,3DLOCAL(J,J)=E11DO I=1,3IF (I.NE.J) DLOCAL(I,J)=E12END DODLOCAL(J+3,J+3)=GSHEAREND DOELSEC----- Cubic materialDO J=1,3DLOCAL(J,J)=PROPS(1)DO I=1,3IF (I.NE.J) DLOCAL(I,J)=PROPS(2)END DODLOCAL(J+3,J+3)=PROPS(3)END DOEND IFELSEC----- Orthotropic metarialDLOCAL(1,1)=PROPS(1)DLOCAL(1,2)=PROPS(2)DLOCAL(2,1)=PROPS(2)DLOCAL(2,2)=PROPS(3)DLOCAL(1,3)=PROPS(4)DLOCAL(3,1)=PROPS(4)DLOCAL(2,3)=PROPS(5)DLOCAL(3,2)=PROPS(5)DLOCAL(3,3)=PROPS(6)DLOCAL(4,4)=PROPS(7)DLOCAL(5,5)=PROPS(8)DLOCAL(6,6)=PROPS(9)END IFELSEC----- General anisotropic materialID=0DO J=1,6DO I=1,JID=ID+1DLOCAL(I,J)=PROPS(ID)DLOCAL(J,I)=DLOCAL(I,J)END DOEND DOEND IFC----- Rotation matrix: ROTA TE, i.e. direction cosines of [100], [010]C and [001] of a cubic crystal in global systemCCALL ROTATION (PROPS(57), ROTA TE)C----- Rotation matrix: ROTD to transform local elastic matrix DLOCAL C to global elastic matrix DCDO J=1,3J1=1+J/3J2=2+J/2DO I=1,3I1=1+I/3I2=2+I/2ROTD(I,J)=ROTATE(I,J)**2ROTD(I,J+3)=2.*ROTATE(I,J1)*ROTA TE(I,J2)ROTD(I+3,J)=ROTA TE(I1,J)*ROTATE(I2,J)ROTD(I+3,J+3)=ROTA TE(I1,J1)*ROTA TE(I2,J2)+2 ROTA TE(I1,J2)*ROTA TE(I2,J1)END DOEND DOC----- Elastic matrix in global system: DC {D} = {ROTD} * {DLOCAL} * {ROTD}transposeCDO J=1,6DO I=1,6D(I,J)=0.END DOEND DODO J=1,6DO I=1,JDO K=1,6DO L=1,6D(I,J)=D(I,J)+DLOCAL(K,L)*ROTD(I,K)*ROTD(J,L) END DOEND DOD(J,I)=D(I,J)END DOEND DOC----- Total number of sets of slip systems: NSETNSET=NINT(PROPS(25))IF (NSET.LT.1) THENWRITE (6,*) '***ERROR - zero sets of slip systems'STOPELSE IF (NSET.GT.3) THENWRITE (6,*)2 '***ERROR - more than three sets of slip systems'STOPEND IFC----- Implicit integration parameter: THETATHETA=PROPS(145)C----- Finite deformation ?C----- NLGEOM = 0, small deformation theoryC otherwise, theory of finite rotation and finite strain, Users C must declare "NLGEOM" in the input file, at the *STEP card CIF (PROPS(146).EQ.0.) THENNLGEOM=0ELSENLGEOM=1END IFC----- Iteration?C----- ITRATN = 0, no iterationC otherwise, iteration (solving increments of stresses andC solution dependent state variables)CIF (PROPS(153).EQ.0.) THENITRATN=0ELSEITRATN=1END IFITRMAX=NINT(PROPS(154))GAMERR=PROPS(155)NITRTN=-1DO I=1,NTENSDSOLD(I)=0.END DODO J=1,NDDGAMOD(J)=0.DTAUOD(J)=0.DGSPOD(J)=0.DO I=1,3DSPNRO(I,J)=0.DSPDRO(I,J)=0.END DOEND DOC----- Increment of spin associated with the material element: DSPIN C (only needed for finite rotation)CIF (NLGEOM.NE.0) THENDO J=1,3DO I=1,3TERM(I,J)=DROT(J,I)TRM0(I,J)=DROT(J,I)END DOTERM(J,J)=TERM(J,J)+1.D0TRM0(J,J)=TRM0(J,J)-1.D0END DOCALL LUDCMP (TERM, 3, 3, ITRM, DDCMP)DO J=1,3CALL LUBKSB (TERM, 3, 3, ITRM, TRM0(1,J)) END DODSPIN(1)=TRM0(2,1)-TRM0(1,2)DSPIN(2)=TRM0(1,3)-TRM0(3,1)DSPIN(3)=TRM0(3,2)-TRM0(2,3)END IFC----- Increment of dilatational strain: DEVDEV=0.D0DO I=1,NDIDEV=DEV+DSTRAN(I)END DOC----- Iteration starts (only when iteration method is used)1000 CONTINUEC----- Parameter NITRTN: number of iterationsC NITRTN = 0 --- no-iteration solutionCNITRTN=NITRTN+1C----- Check whether the current stress state is the initial stateIF (STATEV(1).EQ.0.) THENC----- Initial stateCC----- Generating the following parameters and variables at initialC state:C Total number of slip systems in all the sets NSLPTLC Number of slip systems in each set NSLIPC Unit vectors in initial slip directions SLPDIRC Unit normals to initial slip planes SLPNORCNSLPTL=0DO I=1,NSETISPNOR(1)=NINT(PROPS(25+8*I))ISPNOR(2)=NINT(PROPS(26+8*I))ISPNOR(3)=NINT(PROPS(27+8*I))ISPDIR(1)=NINT(PROPS(28+8*I))ISPDIR(2)=NINT(PROPS(29+8*I))ISPDIR(3)=NINT(PROPS(30+8*I))CALL SLIPSYS (ISPDIR, ISPNOR, NSLIP(I), SLPDIR(1,NSLPTL+1),2 SLPNOR(1,NSLPTL+1), ROTATE)NSLPTL=NSLPTL+NSLIP(I)END DOIF (ND.LT.NSLPTL) THENWRITE (6,*)2 '***ERROR - parameter ND chosen by the present user is less than3 the total number of slip systems NSLPTL'STOPEND IFC----- Slip deformation tensor: SLPDEF (Schmid factors)DO J=1,NSLPTLSLPDEF(1,J)=SLPDIR(1,J)*SLPNOR(1,J)SLPDEF(2,J)=SLPDIR(2,J)*SLPNOR(2,J)SLPDEF(3,J)=SLPDIR(3,J)*SLPNOR(3,J)SLPDEF(4,J)=SLPDIR(1,J)*SLPNOR(2,J)+SLPDIR(2,J)*SLPNOR(1,J)SLPDEF(5,J)=SLPDIR(1,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(1,J)SLPDEF(6,J)=SLPDIR(2,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(2,J) END DOC----- Initial value of state variables: unit normal to a slip planeC and unit vector in a slip directionCSTATEV(NSTATV)=FLOAT(NSLPTL)DO I=1,NSETSTA TEV(NSTA TV-4+I)=FLOAT(NSLIP(I))END DOIDNOR=3*NSLPTLIDDIR=6*NSLPTLDO J=1,NSLPTLDO I=1,3IDNOR=IDNOR+1STA TEV(IDNOR)=SLPNOR(I,J)IDDIR=IDDIR+1STA TEV(IDDIR)=SLPDIR(I,J)END DOEND DOC----- Initial value of the current strength for all slip systemsCCALL GSLPINIT (STATEV(1), NSLIP, NSLPTL, NSET, PROPS(97))C----- Initial value of shear strain in slip systemsCFIX-- Initial value of cumulative shear strain in each slip systemsDO I=1,NSLPTLSTA TEV(NSLPTL+I)=0.CFIXASTA TEV(9*NSLPTL+I)=0.CFIXBEND DOCFIXASTATEV(10*NSLPTL+1)=0.CFIXBC----- Initial value of the resolved shear stress in slip systemsDO I=1,NSLPTLTERM1=0.DO J=1,NTENSIF (J.LE.NDI) THENTERM1=TERM1+SLPDEF(J,I)*STRESS(J)ELSETERM1=TERM1+SLPDEF(J-NDI+3,I)*STRESS(J) END IFEND DOSTA TEV(2*NSLPTL+I)=TERM1END DOELSEC----- Current stress stateCC----- Copying from the array of state variables STA TVE the following C parameters and variables at current stress state:C Total number of slip systems in all the sets NSLPTLC Number of slip systems in each set NSLIPC Current slip directions SLPDIRC Normals to current slip planes SLPNORCNSLPTL=NINT(STA TEV(NSTATV))DO I=1,NSETNSLIP(I)=NINT(STATEV(NSTATV-4+I))END DOIDNOR=3*NSLPTLIDDIR=6*NSLPTLDO J=1,NSLPTLDO I=1,3IDNOR=IDNOR+1SLPNOR(I,J)=STATEV(IDNOR)IDDIR=IDDIR+1SLPDIR(I,J)=STA TEV(IDDIR)END DOEND DOC----- Slip deformation tensor: SLPDEF (Schmid factors)DO J=1,NSLPTLSLPDEF(1,J)=SLPDIR(1,J)*SLPNOR(1,J)SLPDEF(2,J)=SLPDIR(2,J)*SLPNOR(2,J)SLPDEF(3,J)=SLPDIR(3,J)*SLPNOR(3,J)SLPDEF(4,J)=SLPDIR(1,J)*SLPNOR(2,J)+SLPDIR(2,J)*SLPNOR(1,J)SLPDEF(5,J)=SLPDIR(1,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(1,J)SLPDEF(6,J)=SLPDIR(2,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(2,J) END DOEND IFC----- Slip spin tensor: SLPSPN (only needed for finite rotation)IF (NLGEOM.NE.0) THENDO J=1,NSLPTLSLPSPN(1,J)=0.5*(SLPDIR(1,J)*SLPNOR(2,J)-2 SLPDIR(2,J)*SLPNOR(1,J))SLPSPN(2,J)=0.5*(SLPDIR(3,J)*SLPNOR(1,J)-2 SLPDIR(1,J)*SLPNOR(3,J))SLPSPN(3,J)=0.5*(SLPDIR(2,J)*SLPNOR(3,J)-2 SLPDIR(3,J)*SLPNOR(2,J))END DOEND IFC----- Double dot product of elastic moduli tensor with the slipC deformation tensor (Schmid factors) plus, only for finiteC rotation, the dot product of slip spin tensor with the stress:C DDEMSDCDO J=1,NSLPTLDO I=1,6DDEMSD(I,J)=0.DO K=1,6DDEMSD(I,J)=DDEMSD(I,J)+D(K,I)*SLPDEF(K,J)END DOEND DOEND DOIF (NLGEOM.NE.0) THENDO J=1,NSLPTLDDEMSD(4,J)=DDEMSD(4,J)-SLPSPN(1,J)*STRESS(1)DDEMSD(5,J)=DDEMSD(5,J)+SLPSPN(2,J)*STRESS(1)IF (NDI.GT.1) THENDDEMSD(4,J)=DDEMSD(4,J)+SLPSPN(1,J)*STRESS(2)DDEMSD(6,J)=DDEMSD(6,J)-SLPSPN(3,J)*STRESS(2)END IFIF (NDI.GT.2) THENDDEMSD(5,J)=DDEMSD(5,J)-SLPSPN(2,J)*STRESS(3)DDEMSD(6,J)=DDEMSD(6,J)+SLPSPN(3,J)*STRESS(3)END IFIF (NSHR.GE.1) THENDDEMSD(1,J)=DDEMSD(1,J)+SLPSPN(1,J)*STRESS(NDI+1)DDEMSD(2,J)=DDEMSD(2,J)-SLPSPN(1,J)*STRESS(NDI+1)DDEMSD(5,J)=DDEMSD(5,J)-SLPSPN(3,J)*STRESS(NDI+1)DDEMSD(6,J)=DDEMSD(6,J)+SLPSPN(2,J)*STRESS(NDI+1) END IFIF (NSHR.GE.2) THENDDEMSD(1,J)=DDEMSD(1,J)-SLPSPN(2,J)*STRESS(NDI+2)DDEMSD(3,J)=DDEMSD(3,J)+SLPSPN(2,J)*STRESS(NDI+2)DDEMSD(4,J)=DDEMSD(4,J)+SLPSPN(3,J)*STRESS(NDI+2)DDEMSD(6,J)=DDEMSD(6,J)-SLPSPN(1,J)*STRESS(NDI+2) END IFIF (NSHR.EQ.3) THENDDEMSD(2,J)=DDEMSD(2,J)+SLPSPN(3,J)*STRESS(NDI+3)DDEMSD(3,J)=DDEMSD(3,J)-SLPSPN(3,J)*STRESS(NDI+3)DDEMSD(4,J)=DDEMSD(4,J)-SLPSPN(2,J)*STRESS(NDI+3)DDEMSD(5,J)=DDEMSD(5,J)+SLPSPN(1,J)*STRESS(NDI+3) END IFEND DOEND IFC----- Shear strain-rate in a slip system at the start of increment:C FSLIP, and its derivative: DFDXSPCID=1DO I=1,NSETIF (I.GT.1) ID=ID+NSLIP(I-1)CALL STRAINRATE (STATEV(NSLPTL+ID), STATEV(2*NSLPTL+ID),2 STA TEV(ID), NSLIP(I), FSLIP(ID), DFDXSP(ID),3 PROPS(65+8*I))END DOC----- Self- and latent-hardening laws。

Advances in

Advances in Geosciences,4,17–22,2005 SRef-ID:1680-7359/adgeo/2005-4-17 European Geosciences Union©2005Author(s).This work is licensed under a Creative CommonsLicense.Advances in GeosciencesIncorporating level set methods in Geographical Information Systems(GIS)for land-surface process modelingD.PullarGeography Planning and Architecture,The University of Queensland,Brisbane QLD4072,Australia Received:1August2004–Revised:1November2004–Accepted:15November2004–Published:9August2005nd-surface processes include a broad class of models that operate at a landscape scale.Current modelling approaches tend to be specialised towards one type of pro-cess,yet it is the interaction of processes that is increasing seen as important to obtain a more integrated approach to land management.This paper presents a technique and a tool that may be applied generically to landscape processes. The technique tracks moving interfaces across landscapes for processes such as waterﬂow,biochemical diffusion,and plant dispersal.Its theoretical development applies a La-grangian approach to motion over a Eulerian grid space by tracking quantities across a landscape as an evolving front. An algorithm for this technique,called level set method,is implemented in a geographical information system(GIS).It ﬁts with aﬁeld data model in GIS and is implemented as operators in map algebra.The paper describes an implemen-tation of the level set methods in a map algebra program-ming language,called MapScript,and gives example pro-gram scripts for applications in ecology and hydrology.1IntroductionOver the past decade there has been an explosion in the ap-plication of models to solve environmental issues.Many of these models are speciﬁc to one physical process and of-ten require expert knowledge to use.Increasingly generic modeling frameworks are being sought to provide analyti-cal tools to examine and resolve complex environmental and natural resource problems.These systems consider a vari-ety of land condition characteristics,interactions and driv-ing physical processes.Variables accounted for include cli-mate,topography,soils,geology,land cover,vegetation and hydro-geography(Moore et al.,1993).Physical interactions include processes for climatology,hydrology,topographic landsurface/sub-surfaceﬂuxes and biological/ecological sys-Correspondence to:D.Pullar(d.pullar@.au)tems(Sklar and Costanza,1991).Progress has been made in linking model-speciﬁc systems with tools used by environ-mental managers,for instance geographical information sys-tems(GIS).While this approach,commonly referred to as loose coupling,provides a practical solution it still does not improve the scientiﬁc foundation of these models nor their integration with other models and related systems,such as decision support systems(Argent,2003).The alternative ap-proach is for tightly coupled systems which build functional-ity into a system or interface to domain libraries from which a user may build custom solutions using a macro language or program scripts.The approach supports integrated models through interface speciﬁcations which articulate the funda-mental assumptions and simpliﬁcations within these models. The problem is that there are no environmental modelling systems which are widely used by engineers and scientists that offer this level of interoperability,and the more com-monly used GIS systems do not currently support space and time representations and operations suitable for modelling environmental processes(Burrough,1998)(Sui and Magio, 1999).Providing a generic environmental modeling framework for practical environmental issues is challenging.It does not exist now despite an overwhelming demand because there are deep technical challenges to build integrated modeling frameworks in a scientiﬁcally rigorous manner.It is this chal-lenge this research addresses.1.1Background for ApproachThe paper describes a generic environmental modeling lan-guage integrated with a Geographical Information System (GIS)which supports spatial-temporal operators to model physical interactions occurring in two ways.The trivial case where interactions are isolated to a location,and the more common and complex case where interactions propa-gate spatially across landscape surfaces.The programming language has a strong theoretical and algorithmic basis.The-oretically,it assumes a Eulerian representation of state space,Fig.1.Shows a)a propagating interface parameterised by differ-ential equations,b)interface fronts have variable intensity and may expand or contract based onﬁeld gradients and driving process. but propagates quantities across landscapes using Lagrangian equations of motion.In physics,a Lagrangian view focuses on how a quantity(water volume or particle)moves through space,whereas an Eulerian view focuses on a localﬁxed area of space and accounts for quantities moving through it.The beneﬁt of this approach is that an Eulerian perspective is em-inently suited to representing the variation of environmen-tal phenomena across space,but it is difﬁcult to conceptu-alise solutions for the equations of motion and has compu-tational drawbacks(Press et al.,1992).On the other hand, the Lagrangian view is often not favoured because it requires a global solution that makes it difﬁcult to account for local variations,but has the advantage of solving equations of mo-tion in an intuitive and numerically direct way.The research will address this dilemma by adopting a novel approach from the image processing discipline that uses a Lagrangian ap-proach over an Eulerian grid.The approach,called level set methods,provides an efﬁcient algorithm for modeling a natural advancing front in a host of settings(Sethian,1999). The reason the method works well over other approaches is that the advancing front is described by equations of motion (Lagrangian view),but computationally the front propagates over a vectorﬁeld(Eulerian view).Hence,we have a very generic way to describe the motion of quantities,but can ex-plicitly solve their advancing properties locally as propagat-ing zones.The research work will adapt this technique for modeling the motion of environmental variables across time and space.Speciﬁcally,it will add new data models and op-erators to a geographical information system(GIS)for envi-ronmental modeling.This is considered to be a signiﬁcant research imperative in spatial information science and tech-nology(Goodchild,2001).The main focus of this paper is to evaluate if the level set method(Sethian,1999)can:–provide a theoretically and empirically supportable methodology for modeling a range of integral landscape processes,–provide an algorithmic solution that is not sensitive to process timing,is computationally stable and efﬁcient as compared to conventional explicit solutions to diffu-sive processes models,–be developed as part of a generic modelling language in GIS to express integrated models for natural resource and environmental problems?The outline for the paper is as follow.The next section will describe the theory for spatial-temporal processing us-ing level sets.Section3describes how this is implemented in a map algebra programming language.Two application examples are given–an ecological and a hydrological ex-ample–to demonstrate the use of operators for computing reactive-diffusive interactions in landscapes.Section4sum-marises the contribution of this research.2Theory2.1IntroductionLevel set methods(Sethian,1999)have been applied in a large collection of applications including,physics,chemistry,ﬂuid dynamics,combustion,material science,fabrication of microelectronics,and computer vision.Level set methods compute an advancing interface using an Eulerian grid and the Lagrangian equations of motion.They are similar to cost distance modeling used in GIS(Burroughs and McDonnell, 1998)in that they compute the spread of a variable across space,but the motion is based upon partial differential equa-tions related to the physical process.The advancement of the interface is computed through time along a spatial gradient, and it may expand or contract in its extent.See Fig.1.2.2TheoryThe advantage of the level set method is that it models mo-tion along a state-space gradient.Level set methods start with the equation of motion,i.e.an advancing front with velocity F is characterised by an arrival surface T(x,y).Note that F is a velocityﬁeld in a spatial sense.If F was constant this would result in an expanding series of circular fronts,but for different values in a velocityﬁeld the front will have a more contorted appearance as shown in Fig.1b.The motion of thisinterface is always normal to the interface boundary,and its progress is regulated by several factors:F=f(L,G,I)(1)where L=local properties that determine the shape of advanc-ing front,G=global properties related to governing forces for its motion,I=independent properties that regulate and inﬂu-ence the motion.If the advancing front is modeled strictly in terms of the movement of entity particles,then a straightfor-ward velocity equation describes its motion:|∇T|F=1given T0=0(2) where the arrival function T(x,y)is a travel cost surface,and T0is the initial position of the interface.Instead we use level sets to describe the interface as a complex function.The level set functionφis an evolving front consistent with the under-lying viscosity solution deﬁned by partial differential equa-tions.This is expressed by the equation:φt+F|∇φ|=0givenφ(x,y,t=0)(3)whereφt is a complex interface function over time period 0..n,i.e.φ(x,y,t)=t0..tn,∇φis the spatial and temporal derivatives for viscosity equations.The Eulerian view over a spatial domain imposes a discretisation of space,i.e.the raster grid,which records changes in value z.Hence,the level set function becomesφ(x,y,z,t)to describe an evolv-ing surface over time.Further details are given in Sethian (1999)along with efﬁcient algorithms.The next section de-scribes the integration of the level set methods with GIS.3Map algebra modelling3.1Map algebraSpatial models are written in a map algebra programming language.Map algebra is a function-oriented language that operates on four implicit spatial data types:point,neighbour-hood,zonal and whole landscape surfaces.Surfaces are typ-ically represented as a discrete raster where a point is a cell, a neighbourhood is a kernel centred on a cell,and zones are groups of mon examples of raster data include ter-rain models,categorical land cover maps,and scalar temper-ature surfaces.Map algebra is used to program many types of landscape models ranging from land suitability models to mineral exploration in the geosciences(Burrough and Mc-Donnell,1998;Bonham-Carter,1994).The syntax for map algebra follows a mathematical style with statements expressed as equations.These equations use operators to manipulate spatial data types for point and neighbourhoods.Expressions that manipulate a raster sur-face may use a global operation or alternatively iterate over the cells in a raster.For instance the GRID map algebra (Gao et al.,1993)deﬁnes an iteration construct,called do-cell,to apply equations on a cell-by-cell basis.This is triv-ially performed on columns and rows in a clockwork manner. However,for environmental phenomena there aresituations Fig.2.Spatial processing orders for raster.where the order of computations has a special signiﬁcance. For instance,processes that involve spreading or transport acting along environmental gradients within the landscape. Therefore special control needs to be exercised on the order of execution.Burrough(1998)describes two extra control mechanisms for diffusion and directed topology.Figure2 shows the three principle types of processing orders,and they are:–row scan order governed by the clockwork lattice struc-ture,–spread order governed by the spreading or scattering ofa material from a more concentrated region,–ﬂow order governed by advection which is the transport of a material due to velocity.Our implementation of map algebra,called MapScript (Pullar,2001),includes a special iteration construct that sup-ports these processing orders.MapScript is a lightweight lan-guage for processing raster-based GIS data using map alge-bra.The language parser and engine are built as a software component to interoperate with the IDRISI GIS(Eastman, 1997).MapScript is built in C++with a class hierarchy based upon a value type.Variants for value types include numeri-cal,boolean,template,cells,or a grid.MapScript supports combinations of these data types within equations with basic arithmetic and relational comparison operators.Algebra op-erations on templates typically result in an aggregate value assigned to a cell(Pullar,2001);this is similar to the con-volution integral in image algebras(Ritter et al.,1990).The language supports iteration to execute a block of statements in three ways:a)docell construct to process raster in a row scan order,b)dospread construct to process raster in a spreadwhile(time<100)dospreadpop=pop+(diffuse(kernel*pop))pop=pop+(r*pop*dt*(1-(pop/K)) enddoendwhere the diffusive constant is stored in thekernel:Fig.3.Map algebra script and convolution kernel for population dispersion.The variable pop is a raster,r,K and D are constants, dt is the model time step,and the kernel is a3×3template.It is assumed a time step is deﬁned and the script is run in a simulation. Theﬁrst line contained in the nested cell processing construct(i.e. dospread)is the diffusive term and the second line is the population growth term.order,c)doﬂow to process raster byﬂow order.Examples are given in subsequent sections.Process models will also involve a timing loop which may be handled as a general while(<condition>)..end construct in MapScript where the condition expression includes a system time variable.This time variable is used in a speciﬁc fashion along with a system time step by certain operators,namely diffuse()andﬂuxﬂow() described in the next section,to model diffusion and advec-tion as a time evolving front.The evolving front represents quantities such as vegetation growth or surface runoff.3.2Ecological exampleThis section presents an ecological example based upon plant dispersal in a landscape.The population of a species follows a controlled growth rate and at the same time spreads across landscapes.The theory of the rate of spread of an organism is given in Tilman and Kareiva(1997).The area occupied by a species grows log-linear with time.This may be modelled by coupling a spatial diffusion term with an exponential pop-ulation growth term;the combination produces the familiar reaction-diffusion model.A simple growth population model is used where the reac-tion term considers one population controlled by births and mortalities is:dN dt =r·N1−NK(4)where N is the size of the population,r is the rate of change of population given in terms of the difference between birth and mortality rates,and K is the carrying capacity.Further dis-cussion of population models can be found in Jrgensen and Bendoricchio(2001).The diffusive term spreads a quantity through space at a speciﬁed rate:dudt=Dd2udx2(5) where u is the quantity which in our case is population size, and D is the diffusive coefﬁcient.The model is operated as a coupled computation.Over a discretized space,or raster,the diffusive term is estimated using a numerical scheme(Press et al.,1992).The distance over which diffusion takes place in time step dt is minimally constrained by the raster resolution.For a stable computa-tional process the following condition must be satisﬁed:2Ddtdx2≤1(6) This basically states that to account for the diffusive pro-cess,the term2D·dx is less than the velocity of the advancing front.This would not be difﬁcult to compute if D is constant, but is problematic if D is variable with respect to landscape conditions.This problem may be overcome by progressing along a diffusive front over the discrete raster based upon distance rather than being constrained by the cell resolution.The pro-cessing and diffusive operator is implemented in a map al-gebra programming language.The code fragment in Fig.3 shows a map algebra script for a single time step for the cou-pled reactive-diffusion model for population growth.The operator of interest in the script shown in Fig.3is the diffuse operator.It is assumed that the script is run with a given time step.The operator uses a system time step which is computed to balance the effect of process errors with efﬁ-cient computation.With knowledge of the time step the it-erative construct applies an appropriate distance propagation such that the condition in Eq.(3)is not violated.The level set algorithm(Sethian,1999)is used to do this in a stable and accurate way.As a diffusive front propagates through the raster,a cost distance kernel assigns the proper time to each raster cell.The time assigned to the cell corresponds to the minimal cost it takes to reach that cell.Hence cell pro-cessing is controlled by propagating the kernel outward at a speed adaptive to the local context rather than meeting an arbitrary global constraint.3.3Hydrological exampleThis section presents a hydrological example based upon sur-face dispersal of excess rainfall across the terrain.The move-ment of water is described by the continuity equation:∂h∂t=e t−∇·q t(7) where h is the water depth(m),e t is the rainfall excess(m/s), q is the discharge(m/hr)at time t.Discharge is assumed to have steady uniformﬂow conditions,and is determined by Manning’s equation:q t=v t h t=1nh5/3ts1/2(8)putation of current cell(x+ x,t,t+ ).where q t is theﬂow velocity(m/s),h t is water depth,and s is the surface slope(m/m).An explicit method of calcula-tion is used to compute velocity and depth over raster cells, and equations are solved at each time step.A conservative form of aﬁnite difference method solves for q t in Eq.(5). To simplify discussions we describe quasi-one-dimensional equations for theﬂow problem.The actual numerical com-putations are normally performed on an Eulerian grid(Julien et al.,1995).Finite-element approximations are made to solve the above partial differential equations for the one-dimensional case ofﬂow along a strip of unit width.This leads to a cou-pled model with one term to maintain the continuity ofﬂow and another term to compute theﬂow.In addition,all calcu-lations must progress from an uphill cell to the down slope cell.This is implemented in map algebra by a iteration con-struct,called doﬂow,which processes a raster byﬂow order. Flow distance is measured in cell size x per unit length. One strip is processed during a time interval t(Fig.4).The conservative solution for the continuity term using aﬁrst or-der approximation for Eq.(5)is derived as:h x+ x,t+ t=h x+ x,t−q x+ x,t−q x,txt(9)where the inﬂow q x,t and outﬂow q x+x,t are calculated in the second term using Equation6as:q x,t=v x,t·h t(10) The calculations approximate discharge from previous time interval.Discharge is dynamically determined within the continuity equation by water depth.The rate of change in state variables for Equation6needs to satisfy a stability condition where v· t/ x≤1to maintain numerical stabil-ity.The physical interpretation of this is that aﬁnite volume of water wouldﬂow across and out of a cell within the time step t.Typically the cell resolution isﬁxed for the raster, and adjusting the time step requires restarting the simulation while(time<120)doflow(dem)fvel=1/n*pow(depth,m)*sqrt(grade)depth=depth+(depth*fluxflow(fvel)) enddoendFig.5.Map algebra script for excess rainfallﬂow computed over a 120minute event.The variables depth and grade are rasters,fvel is theﬂow velocity,n and m are constants in Manning’s equation.It is assumed a time step is deﬁned and the script is run in a simulation. Theﬁrst line in the nested cell processing(i.e.doﬂow)computes theﬂow velocity and the second line computes the change in depth from the previous value plus any net change(inﬂow–outﬂow)due to velocityﬂux across the cell.cycle.Flow velocities change dramatically over the course of a storm event,and it is problematic to set an appropriate time step which is efﬁcient and yields a stable result.The hydrological model has been implemented in a map algebra programming language Pullar(2003).To overcome the problem mentioned above we have added high level oper-ators to compute theﬂow as an advancing front over a land-scape.The time step advances this front adaptively across the landscape based upon theﬂow velocity.The level set algorithm(Sethian,1999)is used to do this in a stable and accurate way.The map algebra script is given in Fig.5.The important operator is theﬂuxﬂow operator.It computes the advancing front for waterﬂow across a DEM by hydrologi-cal principles,and computes the local drainageﬂux rate for each cell.Theﬂux rate is used to compute the net change in a cell in terms ofﬂow depth over an adaptive time step.4ConclusionsThe paper has described an approach to extend the function-ality of tightly coupled environmental models in GIS(Ar-gent,2004).A long standing criticism of GIS has been its in-ability to handle dynamic spatial models.Other researchers have also addressed this issue(Burrough,1998).The con-tribution of this paper is to describe how level set methods are:i)an appropriate scientiﬁc basis,and ii)able to perform stable time-space computations for modelling landscape pro-cesses.The level set method provides the following beneﬁts:–it more directly models motion of spatial phenomena and may handle both expanding and contracting inter-faces,–is based upon differential equations related to the spatial dynamics of physical processes.Despite the potential for using level set methods in GIS and land-surface process modeling,there are no commercial or research systems that use this mercial sys-tems such as GRID(Gao et al.,1993),and research systems such as PCRaster(Wesseling et al.,1996)offerﬂexible andpowerful map algebra programming languages.But opera-tions that involve reaction-diffusive processing are speciﬁc to one context,such as groundwaterﬂow.We believe the level set method offers a more generic approach that allows a user to programﬂow and diffusive landscape processes for a variety of application contexts.We have shown that it pro-vides an appropriate theoretical underpinning and may be ef-ﬁciently implemented in a GIS.We have demonstrated its application for two landscape processes–albeit relatively simple examples–but these may be extended to deal with more complex and dynamic circumstances.The validation for improved environmental modeling tools ultimately rests in their uptake and usage by scientists and engineers.The tool may be accessed from the web site .au/projects/mapscript/(version with enhancements available April2005)for use with IDRSIS GIS(Eastman,1997)and in the future with ArcGIS. It is hoped that a larger community of users will make use of the methodology and implementation for a variety of environmental modeling applications.Edited by:P.Krause,S.Kralisch,and W.Fl¨u gelReviewed by:anonymous refereesReferencesArgent,R.:An Overview of Model Integration for Environmental Applications,Environmental Modelling and Software,19,219–234,2004.Bonham-Carter,G.F.:Geographic Information Systems for Geo-scientists,Elsevier Science Inc.,New York,1994. Burrough,P.A.:Dynamic Modelling and Geocomputation,in: Geocomputation:A Primer,edited by:Longley,P.A.,et al., Wiley,England,165-191,1998.Burrough,P.A.and McDonnell,R.:Principles of Geographic In-formation Systems,Oxford University Press,New York,1998. Gao,P.,Zhan,C.,and Menon,S.:An Overview of Cell-Based Mod-eling with GIS,in:Environmental Modeling with GIS,edited by: Goodchild,M.F.,et al.,Oxford University Press,325–331,1993.Goodchild,M.:A Geographer Looks at Spatial Information Theory, in:COSIT–Spatial Information Theory,edited by:Goos,G., Hertmanis,J.,and van Leeuwen,J.,LNCS2205,1–13,2001.Jørgensen,S.and Bendoricchio,G.:Fundamentals of Ecological Modelling,Elsevier,New York,2001.Julien,P.Y.,Saghaﬁan,B.,and Ogden,F.:Raster-Based Hydro-logic Modelling of Spatially-Varied Surface Runoff,Water Re-sources Bulletin,31(3),523–536,1995.Moore,I.D.,Turner,A.,Wilson,J.,Jenson,S.,and Band,L.:GIS and Land-Surface-Subsurface Process Modeling,in:Environ-mental Modeling with GIS,edited by:Goodchild,M.F.,et al., Oxford University Press,New York,1993.Press,W.,Flannery,B.,Teukolsky,S.,and Vetterling,W.:Numeri-cal Recipes in C:The Art of Scientiﬁc Computing,2nd Ed.Cam-bridge University Press,Cambridge,1992.Pullar,D.:MapScript:A Map Algebra Programming Language Incorporating Neighborhood Analysis,GeoInformatica,5(2), 145–163,2001.Pullar,D.:Simulation Modelling Applied To Runoff Modelling Us-ing MapScript,Transactions in GIS,7(2),267–283,2003. Ritter,G.,Wilson,J.,and Davidson,J.:Image Algebra:An Overview,Computer Vision,Graphics,and Image Processing, 4,297–331,1990.Sethian,J.A.:Level Set Methods and Fast Marching Methods, Cambridge University Press,Cambridge,1999.Sklar,F.H.and Costanza,R.:The Development of Dynamic Spa-tial Models for Landscape Ecology:A Review and Progress,in: Quantitative Methods in Ecology,Springer-Verlag,New York, 239–288,1991.Sui,D.and R.Maggio:Integrating GIS with Hydrological Mod-eling:Practices,Problems,and Prospects,Computers,Environ-ment and Urban Systems,23(1),33–51,1999.Tilman,D.and P.Kareiva:Spatial Ecology:The Role of Space in Population Dynamics and Interspeciﬁc Interactions.Princeton University Press,Princeton,New Jersey,USA,1997. Wesseling C.G.,Karssenberg, D.,Burrough,P. A.,and van Deursen,W.P.:Integrating Dynamic Environmental Models in GIS:The Development of a Dynamic Modelling Language, Transactions in GIS,1(1),40–48,1996.。

叶轮机械非定常流动及气动弹性计算

中图分类号：V211.3 论文编号：1028701 18-B061 学科分类号：080103博士学位论文叶轮机械非定常流动及气动弹性计算研究生姓名周迪学科、专业流体力学研究方向气动弹性力学指导教师陆志良教授南京航空航天大学研究生院航空宇航学院二О一八年十月Nanjing University of Aeronautics and AstronauticsThe Graduate SchoolCollege of Aerospace EngineeringNumerical investigations of unsteady aerodynamics and aeroelasticity ofturbomachinesA Thesis inFluid MechanicsbyZhou DiAdvised byProf. Lu ZhiliangSubmitted in Partial Fulfillmentof the Requirementsfor the Degree ofDoctor of PhilosophyOctober, 2018南京航空航天大学博士学位论文摘要气动弹性问题是影响叶轮机械特别是航空发动机性能和安全的一个重要因素。

作为一个交叉学科，叶轮机械气动弹性力学涉及与叶片变形和振动相关联的定常/非定常流动特性、颤振机理以及各种气弹现象的数学模型等的研究。

本文基于计算流体力学(CFD)技术自主建立了一个适用于叶轮机械定常/非定常流动、静气动弹性和颤振问题的综合计算分析平台，并针对多种气动弹性问题进行了数值模拟研究。

主要研究内容和学术贡献如下：由于叶轮机械气动弹性与内流空气动力特性密切相关，真实模拟其内部流场是研究的重点之一。

基于数值求解旋转坐标系下的雷诺平均N–S(RANS)方程，首先构造了适合于旋转机械流动的CFD模拟方法。

特别的，针对叶片振动引起的非定常流动问题，采用动网格方法进行模拟，通过一种高效的RBF–TFI方法实现网格动态变形；针对动静叶排干扰引起的非定常流动问题，采用一种叶片约化模拟方法，通过一种基于通量形式的交界面参数传递方法实现转静子通道之间流场信息的交换。

救援后送机器人对人体的作用力分析和新型执行机构结构与运动综合

救援后送机器人对人体的作用力分析和新型执行机构结构与运动综合The Rescue and Evacuation Robot：Human Body Stress Analysis and the Structure, Kinematics Methodology of a Novel HandlingMechanism一级学科：机械工程学科专业：机械制造及其自动化作者姓名：苏卫华指导教师：黄田教授天津大学机械工程学院二○一三年五月独创性声明本人声明所呈交的学位论文是本人在导师指导下进行的研究工作和取得的研究成果，除了文中特别加以标注和致谢之处外，论文中不包含其他人已经发表或撰写过的研究成果，也不包含为获得天津大学或其他教育机构的学位或证书而使用过的材料。

与我一同工作的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示了谢意。

学位论文作者签名：签字日期：年月日学位论文版权使用授权书本学位论文作者完全了解天津大学有关保留、使用学位论文的规定。

特授权天津大学可以将学位论文的全部或部分内容编入有关数据库进行检索，并采用影印、缩印或扫描等复制手段保存、汇编以供查阅和借阅。

同意学校向国家有关部门或机构送交论文的复印件和磁盘。

（保密的学位论文在解密后适用本授权说明）学位论文作者签名：导师签名：签字日期：年月日签字日期：年月日摘要本文以城市、丛林等地域军事行动和大型公共场所生化恐怖袭击、爆炸等突发公共卫生事件中对伤员实施救援所面临的卫勤保障难题为背景，研究救援后送机器人对人体的作用分析和新型执行机构结构与运动综合方法。

全文取得了如下创造性成果：(1) 以分形理论为工具，研究了利用CT扫描影像准确建立人体脊柱有限元模型的快速建模方法。

提出了一种考虑像素邻域特征的最佳阈值求解方法，该方法弥补了二维Otsu法在Mimics中的应用局限以及软件自身在分割阈值选取上的不足且建立的脊柱几何模型为灰度值模型，便于使用密度-CT数的弹性模量模型对骨组织材料进行非线性赋值，提高有限元模型的分析精度。

连续体机器人逆运动学英文

连续体机器人逆运动学英文In the realm of robotics, the inverse kinematics of a serial manipulator is a complex yet fascinating challenge. It involves calculating the joint angles that allow the end-effector to reach a desired position in space.The process begins with a clear understanding of the robot's geometry and the constraints of its joints. Eachjoint's range of motion and the lengths of the links must be meticulously defined to set the stage for the calculations.By applying the forward kinematics equations, we can express the end-effector's position as a function of thejoint variables. However, the inverse kinematics problem is non-trivial, as it requires solving these equations in reverse to find the joint angles for a given end-effector position.One common approach to solving the inverse kinematics problem is the use of Jacobian matrices, which relate the velocity of the end-effector to the joint velocities. This method is particularly useful for real-time applications where speed is crucial.Another technique involves iterative methods, such as the Newton-Raphson method, which incrementally adjusts the joint angles based on the error between the current and desiredend-effector positions.Despite the complexity, the successful implementation of inverse kinematics is essential for tasks that require precision, such as assembly lines or surgical procedures, where the robot must accurately position its tools.In conclusion, the inverse kinematics of serial manipulators is a cornerstone of modern robotics, enabling robots to perform intricate and precise movements that are vital for a wide range of applications.。

雅可比矩阵ppt课件

67
动力学普遍方程的补充：
A
问题的引出
M
m1g m2g
O
BF
m3g
MA
m1g m2g
O
B
m3g
问题1：系统在图示位置平衡，用什么方法求 F与M的关系？
38
机器人的奇异点讨论：
39
斯坦福机械手的运动学奇点：
40
斯坦福机械手的运动学奇点示例（讨论theta 5=0的特殊情况）
（theta 5=0时两轴线重合）
41
通过雅可比矩阵求解平面机械手的奇点分析示例：
42
43
通过雅可比矩阵对斯坦福机械手的奇点分析说明：
44
2.2 机器人静力分析
机器人在工作状态下会与环境之间引起相互作用的力和力矩。机器人各关节的驱动装置提供关节力和力矩，通过连杆传递到末端执行器，克服外界作用力和力矩。关节驱动力和力矩与末端执行器施加的力和力矩之间的关系是机器人操作臂力控制的基础。
velocity of 0P3org. (c) For what joint values is the manipulator at a singularity? What motion is restricted at this singularity?
66
2.3 机器人动力学方程机器人动力学的研究有牛顿-欧拉(Newton-Euler) 法、拉格朗日(Langrange)法、高斯(Gauss)法、凯恩(Kane)法及罗伯逊-魏登堡(Roberon-Wittenburg) 法等。本节介绍动力学研究常用的牛顿-欧拉方程和拉格朗日方程。
2 63
64
2. You are given that a certain RPR manipulator has the following transformation matrices, where {E} is the frame of the end ffector.

《Marc中常用词汇》

网格的划分(MESH GENERATION)与有限元分析相关的常用词：ELEMENT （单元）由多个节点定义的用于分析的最基本区域。

NODE （节点）用于定义单元的点，具体位置由坐标值确定。

与几何实体相关的常用词：POINT （点）描述曲线、曲面的控制点。

CURVE （曲线）线段、圆弧、样条等曲线的统称。

SURFACE （面）四边形面、球面、圆柱面等曲面的统称。

节点的生成：在 MESH ENERATION 菜单下方，有橙色的条目NODE，其右边依次为绿色的ADD、REM、EDIT、SHOW 光钮，分别表示生成、删除、修改、确认节点，选取NODE－ADD 后将<↑>移至格栅中心，按鼠标器左键，则在该点周围有一红色“”表示已将该点生成为节点，注意此时只有格栅点才能被检取生成为节点，同理依次将格栅点（1，0，0）（1，1，0）（0，1，0）生成为节点。

单元的生成当节点已经存在时，选取ELEMS－ADD，用鼠标器按逆时针的顺序检取节点，将<↑>移至节点附近，按<ML>，该节点变为黄色，按<MM>可以取消最近一次检取。

单元几何类型的定义在绿色 ELEMENT CLASS 光钮右边有绿色的“QUAD（4）”表示当前MENTAT 作成的单元几何形状类型为QUAD（4）（四节点四边形单元）。

如果要生成由二节点组成的直线形单元，则先检ELEMENT CLASS，进入下图所示的子菜单，检取LINE （2）并返回到MESH GENERATION 菜单，检取ELEMS－ADD，然后检取节点即可生成LINE（2）型单元。

MENTAT 可以支持以下线类型：LINE（直线）、CUBIC SPLINE（三次样条曲线）（important）、POLY LINE（多折线）、BEZIER（Bezier 曲线）、 NURB （Non Uniform Rational B-spline 非等分B 样条曲线）、INTERPOLATE（插值曲线）、COMPOSITE（复合曲线）、FILLET （倒圆线）、ARC（圆弧）、CIRCLE（圆）面的生成QUAD（四边形面，输入四个控制点的坐标值）、BEZIER（Bezier 曲面）、DRIVEN （驱动曲面，必须指定被驱动的曲线（DRIVEN）及驱动曲线（DRIVE））、NURB（NURB 曲面）、RULER（直纹曲面）、SPHERE（球面）、CYLINDER （圆柱面、圆锥面）、SWEPT（扫描面，输入二条曲线，即扫描线（SWEPT）及轨线（SWEEPING））、COONS（高斯面）、INTERPOLATE（插值面）、SKIN（蒙皮面）几何实体与网格的转换MESH GENERATION——CONVERT中，PONITS TO NODES 将控制几何点转换为单元节点。

Tripod机器人运动学模型仿真

Tripod机器人运动学模型仿真王强;石红瑞【摘要】以Tripod并联机器人为研究对象,通过对Tripod机器人结构的分析建立出运动位移方程,推导出Tripod机器人的逆运动学求解模型.对运动位移方程进行分析,推导出雅可比矩阵并得到速度及加速度正解模型.采用三维实体建模软件(UG)建立出Tripod机器人的实体装配模型,并将三维模型导入机械系统动力学自动分析软件(Adams)中进行运动学仿真,通过给定路径和主动臂旋转角速度两种仿真方案,验证了运动学模型的准确性.%Taking Tripod parallel robot as study object,the motion displacement equation is established by analyzing its structure,and inverse kinematics solution model is derived.Through analyzing motion displacement equation,the Jacobi matrix is deduced,and the forward solution model of velocity and acceleration are obtained.The real assembly model of Tripod robot is constructed with applying 3D entity modeling software UG.3D model is introduced into Adams in mechanic dynamics system for kinematic simulation.The accuracy of the kinematics model is verified by two simulation scenario of specified path and the active arm rotation speed.【期刊名称】《石油化工自动化》【年(卷),期】2018(054)001【总页数】6页(P26-30,68)【关键词】Tripod机器人;运动位移方程;逆运动学;雅可比矩阵;模型【作者】王强;石红瑞【作者单位】东华大学信息科学与技术学院,上海201620;东华大学信息科学与技术学院,上海201620【正文语种】中文【中图分类】TP242机器人的分类多种多样，按照结构可以将机器人分为：串联机器人、并联机器人和混联机器人。

Stewart平台雅可比矩阵分析

Stewart平台雅可比矩阵分析赵慧[1]张尚盈[2][1]武汉科技大学机械自动化学院 430081Email:[2]华中科技大学数字制造及设备技术国家重点实验室 430074Email:摘要：雅可比矩阵是对Stewart平台进行分析时的重要变量，通过对其的分析和计算，可以得到平台速度和液压缸速度之间的关系，得到平台承载与各液压缸出力之间的关系，可以判断液压缸的可控性，可以得到各自由度之间的运动耦合情况。

因此，导出雅可比矩阵，并对其物理意义进行诠释和深刻理解非常重要。

本文通过Stewart平台的运动学分析，推导出雅可比矩阵的公式，并通过仿真结果对其物理意义进行验证。

关键词：Stewart平台，运动学分析，雅可比矩阵1 引言随着科技的发展以及人们对未知世界探索的需求，Stewart平台在飞行模拟器、空中交会对接(RVD)仿真技术[1]、虚拟轴机床、力－扭矩传感器、装配机械手等领域有广泛的应用。

其中液压驱动Stewart平台由于具有快速、高精度、大负载和结构紧凑等特点而受到青睐 [2]。

Stewart平台是一个典型的多变量和本质非线性的复杂系统。

对Stewart平台运动学和动力学进行研究，是设计、分析和控制Stewart平台的基础。

雅可比矩阵是在对Stewart平台进行运动学动力学分析过程中产生和定义的矩阵，具有重要的物理意义，本文将对其实质展开论述，并用仿真结果来验证。

2 Stewart平台描述2.1 坐标系建立如图1所示，Stewart平台的主体部分由上平台（Platform）、下平台（Base）以及六个液压缸组成。

静止不动的下平台与可动作的上平台分别通过上、下胡克铰与液压缸的两端相连。

选取体坐标系{}P—O X Y Z在上平台上，坐p p p p标原点p O 为上铰点的外接圆圆心；惯性坐标系{}G —g g g g O X Y Z 的坐标原点g O 为下铰点的外接圆圆心；坐标轴的方向如图1所示。

数学英汉词典.docx

数学英汉词典-微积分MATHEMATICAL TERMS(Part 1) calculus 微积分definition 定义theorem 定理lemma 引理corollary推论prove 证明proof 证明show 证明solution 解formula 公式if and only if ( iff ) 当且仅当∀x∈X for all x∈X∃x∈X there exists an x∈X such that 使得given 已知set集合finite set有限集infinite set 无限集interval区间open interval开区间closed interval 闭区间neighborhood 邻域number 数natural number 自然数integer 整数odd number 奇数even number 偶数real number 实数rational number 有理数irrational number 无理数positive number 正数negative number 负数mapping 映射function 函数monotone function 单调函数increasing function 增函数decreasing function 减函数bounded function 有界函数odd function 奇函数even function 偶函数periodic function 周期函数composite function 复合函数inverse function 反函数domain 定义域range 值域variable 变量independent variable自变量dependent variable因变量sequence 数列convergent sequence收敛数列divergent sequence 发散数列bounded sequence 有界数列decreasing sequence 递减数列increasing sequence 递增数列limit极限one-sided limit 单侧极限left-hand limit 左极限right-hand limit 右极限The Squeeze Theorem 夹挤定理infinity 无穷大infinitesimal 无穷小equivalent infinitesimal 等价无穷小infinitesimal of higher order 高阶无穷小order of infinitesimal 无穷小的阶infinitesimals of the same order 同阶无穷小increment 增量continuous function 连续函数continuity 连续性f(x) is continuous at x 在x连续f(x) is discontinuous at x 在x间断discontinuity 间断点discontinuity of the first (second) kind 第一（二）类间断点removable discontinuity 可去间断点jump discontinuity 跳跃间断点infinite discontinuity 无穷间断点intermediate value 介值The Intermediate Value Theorem 介值定理zero point 零点The Zero Point Theorem 零点定理root 根equation 方程uniform continuity 一致连续derivative 导数rate of change 变化率velocity 速度instantaneous velocity 瞬时速度tangent (line) 切线normal (line) 法线slope 斜率left-hand derivative 左导数right-hand derivative 右导数f(x) is differentiable at x f(x) 在x处可导(可微)differentiation 求导The Chain Rule 链式法则differentiation formulas 求导公式implicit function 隐函数explicit function 显函数implicit differentiation 隐函数求导logarithm 对数Logarithmic differentiation 对数求导法parameter 参数parametric equation 参数方程parametric curve 参数曲线hyperbolic function 双曲函数hyperbolic sine 双曲正弦hyperbolic cosine 双曲余弦hyperbolic tangent 双曲正切hyperbolic cotangent 双曲余切differential 微分differential quotient 微商approximate value 近似值error 误差relative error 相对误差absolute error 绝对误差invariance of differential form微分形式不变性higher derivative 高阶导数first derivative 一阶导数second derivative 二阶导数third derivative 三阶导数nth derivative n阶导数twice differentiable 二阶可导acceleration 加速度mean value 中值The Mean Value Theorem 中值定理Rolle’s Theorem 罗尔定理Lagrange’s Mean Value Theorem 拉格朗日中值定理Cauchy’s Mean Value Theorem柯西中值定理equality 等式inequality 不等式indeterminate form 不定型(未定式)indeterminate form of type 00(∞∞)0(∞∞)型未定式L’Hospital’s Rule 洛必达法则Taylor’s formula 泰勒公式polynomial 多项式nth-degree polynomial n次多项式remaind 余项Lagrange’s form of remainder拉格朗日型余项Peano’s form of remainder皮亚诺型余项Maclaurin formula 麦克劳林公式Taylor polynomial 泰勒多项式Mauclaurin polynomial 麦克劳林多项式polynomial approximation 多项式逼近accuracy 精确度margin 边际marginal cost 边际成本marginal revenue 边际收益elasticity 弹性density 密度mass 质量extreme value 极值local maximum value 极大值local minimum value 极小值(absolute) maximum 最大值(absolute) minimum 最小值stationary point 驻点（稳定点）critical point 临界点The First (Second) Derivative Test （极值的）一（二）阶判别法convex 凸的convex curve 凸曲线concave 凹的concave curve 凹曲线convex function 凸函数point of inflection 拐点asymptote 渐近线horizontal asymptote 水平渐近线vertical asymptote 垂直渐近线slant asymptote 斜渐近线curve sketching 作图sketch a curve 作图curvature 曲率The bisection method 二分法The secant method 弦位法Newton’s method 牛顿（切线）法The tangent method 切线法differential calculus 微分学integral 积分integral calculus 积分学definite integral 定积分indefinite integral 不定积分partition 分割Riemann sum 黎曼和integral sign 积分符号integrand 被积函数upper (lower) limit of integration 积分上（下）限integration 积分（求积）integrable 可积的f(x) is integrable on [a, b] integrable function 可积函数integrability 可积性sufficient condition 充分条件necessary condition 必要条件piecewise continuous 分段连续property 性质The mean value theorem of integral 积分中值定理The fundamental theorem of calculus 微积分基本定理Newton-Leibniz formula primitive function (anti-derivative) 原函数（反导数）The substitution rule for integration 换元积分法The inverse of the chain rule反链式法（凑微分法）integration by parts 分部积分法rational function 有理函数fraction 分式irreducible fraction 最简分式partial fraction 部分分式partial fraction decomposition部分分式分解vector 矢量free vector 自由矢量zero vector 零矢量magnitude of a vector 矢量的模unit vector 单位矢量scalar product 数量积dot product 点积vector product 矢量积cross product 叉积a is perpendicular (orthogonal) tob a与b垂直coordinate 坐标coordinate system 坐标系coordinate axis 坐标轴x-axis x轴coordinate plane 坐标面direction angle 方向角direction cosine 方向余弦rectangular coordinate system 直角坐标系octant 卦限the first octant 第一卦限variable 变量function of two (three) variables二（三）元函数function of several variables 多元函数independent variable 自变量dependent variable 因变量domain 定义域range 值域set of points 点集neighborhood 邻域interior point 内点boundary point 边界点bound 边界open set 开集closed set 闭集connected set 连通集region 区域open region 开区域closed region 闭区域bounded region 有界区域unbounded region 无界区域cluster point 聚点double limit 二重极限iterated limit 累次极限continuity 连续性increment 增量total increment 全增量partial increment 偏增量partial derivative 偏导数partial derivative of f(x,y) with respect to x ( y ) f(x,y)关于x（y）的偏导数higher partial derivative 高阶偏导数mixed partial derivative 混合偏导数Laplace equation 拉普拉斯方程total differential 全微分differentiable 可微chain rule 链式法则implicit function 隐函数implicit differentiation 隐函数微分法Jacobian determinant 雅可比行列式curve 曲线space curve 空间曲线tangent vector 切矢tangent line 切线normal plane 法平面surface 曲面normal vector 法矢normal line 法线tangent plane 切平面sphere 球面cylinder 柱面cone 锥面directional derivative 方向导数gradient 梯度gradient vector 梯度矢量f del flevel curve 等值线level surface 等值面local extremum 极值local maximum 极大值local minimum 极小值extreme value 最值absolute maximum (minimum)最大（最小）值stationary point (critical point) 驻点（临界点）conditional extremum 条件极值Lagrange multiplier 拉格朗日乘数method of Lagrange multiplier拉格朗日乘数法objective function 目标函数constraint 约束条件method of least square 最小二乘法field 场scalar field 数量场vector field 矢量场gradient field 梯度场potential field 势场potential function 势函数conservative field 保守场gravitational field 引力场force field 力场velocity field 速度场multiple integral 重积分double integral 二重积分iterated integral 累次积分region 区域region of integration 积分区域type X (Y) region X（Y）型区域order of integration 积分秩序reverse the order of integration 交换积分秩序polar coordinates 极坐标double integrals in polar coordinates 极坐标下的二重积分volume 体积lamina 平面薄片mass 质量density 密度moment about x-axis 关于x轴的（静）力矩center of mass 重心moment of inertia 转动惯量surface 曲面area of a surface 曲面的面积triple integral 三重积分rectangle 矩形rectangular coordinates 直角坐标系cylinder 柱面cylindrical coordinates 柱面坐标系sphere 球面spherical coordinates 球面坐标系change of variables in multiple integrals 重积分的变量替换Jacobian determinant 雅可比行列式line integral 曲线积分line integral with respect to arc length 对弧长的曲线积分（第一型）line integral with respect to x ( y ) 对坐标x（y）曲线积分（第二型）line integral of a vector field向量场的曲线积分smooth curve 光滑曲线piecewise smooth curve 逐段光滑曲线oriented curve 有向曲线orientation of a curve 曲线的方向work 功the line integral is independent of path 曲线积分与路径无关connected region 连通区域simply-connected region 单连通区域closed curve 闭曲线Green’s theorem 格林定理（公式）positive orientation of a curve 曲线的正向Fundamental theorem for line integrals 曲线积分的基本定理surface integral 曲面积分surface integral of a scalar field数量场的曲面积分（第一型）surface integral of a vector field向量场的曲面积分（第二型）orientable surface 可定向曲面oriented surface 有向曲面Möbius strip 莫比乌斯带Klein bottle 克莱因瓶one-sided surface 单侧曲面two-sided surface 双侧曲面closed surface 闭曲面flux 流量、通量electric flux 电通量divergence 散度rotation (curl) 旋度Gauss’ theorem 高斯定理（公式）The divergence theorem 散度定理（公式）Stokes’ theorem 斯托克斯定理（公式）curl theorem 旋度定理（公式）circulation of v around L v沿L 的环流量Hamilton operator 哈密顿算子harmonic field 调和场。

雅可比矩阵在对机械手关节速度与力矩分析中的应用

雅可比矩阵在对机械手关节速度与力矩分析中的应用
陈益民（漳州职业技术学院机械与自动化系，福建漳州 363000）
摘要：提出二自由度机械手末端执行器与各关节之间的运动速度分析方程式与静力学分析方程式，通过实例推
导和分析，可以得到机械手运动的瞬时状况和静态下的力或力矩平衡问题，并可说明雅可比矩阵在运动学和力学两
责任编辑：潘伟彬
Research on use of Jacobian matrix on analysis of manipulators joint velocity and torque
CHEN Yi-min （Dept. of Mechanical and Automation, Zhangzhou Institute of Technology，
而t可节表示从尺m空间到jrn空间的一种线性映射建立了机械臂关节空间和操作空间作用力之间的关系这种力的对应关系和微分运动之间的对应关系是类似的
第10卷第4期 2008年12月
闽西职业技术学院学报 Journal of Minxi Vocational and Technical College
Vol.10 No.4 December 2008
118
关于如图 1 末端位置（x，y）和关节位移（θ1，θ2）位置运动学方程通过几何法容易求得为：
x=l1cosθ1+l2cos（θ1+θ2）（1）
y=l1sinθ1+l2sin（θ1+θ2）由此建立了二自由度机械手臂顶端即机械爪运
动位置与各关节角度。
通过微分原理可得：
式中 JT 称为机械手的力雅可比。它表示在静态
平衡状态下，操作力向关节力映射的线性关系。