下载文档原格式
两轮差速驱动移动机器人运动模型研究重庆大学硕士学位论文（学术学位）学生姓名：***指导教师：王牛博士专业：控制科学与工程学科门类：工学重庆大学自动化学院二O一三年四月Motion Modeling of Two-wheel DifferentialDrive Mobile RobotA Thesis Submitted to Chongqing Universityin Partial Fulfillment of the Requirement for theMaster’s Degree of EngineeringByMa QinyongSupervisor by Dr. Wang niuSpecialty: Control Science and EngineeringCollege of Automation ofChongqing University, Chongqing, ChinaApril 2013摘要两轮差速驱动移动机器人结构简单，控制方便，是应用最为广泛的一种移动机器人，其运动模型是移动机器人研究的重要内容。
目前，关于两轮差速驱动移动机器人的相关研究中常常忽略了包括移动机器人底层驱动电机动态性能在内的相关动力学因素，但在实际机器人系统中，由于机器人载重等变化会引起驱动系统负载发生变化，从而影响驱动系统的动态响应过程，导致机器人的运动状态发生改变。
而两轮差速驱动移动机器人是一个多输入多输出的控制系统，其运动学模型具有典型的非完整约束，各驱动回路往往采用内部带有非线性环节的双闭环控制系统，是一种具有多个非线性环节的非线性系统，必须采用非线性建模方法建立模型，因此建立考虑移动机器人动力学行为的两轮差速驱动移动机器人运动模型，对于机器人运动的精确描述和控制具有十分重要的理论和实际意义。
基于以上思考，本文采用“类等效”建模方法建立了两轮差速驱动移动机器人运动模型，主要研究工作如下：①提出了基于“类等效”建模方法的两轮差速驱动移动机器人运动模型。
Package‘statnet’October14,2022Version2019.6Date2019-06-13Title Software Tools for the Statistical Analysis of Network DataDepends R(>=3.5),tergm(>=3.6.1),ergm.count(>=3.3),sna(>=2.4),tsna(>=0.3)Imports ergm(>=3.10.4),network(>=1.15),networkDynamic(>=0.10.0),mon(>=4.2)Suggests ergm.rank(>=1.2.0),ergm.ego(>=0.5),EpiModel(>=1.7.2),degreenet(>=1.3-3),latentnet(>=2.9.0),networksis(>=2.1-3),relevent(>=1.0-4),ndtv(>=0.13.0)BugReports https:///statnet/statnet/issuesDescription Statnet is a collection of packages for statistical network analysis that are designed to work together because they share common data representations and'API'design.They provide an integrated set of tools for the representation,visualization,analysis,and simulation of many different forms of network data.This package is designed to make it easy to install and load thekey'statnet'packages in a single step.Learn more about'statnet'at<>.Tutorials for many packages can be foundat<https:///statnet/Workshops/wiki>.For an introduction to functions in this package, type help(package='statnet').License GPL-3+file LICENSEURL NeedsCompilation noAuthor Mark S.Handcock[aut],David R.Hunter[aut],Carter T.Butts[aut],Steven M.Goodreau[aut],Pavel N.Krivitsky[aut](<https:///0000-0002-9101-3362>),Skye Bender-deMoll[aut],Martina Morris[aut,cre]Maintainer Martina Morris<**************>1Repository CRANDate/Publication2019-06-1408:00:06UTCR topics documented:statnet-package (2)update_statnet (4)Index6statnet-package Easily Install and Load the statnet Packages for Statistical NetworkAnalysisDescriptionstatnet is a collection of software packages for statistical network analysis that are designed to work together,with a common data structure and API,to provide seamless access to a broad range of network analytic and graphical methodology.This package is designed to make it easy to install and load multiple statnet packages in a single step.statnet software implements recent advances in network modeling based on exponential-family random graph models(ERGM),as well as latent space models and more traditional descriptive net-work methods.This provides a comprehensive framework for cross-sectional and dynamic network analysis:tools for description,network visualization model estimation,model evaluation,model-based network simulation.The statistical estimation and simulation functions are based on a central Markov chain Monte Carlo(MCMC)algorithm that has been optimized for speed and robustness.The code is actively developed and maintained by the statnet development team.New function-ality is being added over time.Detailsstatnet packages are written in a combination of R and C It is usually used interactively from within the R graphical user interface via a command line.it can also be used in non-interactive (or“batch”)mode to allow longer or multiple tasks to be processed without user interaction.The suite of packages are available on the Comprehensive R Archive Network(CRAN)at https:// /and also on the statnet project website at / The suite of packages has the following components(those automatically downloaded with the statnet package are noted):For data handling:•network is a package to create,store,modify and plot the data in network objects.The network object class,defined in the network package,can represent a range of relational data types and it supports arbitrary vertex/edge/graph attributes.Data stored as network objects can then be analyzed using all of the component packages in the statnet suite.(automatically downloaded)•networkDynamic extends network with functionality to store information about about evo-lution of a network over time,defining a networkDynamic object class.(automatically down-loaded)For analyzing cross-sectional networks:•sna is a set of tools for traditional social network analysis.(automatically downloaded)•ergm is a collection of functions tofit,simulate from,plot and evaluate exponential random graph models.The main functions within the ergm package are ergm,a function tofit linear exponential random graph models in which the probability of a graph is dependent upon a vector of graph statistics specified by the user;simulate,a function to simulate random graphs using an ERGM;and gof,a function to evaluate the goodness offit of an ERGM to the data.ergm contains many other functions as well.(automatically downloaded)•ergm.count is an extension to ergm enabling it tofit models for networks whose relations are counts.(automatically downloaded)•ergm.ego is an extension to ergm enabling it tofit models for networks based on egocentri-cally sampled network data.(separate download required)•ergm.rank is an extension to ergm enabling it tofit models for networks whose relations are ranks.(separate download required)•latentnet is a package tofit and evaluate latent position and cluster models for statistical networks The probability of a tie is expressed as a function of distances between these nodes in a latent space as well as functions of observed dyadic level covariates.(separate download required)•degreenet is a package for the statistical modeling of degree distributions of networks.It includes power-law models such as the Yule and Waring,as well as a range of alternative models that have been proposed in the literature.(separate download required)For temporal(dynamic)network analysis:•tsna is a collection of extensions to sna that provide descriptive summary statistics for tem-poral networks.(automatically downloaded)•tergm is a collection of extentions to ergm enabling it tofit discrete time models for temporal (dynamic)networks.The main function in tergm is stergm(the“s”stands for separable), which allows the user to specify one ergm for tie formation,and another ergm for tie dissolu-tion.The models can befit to network panel data,or to a single cross-sectional network with ancillary data on tie duration.(automatically downloaded)•relevent is a package providing tools tofit relational event models.(separate download re-quired)Additional utilities:•erterms provides a template for users who want to implement their own new ERGM terms.(separate download required)•networksis is a package to simulate bipartite graphs withfixed marginals through sequential importance sampling.(separate download required)•EpiModel is a package for simulating epidemics(separate download required)statnet is a metapackage;its only purpose is to provide a convenient way for a user to load the main packages in the statnet suite.Those can,of course,also be installed individually.Each package in statnet has associated helpfiles and internal documentation,and additional the information can be found on the statnet project website(/).Tutorials, instructions on how to join the statnet help mailing list,references and links to further resources are provided there.For the reference paper(s)that provide information on the theory and methodol-ogy behind each specific package use the citation("packagename")function in R after loading statnet.We have invested much time and effort in creating the statnet suite of packages and supporting material so that others can use and build on these tools.We ask in return that you cite it when you use it.For publication of results obtained from statnet,the original authors are to be cited as described in citation("statnet").If you are only using specific package(s)from the suite,please cite the specific package(s)as described in the appropriate citation("packgename").Thank you! Author(s)Mark S.Handcock<******************.edu>,David R.Hunter<****************.edu>,Carter T.Butts<**************>,Steven M.Goodreau<***************>,Pavel N.Krivitsky<*************.au>,Skye Bender-deMoll<********************>,Samuel Jenness(for EpiModel)<**************************>,andMartina Morris<**************>Maintainer:Martina Morris<*************>update_statnet Update the Component Packages of the Statnet SuiteDescriptionA wrapper around update.packages to update the component packages of Statnet Suite to theirlatest versions.Usageupdate_statnet(...,ask=FALSE,checkBuilt=TRUE,addURLs=character())Argumentsask,checkBuiltArguments to update.packages documentation.The defaults are different fromthose of that function.addURLs Optional repository URLs in addition to CRAN,such as http://statnet./preview.Defaults to none....Additional arguments to be passed to update.packages.DetailsUpdates the list component packages of Statnet Suite,using setRepositories and update.packages.Since there are no good ways to update packages once they are loaded,this function should be called immediately after restarting R.Valueupdate_statnet returns NULL invisibly.See AlsosetRepositories,update.packages,install.packagesExamples##Not run:#Update from CRANstatnet::update_statnet()#Update from s preview release,taking packages from CRAN#as neededstatnet::update_statnet(addURLs="/preview")##End(Not run)Index∗utilitiesupdate_statnet,4ergm,3gof,3install.packages,5network,2networkDynamic,3setRepositories,5statnet(statnet-package),2statnet-package,2update.packages,4,5update_statnet,46。
No. 2CN 11-5904/U J Automotive Safety and Energy, 2010, Vol. 1 158—162电动轮汽车由于在驱动轮处采用电动轮技术而实现了多电机驱动，代替了传统电动汽车的中央驱动方式。
一般地，电动轮指电机到所驱动的车轮之间的所有部件，最简单的结构就是将电机与车轮组合成为一个整体。
电动轮驱动方式的优点在于，取消了传统汽车的传动轴和差速器等部件，使传动系统简化，不仅可以提高传动效率，而且有利于整车布置，提高车辆的通过性能，非常有利于低地板大客车和军用车辆的设计；由于减速装置布置在车轮附近，而且采用多个电动轮驱动，可以降低车辆对电气系统和机械传动零部件的要求，适合传递大转矩，非常适合于在大型矿用汽车上应用。
2002年，美国通用汽车提出了线控四轮驱动燃料电池概念车Autonomy，2005年推出后轮采用电动轮驱动的燃料电池电动车Sequel，2003年丰田汽车公司在东京国际车展上展示了四轮驱动燃料电池车Fine-S，2006年4月在美国纽约汽车展上又推出四个电动轮驱两后轮驱动的电动轮汽车的动力学建模与仿真分析陈 勇1，陆中奎2，周秋丽1（1.北京信息科技大学，北京 100192；2. 北京福田汽车股份有限公司，北京 102206）摘 要：为分析电动轮汽车的非悬挂质量增加对行驶平顺性、操纵稳定性的影响，建立了两后轮驱动的电动轮汽车整车的11自由度动力学模型。
在MATLAB/Simulink环境下，建立了整车仿真分析模型，采用模拟的路面谱作为路面输入,可实现不同车辆参数、不同控制策略和不同分析目标的仿真,也可分析车轮与路面之间的动载荷、悬架变形和车身姿态（俯仰、侧倾和横摆）的变化。
分析结论对电动轮汽车的开发、悬架的改进以及控制策略的确定具有参考意义。
关键词： 电动汽车；电动轮；控制策略；平顺性；操纵稳定性中图分类号： U469.72Dynamic modeling and simulation analysis of an electricvehicle with two rear hub-motorsCHEN Yong1, LU Zhongkui2, ZHOU Qiuli1(1. Beijing Information & Science Technology University, Beijing 100192, China;2. Beiqi Fonton Motor Co. Lts, Beijing 102206, China)Abstract: An 11 degree-of-freedom dynamic model was constructed for an electric vehicle driven with two rear hub-motors to analyze the infl uence on ride quality and the handling characteristics of unsprung mass increase. A full vehicle simulation model was developed using the MATLAB/Simulink with a simulated road model as input. The simulation model can realize the varies simulations with different vehicle parameters, control strategies and analyzing goals, while it can also determine the changes of dynamic load on tires, suspension defl ection and attitude (including pitch, roll and yaw). The above analyzed conclusions can enhance the development of electric vehicle driven by hub-motors, while they support the design of suspension and control strategies.Key words: electric vehicle; hub-motor; control strategy; ride quality; handling characteristics收稿日期：2010-01-22基金项目：辽宁省科学技术计划项目(2008220025)；辽宁省高等学校优秀人才支持计划项目(RC-05-12)作者简介：陈勇（1966—），男(汉族），辽宁，教授。
Control Parameter(Active for: Harmonic, Spectrum, Modal, Range, and Time History)Currently all of CAESAR II's dynamic analyses act only on linear systems, so any non-linearitiesmust be linearized prior to analysis. This means that one-directional restraints will not lift off andreseat, gaps will not open and close, and friction will not act as a constant effort force. Therefore,for dynamic analyses, all non-linear effects must be modeled as linear - for example, aone-directional restraint must be modeled as either seated (active) or lifted off (inactive), and agap must be either open (inactive) or closed (active). This process is automated when the staticload case is selected here - CAESAR II automatically activates the non- linear restraints in thesystem to correspond to their status in the selected load case (the user may think of this as beingthe loading condition - for example Operating - of the system at the time at which the dynamicload occurs). It must be noted that this automated linearization does not always provide anappropriate dynamic model, and it may be necessary to select other static load cases or even tomanually alter the restraint condition in order to simulate the correct dynamic response.A static load case must precede the dynamics job whenever:1) There are spring hangers to be designed in the job. The static runs must be made in order todetermine the spring rate to be used in the dynamic model.2) There are non-linear restraints, such as one-directional restraints, large-rotation rods,bi-linear restraints, gaps, etc. in the system. The static analysis must be made in order todetermine the active status of each of the restraints for linearization of the dynamic model.3) There are frictional restraints in the job, i.e. any restraints with a nonzero (mu) value.(Active for: Harmonic, Spectrum, Modal, Range, and Time History)All of CAESAR II's dynamic analyses are currently linear, so non-linear effects must be linearized. Modeling of friction in dynamic models presents a special case, since friction actually impacts the dynamic response in two ways - static friction (prior to breakaway) affects the stiffness of the system, by providing additional restraint, while kinetic friction (subsequent to breakaway) actually affects the damping component of dynamic response; due to mathematical constraints, damping is ignored for all analyses except time history (for which it is only considered on a system-wide basis). CAESAR II allows friction to be taken into account through the use of this Friction Stiffness Factor. CAESAR II approximates the restraining effect of friction on the pipe by including stiffnesses transverse to the direction of the restraint at which friction was specified. The stiffness of these "frictional" restraints is computed as:Kfriction = (F)*()*(Fact)Where:Kfriction = stiffness of frictional restraint inserted by CAESAR IIF = the force at the restraint taken from the static solution= mu, friction coefficient at restraint, as defined in the static model Fact = Friction Factor from the control spreadsheetThis factor should be adjusted as necessary in order to make the dynamic model simulate the system's actual dynamic response (note that use of this factor does not correspond to any actual dynamic parameter, but is actually a "tweak" factor to modify system stiffness). Entering a friction factor greater than zero causes these friction stiffnesses to be inserted into the dynamics job. Increasing this factor correspondingly increases the effect of the friction. Entering a frictionMAX. NO. OF EIGENV ALUES CALCULATED (0-NOT USED)(Active for: Spectrum, Modal, and Time History)The first stage of the Spectrum, Modal, and Time History analyses, is the use of the Eigensolver algorithm to extract the piping system's natural frequencies and mode shapes. For the Spectrum and Time History analyses, the response under loading is calculated for each of the modes, with the system response being the sum of the individual modal responses. Obviously, the more modes that are extracted, the more the sum of those modal responses resembles the actual system response. The problem is that this algorithm uses an iterative method for finding successive modes, so extraction of a large number of modes usually requires much more time than does a static solution of the same piping system. The object is to extract sufficient modes to get a suitable solution, without straining computational resources. CAESAR II permits the user to specify - either through a mode number cutoff or a frequency cutoff - the number of modal responses to be included in the system results. This parameter is used, in combination with theFrequency Cutoff described below, to limit the maximum number of modes of vibration to be extracted during the dynamic analysis. If this parameter is entered as 0, the number of modes extracted is limited only by the frequency cutoff (and potentially, the number of degrees-of-freedom in the system model).If the analyst is more interested in providing an accurate representation of the system displacements, it may only be necessary to request the extraction of a few modes, allowing a rapid calculation time. However, if an accurate estimate of the forces, stresses, etc. in the system is the objective, calculation time grows as it becomes necessary to extract far more modes. This is particularly true in the case when solving a fluid hammer problem in the presence of axial restraints; often modes with natural frequencies of up to 300 Hz can be large contributors to the solution.The usual procedure for determining how many modes are sufficient is to extract a certain number of modes and review the results; then to repeat the analysis while extracting 5 to 10 additional modes, and comparing the new results to the old. If there is a significant change between the results, a new analysis is made, again extracting 5 to 10 more modes above those that were extracted for the second analysis. This iterative process continues until the results taper off, becoming asymptotic. This procedure has two drawbacks, the first one obvious - the time involved in making the multiple analyses, as well as the time involved in extracting the potentially large number of modes. The second drawback, occurring with Spectrum analysis, is less obvious - a degree of conservatism is introduced when combining the contributions of the higher order modes. Possible spectral mode summation methods include SRSS, ABSOLUTE, and GROUP - all methods that combine modal results as same-sign (positive) values. In reality, theory states that the rigid modes actually act in phase with each other, and should therefore be combined algebraically, thus permitting the response of some rigid modes to cancel the effect of other rigid modes (this is actually what occurs in a time history analysis). Because of this conservatism, it is actually possible to get results which exceed twice the applied load, despite the fact that the Dynamic Load Factor (DLF) of an impulse load cannot be greater than 2.0.An alternative method of ensuring that sufficient modes are considered in the dynamic model is through the use of the Included Mass Data Report. This report (available from the Dynamic Output Screen) is compiled for all spectrum and time history shock cases, whether missing mass (see description in the section Include Missing Mass Components) is to be included or not. It displays the percent of system mass along each of the three global axes, as well as the percent of total force, which has been captured by the extracted modes. The percent of system mass active along each of the three global axes (X-, Y-, and Z-) is calculated by summing the modal mass (corresponding to the appropriate directional degree-of-freedom) attributed to the(Active for: Spectrum, Modal, and Time History)As noted above, CAESAR II permits the user to specify either a number of modes or a frequency cutoff for extracting modes to be considered in the dynamic analysis. Modal extraction ceases when the Eigensolver extracts either the number of modes requested, or extracts a mode with a frequency above that of the Frequency Cutoff, whichever comes first. One recommendation for selection of a frequency cutoff point is that the user extract modes up to, but not far beyond, arecognized "rigid" frequency, and then include the missing mass correction (discussed in thesection Include Missing Mass Components). Choosing a cutoff frequency to the left of theresponse spectrum's resonant peak will provide a non-conservative result, since resonantresponses may be missed. During spectrum analysis, using a cutoff frequency to the right of thepeak, but still in the resonant range, will yield either overly- or underly-conservative results,depending upon the method used to extract the ZPA from the response spectrum. (In the case oftime history analysis, selecting a cutoff frequency to the right of the peak, but still in the resonantrange, will probably yield non-conservative results, since the missing mass force is applied with adynamic load factor of 1.0). Extracting a large number of rigid modes for calculation of thedynamic response may be conservative in the case of Spectrum analysis, since all spectral modal combination methods (SRSS, GROUP, ABS, etc.) give conservative results versusthe algebraic combination method (always used during time history analysis), which gives a morerealistic representation of the net response of the rigid modes. When the analysis type is SPECTRUM, MODES, or TIMEHIST, either this parameter or the(Active for: Spectrum/GROUP and Time History) This parameter does double duty, depending upon the analysis type. For a Spectrum analysis typewith GROUP modal Combination Method (as defined by USNRC Regulatory Guide 1.92), thisparameter specifies the frequency spacing defining each modal group - i.e., the percent (ofthe base frequency) between the lowest and highest frequency of the group. Regulatory Guide1.92 specifies the group spacing criteria as 10% (entered here as 0.1), so it is unlikely that the userwould ever wish to change the Closely Spaced Mode Criteria from the CAESAR II default valueof 0.1. For a Time History analysis type, this parameter is used to enter the length of the time slice,in milliseconds, to be used by the program during its step-by-step integration of the equations ofmotion for each of the extracted modes (CAESAR II uses theunconditionally stable Wilson qintegration method, so any size time step will provide a solution, with a smaller step providing greater accuracy - and more strain on computational resources). The time step should besufficiently small that it can accurately map the force vs. time load profile (i.e., the time stepshould be smaller than typical force ramp times). Additionally, the time step must be smallenough that the contribution of the higher order modes is not filtered from the response. For thisreason, it is recommended that the time step should be selected such that Time Step (in seconds)times Maximum Modal Frequency (in Hz) be less than 0.1. For example, if the modal frequencycutoff is set to 50 Hz, the time step should be set to a maximum of milliseconds: 0.002 sec x 50(Active for: Spectrum and Time History) When repeating a dynamic analysis, this parameter may be set to "Yes", causing CAESARII toskip the eigensolution (reusing the results of the earlier analysis), and only perform the computations for displacements, reactions, forces, and stresses. Activating this option is only valid after an initial eigensolution has been performed and is still available. Additionally, the mass and stiffness parameters of the model must be unchanged or the previous eigensolution isSpatial or Modal Combination First(Active for: Spectrum)This directive tells CAESAR II whether to combine the Spatial components or the Modal components of the load case first. When performing a spectrum analysis, each of the modal responses must be summed. In addition, if multiple shocks have been applied to the structure in more than one direction, the results from different directions must be combined - for example, spatially combining the X-direction, Y-direction, and Z-direction results. The question arises as to whether the spatial summations should precede or follow the modal summations. A difference in the final results (of Spatial first vs. Modal first) arises whenever different methods are used for the spatial and modal combinations. The combination of Spatial components first implies that the shock loads are dependent, while the combination of Modal components first implies that the shock loads are independent. Dependent and Independent refer to the time relationship between the X, Y, and Z components of the earthquake. With a dependent shock case, the X, Y, and Z components of the earthquake have a direct relationship - a change in the shock along one direction produces a corresponding change in the other directions. For example, this would be the case when the earthquake acts along a specific direction having components in more than one axis - such as when a fault runs at a 30o angle between the X- and Z-axes. In this case, the Z-direction load would be a scaled (by a factor of tan 30o), but otherwise identical version of the X-direction load. In this case, spatial combinations should be made first.An Independent shock is one where the X, Y, and Z time histories produce related frequency spectra but have completely unrelated time histories. It is the Independent type of earthquake that is far more common, and thus in most cases the modal components should be combined first. For example, IEEE 344-1975 (IEEE Recommended Practices for Seismic Qualification of Class 1E Equipment for Nuclear Power Generating Stations) states: "Earthquakes produce random ground motions which are characterized by simultaneous but statistically INDEPENDENT horizontal and vertical components." This is usually less of an issue for force spectrum combinations, since normally there are no separate spatial components to combine - i.e., there are not X-, Y-, and Z-shocks acting simultaneously. However, in the event that there is more than one potential force load (such as when there is a bank of relief valves that can fire individually or in combination), the spatial combination method may be used to indicate the independence of the loadings. For example, if two relief valves may or may not fire simultaneously (i.e., they are independent), the two shocks should be defined as being in different directions (for example, X- and Y-), and the combination method selected should be "Modal before Spatial". If under certain circumstances, the two valves will definitely open simultaneously (i.e., the loadings are dependent), the combination method should be "Spatial before Modal". (Otherwise, the direction defined for a force spectrum loading has no particular meaning.)Nuclear Regulatory Guide 1.92 (published in February, 1976) describes the requirements for combining spatial components when performing seismic response spectra analysis for nuclear power plants. Note that since all Time History combinations are done algebraically (in-phase)Spatial Combination Method (SRSS/ABS)(Active for: Spectrum)This parameter is used to define the method for combining the spatial contributions of the shocks in a single spectrum load case. This option is only used for spectrum runs with more than a single excitation direction. Since directional forces are usually combined vectorially, this points to a Square Root of the Sum of the Squares (SRSS) combination method as being most appropriate. An Absolute method is provided for additional conservativism.Note that since all Time History combinations are done algebraically (in-phase) this parameter(Active for: Spectrum)During a spectrum analysis, responses are calculated for each of the individual modes; these individual responses are then combined to get the total system response. Considering that the response spectrum yields the maximum response at any time during the course of the applied load, and considering that each of the modes of vibration will probably have different frequencies, it is probable that the peak responses of all modes will not occur simultaneously. Therefore an appropriate means of summing the modal responses must be considered. Nuclear Regulatory Guide 1.92 (published in February, 1976) defines the requirements for combining modal responses when performing seismic response spectra analysis for nuclear power plants. The four options presented there are also available, along with one other, for modal combinations under non-nuclear seismic and force spectrum analyses. There are 5 available modal combination methods:1) Grouping Method: This method is defined in USNRC Regulatory Guide 1.92. The Grouping Method attempts to eliminate the drawbacks of the Absolute and SRSS methods (see below) by assuming that modes are completely correlated with any modes with similar (closely spaced) frequencies, and are completely uncorrelated with those modes with widely different frequencies. Effectively, this method dictates that the responses of any modes which have frequencies within 10% of each other first be added together absolutely, with the results of each of these groups then combined with the remaining individual modal results using the SRSS method. Note that the 10% figure controlling the definition of a group may be changed by using the Closely Spaced Mode Criteria/Time History Time Step (ms) parameter.2) Ten Percent Method: This method is defined in USNRC Regulatory Guide 1.92. The Ten Percent Method is similar to the Grouping Method in that it assumes that modes are completely correlated with any modes with similar (closely spaced) frequencies, and are completely uncorrelated with those modes with widely different frequencies. The difference between this one and the preceding method is that the Grouping Method assumes that modes are only correlated with those that fall within the group - i.e., are within a 10% band, while this method assumes that modes are correlated with those that fall within 10% of the subject mode - effectively creating a 20% band - 10% up and approximately 10% down. Note that the 10% figure controlling the definition of closely spaced frequencies may be changed by using the Closely Spaced Mode Criteria/Time History Time Step (ms) parameter.3) Double Sum Method (DSRSS): This method is also defined in USNRC Regulatory Guide 1.92. This combination method is the most technically correct for earthquake loads, in that an attempt is made to estimate the actual intermodal correlation coefficient based upon empirical data. Note that the load duration (td) and the damping ratio () may be specified by using the Load Duration (Time History or DSRSS method) (sec.) and Damping (Time History or DSRSS) (ratio of critical) parameters described in the corresponding sections found earlier in this chapter.4) Absolute: This method states that the total system response is equal to the sum of the absolute values of the individual modal responses. (This is effectively the same as using the DSRSS method with all correlation coefficients equal to 1.0, or the Grouping method, with all modes being closely spaced.) This method gives the most conservative result, since it assumes that the all maximum modal responses occur at exactly the same time during the course of the applied load. This is usually overly-conservative, since modes with different natural frequencies will probably experience their maximum DLF at different times during the load profile.5) Square Root of the Sum of the Squares (SRSS): This method states that the total system response is equal to the square root of the sum of the squares of the individual modal responses. (This is effectively the same as using the DSRSS method with all correlation coefficients equal to 0.0, or the Grouping method, with none of the modes being closely spaced.) This method is based upon the statistical assumption that all modal responses are completely independent, with the maxima following a relatively uniform distribution throughout the duration of the applied load. This is usually non-conservative, especially if there are any modes with very close frequencies, since those modes will probably experience their maximum DLF at approximately the same time during the load profile. Note that since all Time History combinations are done algebraicallyInclude Missing Mass Components (Y/N)(Active for: Spectrum and Time History)During spectrum (either seismic or force spectrum) or time history analyses, the response of a system under a dynamic load is determined by superposition of modal results. One of the advantages of this type of modal analysis is that usually only a limited number of modes are excited and need to be included in the analysis. The drawback to this method is that although displacements may be obtained with good accuracy using only a few of the lowest frequencymodes, the force, reaction, and stress results may require extraction of far more modes (possibly far into the rigid range) before acceptable accuracy is attained. CAESAR II provides a feature, called the "Missing Mass Correction", which helps solve these problems. This feature offers the ability to include a correction which represents the contribution of the higher order modes not explicitly extracted for the modal/dynamic response, thus providing greater accuracy without additional calculation time. When this option is activated (by entering Yes for this parameter), the program automatically calculates the net (in-phase) contribution of all non-extracted modes and combines it with the modal contributions - avoiding the long calculation time associated with the extraction of the high order modes and the possible excessive conservativism of the summation methods. When considering force spectrum analyses for loads such as fluid hammer, relief valve, slug flow, etc., the applied load is not a function of the system mass, so the correction is truly a "missing force" correction. Including the missing mass correction is especially effective when performing force spectrum calculations. This is because force spectrum loads due to causes such as fluid hammer and slug flow act axially along the pipe, and therefore excite most highly the extremely rigid axial extension modes. Therefore the rigid response represented by the missingMissing Mass Combination Method (SRSS/ABS)(Active for: Spectrum)This directive defines the method used to combine the missing mass/force correction components (see description in an earlier section, Include Missing Mass Components) with the modal (dynamic) results. Research suggests that the modal and rigid portions of the response are statistically independent, so the SRSS combination method (CAESAR II's default) is usually most accurate. The Absolute combination method provides a more conservative result, based upon the assumption that the modal maxima occur simultaneously with the maximum ground acceleration. Missing mass components are combined following the modal combination.Note that even though missing mass components may be included during Time History analyses, all Time History combinations are done algebraically (in-phase), so this parameter has no effect onMissing Mass Combination Method (SRSS/ABS)(Active for: Spectrum)This directive defines the method used to combine the missing mass/force correction components (see description in an earlier section, Include Missing Mass Components) with the modal (dynamic) results. Research suggests that the modal and rigid portions of the response are statistically independent, so the SRSS combination method (CAESAR II's default) is usually most accurate. The Absolute combination method provides a more conservative result, based upon the assumption that the modal maxima occur simultaneously with the maximum ground acceleration.Missing mass components are combined following the modal combination.Note that even though missing mass components may be included during Time History analyses, all Time History combinations are done algebraically (in-phase), so this parameter has no effect onDirectional Combination Method (SRSS/ABS)(Active for: Spectrum)This directive specifies the method used for combining shock components acting in the same direction. This directive is used most typically with Independent Support Motion load cases, where it defines the way in which responses from different support groups caused by excitation in the same direction are combined. Additionally, if there are multiple uniform shock spectra acting in the same direction (although this is unusual), this directive would govern their combination. In general, directional combinations should be made using the absolute method. (As noted in the earlier section, Include Pseudostatic (Anchor Movement) Components, this is the USNRC recommendation for directional combination of pseudostatic responses.) However, in the case of force spectrum loads, if several loads (for example, several relief valve loads) are all defined with the same "shock direction", using an SRSS combination method would be a way of modeling these as independent loads, while using the Absolute method would model these as dependent loads.Note that since all Time History combinations are done algebraically (in-phase) this parameterMASS MA TRIX FormulationSelect either "Consistent" or "Lumped".The consistent mass model is a fairly well documented concept in the world of dynamic analysis. Most texts on the subject (for example, "Structural Dynamics - Theory and Computation" by Mario Paz) describe how to build the mass matrix. A lumped mass model makes very coarse simplifications (which often results in correspondingly coarse results) - the benefit of the lumped mass model being that it does not require a lot of memory for the data storage. The consistent mass matrix requires a lot more data storage, but it takes into consideration the effects of bending and other rotational effects of the beam on its mass distribution, giving a more realistic result. Note:1) If mass is added at a DOF, CAESAR II assumes that it was a concentrated mass, and put it on the on-diagonal term, effectively treating it as a lumped mass.2) If the user zeroes the mass at a DOF, we assume that he/she wants to eliminate consideration of that DOF, so we go and zero out all elements on that row/column.(Active for: Spectrum, Modal, and Time History)In almost all cases, the eigensolver will detect modal frequencies from the lowest frequency to the highest. Sometimes, when there is some strong directional dependency in the system, the modes may converge in the wrong order. This could cause a problem if the eigensolver reaches the cutoff number of modes (i.e., 20), but has not yet found the 20 modes with the lowest frequency (it may have found modes 1 through 18, 20, and 21, and would have found number 19 next). CAESARII checks for this anomaly using the Sturm Sequence calculation. This procedure determines the number of modes that should have been found between the highest and lowest frequencies found, and compares that against the actual number of modes extracted. If those numbers are different, the user is given a warning. For example, if 22 natural frequencies are extracted for a particular system, and if the highest natural frequency is 33.5 Hz, the Sturm Sequence check makes sure that there are exactly 22 natural frequencies in the model between zero and 33.5+p Hz, where p is a numerical tolerance. The Sturm Sequence check would fail in the case where there are two identical frequencies at the last frequency extracted. The significance of this failure can only be estimated by the user. For example, consider a system with the following natural frequencies: 0.6637 1.2355 1.5988 4.5667 4.5667If the user asks for only the first four natural frequencies, a Sturm Sequence failure would occur because there are five frequencies, rather than four, which exist in the range between 0.0 and4.5667 + p (where p calculates to 0.0041). To correct this problem the user can either:a)Increase the frequency cutoff by the number of frequencies not found. (This number is reported by the Sturm Sequence Check.)b)Increase the cutoff frequency some small amount, if the frequency cutoff terminated the eigensolution. This will usually allow the lost modes to fall into the solution frequency range.c)Fix the subspace size at 10 and rerun the job. Increasing the number of approximation vectors improves the possibility that at least one of them will contain some component of the missing modes, allowing the vector to properly converge. The default here is "Yes", and should be left alone unless the user has some specific reason for deactivating the check.。
第27卷第5期 水下无人系统学报 Vol.27No.52019年10月JOURNAL OF UNMANNED UNDERSEA SYSTEMS Oct. 2019收稿日期: 2018-10-25; 修回日期: 2018-12-10.基金项目: 国家重点研发计划重点专项(2017YFC0305902); 青岛海洋科学与技术国家实验室“问海计划：项目(2017WHZZB0101); 天津市自然科学基金重点基金(18JCZDJC40100).作者简介: 孙秀军(1981-), 男, 教授, 主要研究方向为海洋移动观测平台技术.[引用格式] 孙秀军, 王雷, 桑宏强. Petrel-II 200水下滑翔机动力学建模及仿真[J]. 水下无人系统学报, 2019, 27(5): 480-487.Petrel-II 200水下滑翔机动力学建模及仿真孙秀军1,2,3, 王 雷1, 桑宏强4(1. 河北工业大学 机械工程学院, 天津, 300130; 2. 中国海洋大学 物理海洋教育部重点实验室, 山东 青岛, 266100; 3. 青岛海洋科学与技术试点国家实验室 海洋动力过程与气候功能实验室, 山东 青岛, 266237; 4. 天津工业大学 机械工程学院, 天津, 300387)摘 要: 国内对水下滑翔机动力学行为的研究多针对深海型及横滚转向机制, 而对浅海型、尾舵转向类型研究较少。
基于此, 文中以Petrel-II 200浅海型水下滑翔机作为模型进行动力学建模及运动仿真分析, 并引入海流等干扰因素, 为浅海型水下滑翔机的运动形式提供参考。
根据Petrel-II 200三维模型中各质量的相对运动关系将其质心简化为由多个质点组成的多刚体系统, 构造质心与质点之间的关系式; 基于动量及动量矩定理对Petrel-II 200进行动力学分析, 对水下滑翔机所受重力、驱动浮力以及水动力进行体坐标系转换, 并推算出完整的动力学方程, 明确了升阻比与回转半径表达式; 通过选取水下滑翔机物理、水动力学参数, 对锯齿、螺旋等典型运动进行仿真试验。
Dynamic network modeling of two-phase drainage in porous mediaMohammed S.Al-Gharbi and Martin J.BluntDepartment of Earth Science and Engineering,Imperial College,London SW72AZ,United Kingdom͑Received21January2004;published13January2005͒We present a dynamic network model for modeling two-phaseflow.We account for wetting layerflow,meniscus oscillation,and the dynamics of snapoff.Interfaces are tracked through pore elements using amodified Poiseuille equation for the equivalent hydraulic resistance of thefluids between the pore elementcenters.The model is used to investigate the effects of capillary number and viscosity ratio on displacementpatterns and fractionalflow in primary drainage.We show that the amount of snapoff increases with increasingcapillary number and decreasing wetting phase viscosity.For capillary numbers lower than approximately10−5,the pore-scalefluid distribution and fractionalflow are similar to those obtained using a quasistatic modelthat ignores viscous forces.The contribution of oil transport from ganglia,formed by snapoff,is negligibleexcept for very large capillary numbers,greater than around0.1.DOI:10.1103/PhysRevE.71.016308PACS number͑s͒:47.55.Mh,47.55.Dz,47.55.KfI.INTRODUCTIONUnderstanding multiphaseflow in porous media is impor-tant for a number of practical applications,including theimplementation of improved recovery schemes from hydro-carbon reservoirs,contaminant cleanup,and designing un-derground nuclear waste repositories.Fundamentally,the ge-ometry of the void space of a porous medium and theinteractions of the multiple phases with the solid determinemacroscopic properties such as porosity,relative permeabil-ity,capillary pressure,and resistivity index͓1,2͔.One ap-proach to predicting these quantities is through pore-scalemodeling that requires a detailed understanding of the physi-cal processes occurring at the pore scale and a complete de-scription of the morphology of the pore space.The pore-scale network is a representation of the voidspace of the reservoir rock.Wide void spaces are representedby pores that are connected by narrower regions calledthroats.Then using rules for determiningfluid configurationsand appropriateflow and transport equations,multiphaseflow is computed in the network.An excellent description ofdifferent types of pore-scale network models can be found inthe classic text by Dullien͓2͔,while more recent reviews areto be found in Refs.͓3,4͔.The majority of the existing pore network models are qua-sistatic͓1,4–11͔assuming that capillary forces alone controlthefluid configuration in the pore space:capillary pressure isimposed on the network and thefinal,static position of allfluid-fluid interfaces is determined,ignoring the dynamic as-pects of pressure propagation and interface dynamics due toviscous forces.There are many important circumstances where the as-sumption of purely capillary-controlled displacement at thepore scale is not appropriate.The ratio of viscous to capillarypressure is conventionally defined using a capillary numberCa͓2͔,Ca=q,͑1͒whereis the viscosity of the displacing phase,q is the total Darcy velocity͑volumeflowing per unit area per unit time͒,andis the interfacial tension between the twofluid phases. For typical oil/waterflows in reservoirs,is around 10−3Pa s,is0.05N m−1,and q is10−5ms−1or lower, giving CaӍ2ϫ10−7.Viscous forces become significant at the pore scale only when Ca is the range10−4–10−3or larger ͓2͔.Hence quasistatic models are often accurate and can pre-dict a variety of low-flow-rate experiments successfully͓1,4͔.However,if the viscosity is high͑polymerflooding͒, theflow rate is very large͑for instance,in fractures,gas reservoirs,and near well bores͒or the interfacial tension is low͑surfactantflooding,near-miscible gas injection,and gas condensate reservoirs͒,viscous and capillary forces can be comparable at the pore scale.In these cases,a number of effects,such as the movement of disconnected ganglia of oil and the simultaneousfilling of neighboring pores become significant͓12͔.Dynamic network models explicitly account for viscous forces:a specified inflow rate for one of thefluids is imposed and the subsequent transient pressure response and the asso-ciated interface positions are calculated.Koplik and Lasseter ͓13͔simulated primary drainage in networks of spherical pores connected to cylindrical pore throats.Toubou et al.͓14͔and Blunt et al.͓15͔used a simplification of the model of Koplik,assuming that the pores have volume but no re-sistance toflow and the throats have resistance toflow but no volume.Payatakes and co-workers͓16–20͔simulatedflow in a network of spherical chambers connected through long cy-lindrical throats with a sinusoidally varying cross section.We will use a similar geometry in our work.They concluded that a significant fraction of the oilflow even at reservoir rates is accommodated through the movement of disconnected gan-glia.Celia and co-workers͓21–24͔extended the model of Blunt et al.͓15͔to study the effect of material heterogene-ities on the capillary pressure-saturation relationship͓21͔, effect of nonzero stress at thefluid-fluid interface͓22͔,inter-facial area and its relation to capillary pressure͓23͔,and interfacial velocity͓24͔.Recently Hansen and co-workers ͓25–27͔developed a dynamic network to studyfluid move-ment in drainage and imbibition.They used a network of tubes that connected to each other through volumeless nodes.PHYSICAL REVIEW E71,016308͑2005͒The capillary pressure within the tube was a function of the location of the interface.All these previously developed models except͓19͔only allowed a single phase to be present in any cross section through a throat—that is,they ignored contributions toflow through wetting layers that occupy the roughness and crev-ices of the pore space even when the center of the pore or throat isfilled by a nonwetting phase͓4͔.Blunt Scher͓28͔and Hughes and Blunt͓29͔did accommodate wetting layer flow in a perturbative model where the wetting layers were all assigned afixed conductance and where simultaneousfill-ing of multiple pores was disallowed.They showed that rate effects can have a significant effect on the trapping of the nonwetting phase͑oil͒during imbibition for capillary num-bers as low as10−6,since the ratio of the viscous to capillary force can be high,even at lowflow rates,ifflow occurs through low conductivity layers.Mogensen and Stenby͓30͔also assumed afixed conductance to the wetting layers in their dynamic model of imbibition.They concluded that the capillary number,aspect ratio͑ratio of average pore diameter to average throat diameter͒,and contact angle all have a significant influence on the competition between pistonlike advance͑frontal movement of a phase displacing another in a throat͒and snapoff͑where wetting phaseflows through layers andfills narrow regions of the pore space in advance of a connected front͒.Singh and Mohanty͓31͔developed a dynamic model to simulate two-phaseflow.They used a cu-bic network with cubic pores and throats of square cross section.A pseudopercolation model was included in the model for simulating low capillary numberflow.In their model,the volumeflowing in the wetting layers was set to be 1%of the volumeflowing through the bulk.The model was used to study primary drainage with constant inletflow rate. Saturation and relative permeability were computed as a function of capillary number,viscosity ratio,and pore-throat size distribution.Despite this extensive literature on dynamic pore-scale modeling,a number of key physical effects have yet to be captured accurately.In particular,previous work did not ac-count for the swelling of wetting layers in both drainage and imbibition that allows snapoff,as observed in micromodelexperiments͓14͔.Snapoff is the key process by which thenonwetting phase becomes trapped,and determines,for in-stance,the amount of oil that is left unrecovered after waterflooding.Related to this lack of physical realism is a contro-versy in the literature over the generic nature of multiphaseflow in porous media.The conventional picture,basedlargely on quasistatic approaches to modeling,assumes thatfor typical reservoir displacements each phase occupies itsown connected subnetwork through the porous medium.Thehydraulic conductance of these subnetworks determine themultiphaseflow properties,specifically the relative perme-ability͓2,4͔.Disconnected regions do notflow unless thecapillary number is very high.In contrast,Payatakes et al.͓16–20͔suggest based on micromodel experiments and net-work modeling that the typical scenario for multiphaseflowis significant transport via disconnected ganglia even at res-ervoirflow rates.However,some of their work can be criti-cized for not accommodating wetting layerflow and conse-quently substantially restricting the connectivity of thewetting phase.In this work we introduce a conceptually simple dynamic model that explicitly simulates the dynamics of wetting layer swelling and snapoff.We are then able to address whether or not multiphaseflow involves significant transport of the dis-connected nonwetting phase,even at typical reservoirflow rates,or whether this phenomenon is restricted to high cap-illary numbers.II.PRINCIPLES OF THE MODELThe model is based on three main principles.͑1͒The amount of each phase in each pore or throat is known at each time step.The volume of each phase͑with the contact angle͒controls the configuration offluids.This in turn determines the curvature of the oil-water interface and the pressure difference between these twofluids in each pore or throat.͑2͒By using equivalent networks of electrical resistors, the hydraulic resistances of thefluids between pore and throat centers are calculated and used in a volume balance equation to obtain thefluid pressures at pore and throat ing an equivalent resistor network simplifies the problem and makes it no more complex than solving the material balance equations for single-phaseflow.͑3͒The pressure difference between the pore and throat centers and the previously computed hydraulic resistances of each phase are used to move phases between pores and throats and hence to update thefluid volumes.The simula-tion then returns to step͑1͒.WORK AND PORE GEOMETRY The porous medium is represented by a square lattice of pores and throats.In cross section each pore or throat is a scalene triangle.The inscribed radius of a pore or throat varies sinusoidally,as shown in Fig.1͓9͔.Each pore is di-vided into several branches͑equal to the number of con-nected throats͒,which are considered as extensions of the throats they are connected to.All the pore branches meet at the center of the pore that is treated as a volumeless joining point.This feature is introduced for two reasons.First,we believe it is more realistic than assuming a straight channel ͑i.e.,uniform inscribed radius along the length͒.Second,it allows us to assign a unique and continuously varying cap-illary pressure as the interface͑meniscus͒moves in a pore or throat.The inscribed radius͑R͒at any point between the pore and the connecting throat centers is given byR=ͩR p+R t2ͪ+ͩR p−R t2ͪcosͩ2x l p+l tͪ,͑2͒where R p and R t are the pore and throat center radii,respec-tively,l p and l t are the pore and throat lengths,respectively. x=0at the pore center and x=͑l p+l t͒/2at the throat center.A.Selection of pore and throat sizesThe inscribed radius of the center of any throat is assigned at random according to a truncated Weibull distributionM.S.AL-GHARBI AND M.J.BLUNT PHYSICAL REVIEW E71,016308͑2005͒R t =͑R t ,max −R t ,min ͕͒−␦ln ͓z ͑1−e −1/␦͒+e −1/␦͔͖1/ +R t ,min ,͑3͒where R t is a radius of a throat center and z is a random number between 0and 1.The parameters used in the distri-bution are shown in Table I.The pore radius at its center must be greater or equal to the maximum radius of the con-necting throats.Therefore the pore radius is given by the following expression:R p =max ͕R ti ͉i =1,...,n ͖a ,͑4͒where n is the number of the connecting throats,and the aspect ratio ͑a ͒is the ratio between the pore radius and the maximum radius of the connecting throats.The aspect ratio is obtained by using Eq.͑3͒,replacing R t ,max and R t ,min by the maximum and minimum aspect ratios provided in Table I.While we use a topologically square lattice,we do allow the pore and throat lengths to vary—this corresponds physically to a distorted lattice,although we do not check if the network is physically realizable in two-dimensional space.We used the same parameters for selecting the pore and throat lengths,replacing R t ,max and R t ,min by l max and l min provided in Table I in Eq.͑3͒.Once l p ,l t ,R p ,and R t are known,Eq.͑2͒is used to determine the inscribed radius as a function of distance between the pore and throat centers.B.Determination of the pore cross sectionand corner half anglesEach pore and throat has a scalene triangular cross section with corner angles selected at random.The triangular cross section of an element is determined through two parameters:the inscribed radius ͑described above ͒,and the shape factor F =A /P 2,where A is the cross sectional area and P is the perimeter ͓32͔.The shape factor is used to determine the corner half angles for a triangle.In our model,we assume that the shape factor for each pore and throat is chosen ac-cording to the truncated Weibull distribution ͓Eq.͑3͔͒for a range of shape factor between zero and ͱ3/36͑equilateral ͒,F =A P 2=14͚i =13cot ͑␣i ͒=14tan ␣1tan ␣2cot ͑␣1+␣2͒,͑5͒where ␣1and ␣2are the two corner half angles subtended at the longest sides of the triangle.It is clear from the above equation that for a single value of the shape factor,there are ranges of corner half angles and triangular shapes.We follow the procedure of Patzek ͓10͔to select a nonunique solution for corner half angles.͑1͒Select the upper and lower limits of the second largest corner half-angle.These two limits are given according to the following equations␣2,min =arctanͫ2ͱ3cos ͩarccos ͑−12ͱ3F ͒3+43ͪͬ,͑6͒␣2,max =arctanͫ2ͱ3cos ͩarccos ͑−12ͱ3F ͒3ͪͬ.͑7͒͑2͒Pick randomly a value of ␣2between the two limits,␣2=␣2,min +͑␣2,max −␣2,min ͒ϫz ,where z is a random number between 0and 1.͑3͒The corresponding value of the largest corner half angle ␣1can be found from␣1=−12␣2+12arcsinͩtan ␣2+4F tan ␣2−4Fsin ␣2ͪ.͑8͒TABLE I.Parameters used to determine pore and throat geometries.Parameter Value R t ,min 0.2m R t ,max 100m l min 1m l max 50m a max 2.2a min2.0␦0.81.6FIG.1.Schematic diagram of pore and throat geometries.͑a ͒Side view:the inscribed radius of the pore element varies sinusoi-dally with the length.͑b ͒Cross-sectional view:the pore element has a triangular cross section.␣is the corner half angle.DYNAMIC NETWORK MODELING OF TWO-PHASE …PHYSICAL REVIEW E 71,016308͑2005͒͑4͒The smallest corner half angle is then obtained from ␣3=/2−͑␣1+␣2͒.C.Restricting snapoff to throatsIn theory it is possible for snapoff to occur in pores.Con-sider a pore that is connected to two throats.One branch has a wide pore/throat boundary radius and the other has a small pore/throat boundary radius.It is possible that a meniscus passing the wide boundary will cause snapoff at the narrow boundary,as we will describe later in the context of snapoff in throats.This is an unphysical effect.Therefore we limit snapoff to the throats.This can be done by ensuring that the radius at the pore/throat boundary is greater than half the radius at the pore center,Eq.͑9͒,R boundary Ͼ12R P .͑9͒Substituting Eq.͑2͒in Eq.͑9͒,ͫͩR p +R t 2ͪ+ͩR p −R t 2ͪcos ͩl pl p +l tͪͬϾ12R P .͑10͒Rearranging,1arccos ͩR t R t −R p ͪϾl pl p +l t.͑11͒If “=”sign is used in Eq.͑11͒instead of “Ͼ,”the final pore length can be given by the following expression:l p −final =͑l p +l t ͒ͫ1arccos ͩR t R t −R pͪͬ,͑12͒where ͑l p +l t ͒represents the spacing between the pore and the throat.If Eq.͑11͒is not satisfied,the pore length is re-duced using Eq.͑12͒with ͑l p +l t ͒held constant.IV .FLUID FLOW THROUGH THE NETWORKInitially,the system is assumed to be completely filled with a defending,wetting fluid ͑water ͒with viscosity 2.The invading nonwetting fluid ͑oil ͒with viscosity 1is in-jected into the system from the inlet side with a constant injection rate.The fluids are assumed to be immiscible and incompressible.We assume we have a water-wet system with an oil/water contact angle =0.A.Determination of fluid configurationsWe start by assuming the volume of each phase in each pore or throat is known.Then,with known contact angles,the fluid configuration and the local capillary pressure can be computed.The total cross-sectional area ͑A t ͒for an element ͑pore or throat ͒filled with a single phase isA t =R 2͚i =1ncot ͑␣i ͒=CR 2,͑13͒where n represents the number of the corners,␣is the cornerhalf angle,and C is a constant defined by Eq.͑13͒.There-fore,the volume for a fluid that occupies the whole cross section can be given by the following equation:V =C͵x 1x 2R 2͑x ͒dx ,͑14͒where R ͑x ͒is given by Eq.͑2͒.Here,x 1and x 2are the limits of the integration that are determined by the location of the fluid interfaces.For example in Fig.2͑a ͒,the water fills the pore from the center up to the first oil-water interface,so here,x 1=0.0and x 2=the location of the interface.In the case of water in wetting layers and the oil in the center,the wetting layer volume is equivalent to the volume of a cylinder that has cross-sectional area equal to the wet-ting layer cross-section and length equal to the difference of the location of the fluidinterfaces,FIG.2.Illustration of a fluid configuration:oil occupies thepore/throat boundary.͑a ͒Side view:the inscribed radius of the pore element ͓R ͑x ͔͒varies sinusoidally with the length.͑b ͒Cross-sectional view at the pore/throat boundary:the pore element has a triangular cross section.Oil is in the center and water is in the corners ͑wetting layers ͒.is the contact angle and r is the radius of curvature of the wetting layers.We assume that r is same in all the three corners and constant in each pore or throat,but varies between pores and throats and over time.M.S.AL-GHARBI AND M.J.BLUNT PHYSICAL REVIEW E 71,016308͑2005͒A ci =r 2ͩcos ͑cot ␣i cos −sin ͒++␣i −2ͪ,͑15͒where A ci is the cross sectional area of the i th wetting layer,r is the radius of the curvature of the wetting layer,and is the contact angle.The total cross sectional area for all the wetting layers isA c =͚i =1nA ci ,͑16͒where n represents the number of the corners.The oil volume will be the difference between the total volume given by Eq.͑14͒and the wetting layer volume.For example,in Fig.2͑a ͒,the wetting layer volume between the interface ͑x =x 2͒and the pore-throat boundary is obtained by computing the wetting layer cross-sectional area at the loca-tion of the interface ͑x =x 2͒using Eq.͑16͒and multiplying it by the difference between the location of the interface and the pore-throat boundary.Equation ͑14͒is used to find the volume between the interface ͑x =x 2͒and the pore-throat boundary.The oil volume then will be the difference be-tween these two volumes.For simplicity,in the volume cal-culations we assume that fluid interfaces,except in the wet-ting layers,are flat.However,interfacial curvature is accounted for in the computation of capillary pressure.puting the fluid resistanceCalculating the fluid resistance is potentially a compli-cated problem,but using an equivalent electrical resistors network helps to simplify the computations.Before giving a detailed description of how the equivalent hydraulic resis-tance is obtained,it is worthwhile to state the expressions used in finding the resistance of each phase.A phase may occupy:the whole element cross section,the corners with the other ͑nonwetting ͒phase in the center,or the center of the element with the other ͑wetting ͒phase in the corners.The general form of the fluid hydraulic resistance ͑W ͒for the three regions isW =f͵x 1x 2dxG ͑x ͒region,͑17͒where the subscript region stands for the three regions:whole cross section,center,or wetting layers;x 1and x 2are the location of the interfaces;G ͑x ͒is the fluid conductance per unit length,and f is the fluid viscosity.The conductance per unit length of a fluid occupying the whole cross section is given by the following approximation based on Poiseuille’s law for flow in a circular cylinder ͓7͔:G =128ͩͱA t+R ͪ4,͑18͒where A t is the cross-sectional area given by Eq.͑13͒and R is the inscribed radius given by Eq.͑2͒.Since R in Eq.͑18͒is function of x ,the integration in Eq.͑17͒is evaluated nu-merically.In the case where fluid occupies the corners with +␣Ͻ/2,the conductance per unit length is given by the fol-lowing approximation ͓33͔:G =͚i =1nͩA ci 2͑1−sin ␣i ͒2͑2cos −1͒3212͑sin ␣i ͒2͑1−3͒2͑2+f 1͒2ͪ,͑19͒where 1=/2−␣i ,2=cot ␣i cos −sin ,3=͑/2−␣i ͒tan ␣i ,and A ci is the corner area given by Eq.͑15͒.f is used to indicate the boundary condition at the fluid/fluid in-terface,f =1represents a no-flow boundary,while f =0is a free boundary.In Eq.͑19͒we assume f =1.The curvature of the wetting layer is assumed to be constant in a single ele-ment ͑i.e.,G is not function of x ͒,although it varies between pores and throats.Thus we can write the hydraulic wetting layer flow resistance ͑W l ͒asW l =f ͑x 2−x 1͒G,͑20͒ignoring the curvature of the pores and throats in the x di-rection.When fluid occupies the center of an element with wetting phase in the corners,the conductance per unit length is given by Eq.͑18͒replacing A t by A cen ,where A cen =A t −͚i =1n A ci ,A ci is the corner area,n represents the number of ing these formulas enables the model to handle any number of fluid interfaces between pore and throat centers.The use of the electrical resistors diagram and equivalent hydraulic resistance simplifies and clarifies this approach.For example,Fig.2͑a ͒shows a fluid configuration where oil occupies the center at the pore/throat boundary.Its equiva-lent electrical resistors diagram is given in Fig.3.Thus the equivalent hydraulic resistance of the fluids between the in-terfaces ͑W eq interface ͒is given by the following equation:1W eqinterface=1W P -l w +W T -lw+1W P -c o +W T -co,͑21͒where W P -l w is the pore wetting layer resistance,W T -l wis the throat wetting layer resistance,W P -c o is the pore oil resistance,and W T -c ois throat oil resistance.Then,the equivalent hydrau-lic resistance for the fluids between pore and the throat cen-ters ͑W eq ͒can be obtained through the following expression:W eq =W P -c w +W eqinterface +W T -c w,͑22͒where W P -c w is the pore water resistance and W T -c wis the throat water resistance.This is only one example of how the equivalent hydraulic resistance iscomputed.FIG.3.The equivalent electrical resistors diagram for Fig.2͑a ͒used to compute hydraulic resistance ͑W ͒.DYNAMIC NETWORK MODELING OF TWO-PHASE …PHYSICAL REVIEW E 71,016308͑2005͒The Appendix lists the expressions used for the different fluid configurations considered in the model.We only allow up to two menisci to be present in each throat and one me-niscus in each pore branch.V .SOLVING FOR THE FLUID PRESSURESince the equivalent hydraulic resistance is used in thevolume conservation equations,the model can be considered as solving a single-phase flow problem in which the phase conductance between the pore and throat centers is found from the equivalent hydraulic resistance calculated in the previous step.Therefore,the total flow rate between the pore center and throat center will beQ total =P p −P t +P cW eq,͑23͒where P p and P t are the pore center pressure and throat cen-ter pressure,respectively,P c is the sum of the capillary pres-sures of the menisci between the pore and throat centers.For instance,in Fig.2͑a ͒,P c =P c 1+P c 2,where,P c 1is the capil-lary pressure at the first interface ͑x =x 2͒and P c 2is the cap-illary pressure at the second interface ͑x =x 3͒.The absolute value of the capillary pressure at any meniscus is given by the following expression:P c =2cos ͑+␥͒R ͑x ͒,͑24͒where ␥is the inclination angle;tan ␥=dR ͑x ͒/dx .The sign of the capillary pressure at any meniscus depends on the location of the nonwetting fluid ͑oil ͒.If it is on the right of the meniscus,the sign is positive,otherwise it will be nega-tive.For a system of m throats and n pores,we have m +n unknowns.These unknowns are determined by applying vol-ume conservation for both the pores and throats.The conser-vation equation for pore j with n connecting throats,labeled,i is͚i =1nQ total ij=͚i =1nP p j−P t i +P cijW eqij=0.͑25͒For a throat i ,where R and L label the left and right pores,P p R −P t i +P cRiW eqRi+P p L −P t i +P cLiW eqLi=0.͑26͒We use Eq.͑26͒to find the throat pressures that are then put in Eq.͑25͒to obtain a series of linear equations for pore pressures only.Once the pore pressures have been found,Eq.͑26͒can be used to compute the throat pressures.Equation ͑25͒is solved using a standard iterative matrix solver.The pressures are used to compute the phase flow rates across the pore/throat boundaries.For example,from the equivalent resistors diagram of Fig.2͑a ͒,the water flow rate across the pore/throat boundary is the wetting layer flow rate and it is given through the following equation:Q water =͑P p −P t +P c ͒W eq ͑W P -lw +W T -l w ͓͒W eq −͑W P -c w +W T -c w͔͒.͑27͒The oil flow rate isQ oil =Q total −Q water .͑28͒Solving for constant injection rateWhile we set up the pressure equation for constant inlet and outlet pressures we want to simulate flow with a constant injection rate.Aker et al.͓25͔have shown how to achieve a constant injection rate in a dynamic network model.How-ever,their method involves solving for the pressure field twice at each time step.Here we use an approximate tech-nique that only requires a single pressure solution.Over time the pressure drop across the network changes.We assume that between time steps the change in pressure drop neces-sary to maintain a constant injection rate is small.Hence we simply adjust the pressure drop at each time step to maintain a constant value of Q as follows.͑a ͒For the n th time step,the pore and throat pressures are computed as described above,with a pressure drop ⌬P n .͑b ͒The total injection rate ͑Q n ͒is then obtained by summing all the flow rates between the inlet throats and their connected pores.͑c ͒For the next pressure solve,for the ͑n +1͒th time step we use a pressure drop⌬Pn +1=⌬P nϫͫ1+ͩQ desired −Q nQ desiredͪͬ,͑29͒where Q desired is the desired,target injection rate and is a constant parameter,which we set to 0.5.This method maintained Q n to within 0.5%of Q desired for the cases we studied.VI.SELECTION OF THE TIME STEPWe choose a time step ͑⌬t ͒according to the following formula:⌬t =min ͭ5ϫ10−5s,minͭV i2Q iͯi =1,...,n ͮ,͑30͒where n is the number of the elements in the pore network,V i is the i th element volume ͑i.e.,the maximum amount of fluid that can be held in the i th element ͒,and Q i is the total flow rate into the i th element.Equation ͑30͒ensures that an element cannot be completely filled in a single time step.The time step value of 5ϫ10−5s ensures that in most cases only a small fraction of a volume of a pore or throat is filled with invading fluid in a time step.VII.UPDATING FLUID VOLUMESOnce the fluid volumes are determined,the configuration of the phases in the elements is adjusted.There are two steps in this process.First,the pressure computation only deter-mines the total flow of oil and water between pore and throatM.S.AL-GHARBI AND M.J.BLUNT PHYSICAL REVIEW E 71,016308͑2005͒。
