PARALLEL FINITE ELEMENT METHODS FOR LARGE-SCALE

有限元法释义finit element method[计] 有限元法;finite element method有限元素法;点击人工翻译，了解更多人工释义实用场景例句全部The finite - element method is the most versatile.有限元法是最有用的方法.辞典例句The rectangular groove guide ( RGG ) is analyzed by finite element method ( FEM ) .本文用有限元法分析了矩形槽波导.互联网According to axi - symmetry of shaft, the semi - analytical finite element method is used.考虑到结构的轴对称性质, 分析时采用了半分析有限元法.互联网FEM FCT is used to solve three dimensional hypersonic inviscid flow.从Euler方程出发,利用流量修正有限元法(FEM?FCT)求解三维无粘流动的高速流场.互联网The result are compared with the finite element method's, and anastomosed. "本文计算结果与有限元法分析结果作了比较, 结果吻合较好.互联网Method: 3 - D finite element modeling was computed and analyzed.方法: 采用三维有限元法建立模型并计算、分析.互联网The cloth draping property is studied by using finite element method.以梁单元为模型,运用有限元法研究织物的悬垂性问题.互联网Therefore, the core question of rigid - viscoplastic finite element method has been solved.从而解决了刚粘塑性有限元法核心问题.互联网A permanent magnetic field for mono - crystal furnace was designed.采用有限元法设计了硅单晶炉用永磁磁体.互联网The edge finite element interpolation function of 1 - forms for prism is derived.就三梭柱单元导出了棱边有限元法的1 -形式线性插值基函数.互联网Methods: Three - dimensional finite element analysiswas adopted.方法: 采用三维有限元法.互联网A mixed method - NES FEM is systematically illustrated.提出了新型等效源法与有限元法的一种新耦合算法.互联网The finite element method was used for analysis of heat stress in chip on board ( COB ).本文采用有限元法分析了板上芯片( COB ) 的热应力分布.互联网Finite element method ( FEM ) occupies an important part in Computer Aided Engineering ( CAE ) methods.有限元法,也称有限单元法或有限元素法,在计算机辅助工程CAE 中占有重要的位置.互联网A finite element model for drawing process of high carbon wire with inner micro - defects was built.通过对高碳盘条内部缺陷的假设,采用有限元法研究了高碳盘条拉拔过程中,工艺参数对裂纹扩展情况的影响.互联网。

高性能计算在目标电磁散射特性分析中的应用刘阳;周海京;郑宇腾;陈晓洁;王卫杰;鲍献丰;李瀚宇【摘要】基于高性能计算的电磁数值模拟在目标电磁散射特性分析中发挥着越来越重要的作用.由于任一种数值方法都有一定的适用范围,不能高效处理所有问题,因此,有必要发展和集成多种数值方法,形成能够为不同类型问题的雷达散射截面(radar cross section,RCS)计算提供高效解决途径的软件系统.文中在并行自适应结构/非结构网格应用支撑软件框架之上,充分考虑数值方法的可扩展性以及物理个性的可分离性,通过基于机理、数据的混合可计算建模和接口设计,以及算法的模块化开发,发展了多种用于RCS计算的数值方法,并将其集成到高性能电磁数值模拟软件系统JEMS中.数值算例表明了JEMS具有高效分析多种目标电磁散射特性的能力,并在大规模并行计算方面具有显著优势.%The electromagnetic numerical simulation based on high performance computing gains more and more attention in analyzing the electromagnetic scattering characteristics of targets to meet the engineering increasing requirements. Since each method has its own advantages and disadvantages, and there is no one method which can deal with all problems, it is necessary to develop multi approach for integrating the software system, which can provide efficient means to analyze the electromagnetic scattering characteristics of different targets. Considering scalability of algorithms and separability of physical characteristics, based on parallel adaptive structured/unstructured mesh applications infrastructure, several numerical methods are developed and integrated into the electromagnetic numerical simulation software system, JEMS, with studying computable modeling, interface design andmodularized realization of algorithms. Some numerical examples illustrate JEMS has the capability in efficient solving the radar cross sections of different targets, and has advantages in large-scale parallel computing.【期刊名称】《电波科学学报》【年(卷),期】2019(034)001【总页数】9页(P3-11)【关键词】电磁散射;雷达散射截面;高性能计算;数值方法;并行支撑框架【作者】刘阳;周海京;郑宇腾;陈晓洁;王卫杰;鲍献丰;李瀚宇【作者单位】北京应用物理与计算数学研究所, 北京 100094;北京应用物理与计算数学研究所, 北京 100094;北京应用物理与计算数学研究所, 北京 100094;中物院高性能数值模拟软件中心, 北京 100088;中物院高性能数值模拟软件中心, 北京100088;中物院高性能数值模拟软件中心, 北京 100088;北京应用物理与计算数学研究所, 北京 100094【正文语种】中文【中图分类】O441引言目标电磁散射特性在雷达技术、目标识别、隐身与反隐身技术等应用中都有重要意义[1-4]. 电子技术的不断发展使它在军事和民用领域的应用日益拓展,以致目标电磁散射特性的数据获取与分析评估一直备受瞩目,建立在计算电磁学基础上的数值模拟技术为其提供了强有力的研究手段. 同时,各应用领域不断提高的实际工程需求,也为目标电磁散射特性的数值模拟提出了许多具有挑战性的问题,如超电大尺寸、复杂结构(包括深腔、缝隙、尖劈等)、复杂材质(非线性、各向异性、色散、时变媒质等)、宽频谱等[5-7]. 这些问题的求解不仅需要从数值算法设计的角度提高计算效率和精度,还需要从计算资源和并行技术的角度来增强对大规模计算的支撑. 近年来,计算机集群技术和并行计算技术的进步,促进电磁场问题的并行计算技术蓬勃发展,使基于高性能计算的电磁场数值模拟在实际工程应用中发挥着越来越重要的作用[8-10]. 许多商业软件,如CST、FEKO、HFSS等均提供并行版本,国内外很多科研团队也都针对不同的数值方法发展了各自的并行程序,有的甚至已形成了较为成熟的软件,如美国伊利诺伊大学的W. C. Chew教授的团队[8]、美国俄亥俄大学的J. F. Lee教授的团队[9],国内电子科技大学聂在平教授的团队[10]、北京理工大学盛新庆教授团队[11]、西安电子科技大学张玉教授的团队[12]等.由北京应用物理与计算数学研究所研制的并行自适应结构/非结构网格应用支撑软件框架(JASMIN/JAUMIN/JCOGIN)是针对科学计算中的结构/非结构网格应用,将高性能的数据结构进行了封装、并屏蔽了大规模并行和网格自适应的计算技术,能够支撑物理建模、数值方法、高性能算法的创新研究,可有效缩短基于现代高性能计算机的并行计算应用程序的研制时间[13]. 在该框架基础上,我们发展了高性能计算软件系统JEMS(J electro magnetic solver),用于多种电磁场问题的高效数值模拟. 本文将主要介绍JEMS中可用于目标电磁散射特性计算方面的内容,从各种数值算法及适用问题展开阐述,并通过介绍JEMS中针对不同类型问题的雷达散射截面计算的数值方法的研究进展和一系列数值算例,展示了JEMS具有高效分析多种目标电磁散射特性的能力,及其在大规模并行计算方面具有的优势.1 电磁散射的数值计算方法雷达散射截面[5](radar cross-section, RCS)是度量目标对电磁波散射能力的一个重要量化指标. RCS的定义为单位立体角内目标朝接收方向散射的功率与从给定方向入射于该目标的平面波功率密度之比的4π倍. 快速和精确获取目标的RCS成为衡量用于目标电磁散射特性研究数值方法有效性的关键.用于RCS计算的方法大致分为三类. 一类是解析方法,如Mie级数方法. 这类方法效率高且可得到问题的准确解,便于分析问题的物理本质,但适用范围太窄,不能满足复杂目标的分析需求.另一类是高频近似方法,如物理光学(physical optics, PO)、几何光学(geometrical optics, GO)、几何绕射理论(geometrical theory of diffraction, GTD)和物理绕射理论(physical theory of diffraction, PTD)等[14-16]. 高频近似方法计算速度快且对存储需求不高,特别在对电大尺寸目标的RCS计算中具有明显优势,能满足一定的工程需要. 然而对目标隐身与识别等应用,特别是含复杂结构或复杂材质的工程问题来说,该类方法的精度不够或无法求解.第三类是全波方法. 这类方法是目前计算电磁学的主流研究方向,如矩量法(method of moments, MoM)及其加速算法、有限元方法(finite element method, FEM)、时域有限差分法(finite difference time domain, FDTD)等[17-18],多用于处理电小或电中尺寸问题. 这类方法能够处理复杂目标,且给出较精确的数值解. MoM是基于积分方程的数值方法,积分方程中格林函数的使用,使无穷远处的辐射条件能够自然满足,场在数值网格中的传播过程得到精确描述,因此该方法的数值色散误差很小. 此外,MoM未知量数目较少且阻抗矩阵条件数较好. 然而,生成的阻抗矩阵是稠密的,造成矩阵元素的计算和存储以及矩阵方程的求解成为影响MoM求解能力的关键因素. 因此,其快速算法成为MoM重要的研究方向,如基于快速傅里叶变换的方法(CG-FFT、IE-FFT、AIM等)[19-20]、基于低秩矩阵压缩的纯代数方法(ACA、MLMDA等)[21-22]和基于快速多极子的方法(MLFMA)[23],有效解决了MoM的上述问题,使其在RCS计算中得到广泛使用. FEM[24]和FDTD[25]均是基于微分方程的方法. 这类方法通常算法简单,易于编程实现和程序并行化. 而且,FEM通用性强,可以处理复杂材质和结构,生成的矩阵具有稀疏性,但矩阵条件数较差. FDTD 方法是计算电磁学中被广泛使用的时域方法,具有宽频带瞬变电磁场分析计算的能力,适用于对宽带RCS的计算需求. 然而,这类方法在求解开的或无限大区域的问题时,需要辅以截断边界. 由于这类方法的未知量分布在整个传播空间,且为了保证所需的计算精度,在处理大尺寸和复杂结构时,往往需要较大的截断区域和精细的网格,从而造成巨大的未知量数目,导致其对计算机资源需求很大. 偏微分方程的局域性还造成这类方法中电磁场在数值网格的传播过程中形成较大的色散误差,导致其计算精度较差. 由于每种数值方法各具优点和劣势,因此将多种数值方法有效结合,取长补短发挥各自的优势,更好地高效求解RCS成为目前的研究热点之一．如全波方法之间的一种混合,即有限元边界积分(finite element boundary integral, FEBI)方法,它是有限元方法和积分方程方法的结合,能够有效消除FEM的截断误差,实现计算区域的最小化,同时具有处理复杂结构和材质的能力,其很强的实用性使其得到了深入发展. 此外,FEM和MoM的许多研究成果都能够应用到FEBI中[26]. 虽然在近几十年全波方法得到了系统的发展,各种快速算法、并行技术、矩阵求解加速技术等不断拓展了全波方法的求解能力,但是仍然有许多实际工程问题是全波方法无法有效或独立解决的. 因此,全波方法与高频方法的混合技术不可避免也成为一个备受关注的发展方向[5,27],包括MoM与PO、MoM与PTD、FEM与PO等,这类混合虽然由于高频近似方法的使用在一定程度上损失了计算精度,但是,它们不仅能够刻画电大目标上电小复杂结构,而且实现了较高的计算效率和较低的内存需求,在解决一些实际工程问题中成为能够折中考虑精度和效率的有效方法.综上所述,各种数值方法都有一定的适用范围,可以高效地求解一些问题. 然而,至今还未有哪种方法可以高效地处理所有问题,因此,有必要发展和集成用于RCS计算的多种数值方法,形成能够为不同类型问题的RCS计算提供高效解决途径的软件系统.2 电磁数值模拟软件系统JEMS目前,国防和高端商用领域迫切需要解决的复杂电磁工程问题,常常具有超电大尺寸、多尺度、多介质或复杂介质、多物理等特性. 基于高效能计算环境和并行支撑软件框架,我们将多种数值方法有机集成,发展了JEMS软件系统,用于电磁场问题的高效数值模拟. JEMS软件系统的设计,充分考虑了保持计算方法的持续可扩展性,并基于机理、数据的混合可计算建模及接口设计,保持物理个性的可分离性及可扩展性. 此外,由于并行支撑软件框架支持基于分布式内存和共享式内存的高性能计算,因此在该框架上发展的JEMS软件系统也支持上述两种高性能计算模式.JEMS软件系统的数值模拟能力并不仅限于目标散射特性分析,因而,本文在简单地整体回顾JEMS软件系统之后,将着重介绍JEMS中针对不同类型问题的RCS计算的解决方案和一系列数值算例,展示JEMS在大规模并行计算方面的优势.2.1 JEMS软件系统简介JEMS软件系统是基于并行自适应结构/非结构网格应用支撑软件框架(JASMIN/JAUMIN/JCOGIN)以算法模块联合研究的形式,与国内优势高校合作,充分发挥国内优势高校的研究力量,将国内外许多最新成果持续融入到软件平台的设计和研制中.综合考虑电磁场问题物理问题的特性、所关注的具体物理量,以及不同物理层次所需的模拟软件算法的共性基础构架的不同,发展的JEMS软件系统的软件体系如图1所示. 该软件系统的总体目标是通过突破在并行支撑框架上高效并行实现电磁脉冲源模拟、区域级/场景电磁模拟、电大多尺度结构全波电磁模拟以及多物理电磁计算等关键技术,在高性能计算环境中构建能力型电磁数值模拟软件系统,为具有明确应用牵引的高价值目标提供基于高性能计算的复杂电磁系统分析、优化及评估解决方案,为国内重大电磁工程问题快速定制高端专用计算平台[28].图1 电磁数值模拟软件系统JEMS体系图Fig.1 System diagram of electromagnetic numerical simulation software system JEMS用于目标电磁散射特性分析的多种数值方法属于平台级全波电磁模拟软件. 该软件包括时域和频域两部分内容,时域部分发展了基于HPA-adaptive模式的时域多算法求解技术,频域部分则采用基于非重叠区域分解的多种频域全波方法的混合集成技术,此外还发展了并行网格剖分技术、基于耦合波方法的电大馈线系统的快速计算技术以及电磁场/电路协同计算技术. 为典型的平台级目标(如飞行器等)构建了精确建模和电磁模拟能力,可实现目标近场和远场的多种电磁特性仿真数据. 此外,JEMS还包括电磁脉冲源模拟软件、区域级电磁模拟软件,以及器件级多物理电磁模拟软件.由于不同数值方法所需要的输入数据形式迥异,如网格数据、模型参数等,JEMS目前对基于不同数值方法发展的求解器的输入数据未做统一. 然而,JEMS中多种数值方法所需的网格数据均可由前处理引擎SuperMesh产生.2.2 用于RCS计算的不同数值方法的研究进展实际应用中需要进行电磁散射特性分析的目标从电尺寸、结构复杂度、材质以及频谱范围等方面都不尽相同,为从精度和效率两方面满足不同应用需求,JEMS软件系统提供多种算法供实际计算选择,包括MLFMA、FEM、PTD、FEBI-MLFMA-PO 以及FDTD等. 下面将逐一对其特点和适用范围进行介绍.2.2.1 多层快速多极子方法JEMS中的平台级频域全波电磁模拟软件JEMS-FD提供了基于组合场积分方程的MLFMA. 特别地,该方法通过高阶奇异值提取技术保证了算法的数值精度和计算稳定性,并提供块对角、稀疏近似逆等预条件技术保证超电大含腔目标的求解稳定性,可满足电大尺寸金属目标对应千万自由度矩阵方程的RCS高效求解. 算例1和算例2分别是使用JEMS中MLFMA对不同频率下F117隐身战机和含腔超电大目标的电磁散射特性分析.算例1 F117隐身战机不同频率下的电磁散射特性分析．模型如图2所示,入射平面波频率为1.5 GHz,入射方向沿机头正入射且采用垂直极化,模型电尺寸为88.8λ×60λ×10.6λ,λ为波长. 表面剖分的三角形网格数目97.6万,未知量数目146.5万,使用16个CPU核并行计算,计算时间为2.27 h,内存需求为7.9 GB,该频率下F117隐身战机的双站RCS如图3所示,与商业软件FEKO的结果吻合很好.当入射平面波频率为5.0 GHz时,模型电尺寸为310.8λ×210λ×37.1λ. 表面剖分的三角形网格数目为996.8万,未知量数目 1 495.2万,使用10个CPU核并行计算,计算时间约5.5 h,内存需求约为84.3 GB,图4给出该频率下F117隐身战机的双站RCS的模拟结果.图2 F117隐身战机模型Fig.2 F117 model图3 频率1.5 GHz时F117的双站RCSFig.3 Bistatic RCS of F117model(frequency=1.5 GHz)图4 频率5 GHz时F117的双站RCSFig.4 Bistatic RCS of F117model(frequency=5 GHz)算例2 含腔超电大目标的电磁散射特性分析. 模型如图5所示,入射平面波频率为0.9 GHz,入射方向沿机头正入射且采用垂直极化,模型电尺寸为66λ×48λ×20λ. 网格剖分的未知量数目约118万,计算时间13 181 s,内存需求为6.7 GB,此含腔超电大目标的双站RCS如图6所示.图5 含腔超电大目标模型Fig.5 Model for the electrical large target with a cavity图6 频率0.9 GHz时含腔超电大目标的双站RCSFig.6 Bistatic RCS of the electrical large target with a cavity(frequency=0.9 GHz)2.2.2 有限元方法在频域全波方法中,还发展了针对复杂多尺度、多材料(包括介质、金属、吸波材料、频变材料、各项异性材料等)结构的FEM,可支持多种激励源(如平面波、高斯波束、点源、波导激励源、电压/电流源等),采用非结构网格并行自适应加密技术和区域分解求解技术,具有数万CPU核的并行扩展能力,可实现对数亿网格规模复杂目标的RCS分析. 算例3和算例4分别是使用JEMS中FEM分析频率选择表面和舰船模型的电磁散射特性.算例3 频率选择表面的电磁散射特性分析. 模型如图7所示含1 000个单元. 入射平面波频率0.3 GHz,入射方向沿-Z轴(即垂直于频率选择表面),极化方向沿+X轴. 模型电尺寸为数十个波长,四面体网格数目为414万,采用8个CPU核并行,区域分解迭代步数为8. 如图8中所示,JEMS中FEM获得的该模型的双站RCS计算结果与商业软件HFSS的一致.图7 频率选择表面的模型Fig.7 Model for frequency selective surface图8 频率0.3 GHz时频率选择表面的双站RCSFig.8 Bistatic RCS of the frequency selective surface(frequency=0.3 GHz)算例4 舰船模型的电磁散射特性分析. 模型如图9所示,尺寸为130.8 m×20m×23.1 m. 入射平面波频率为1 GHz,入射方向的俯仰角为45°,方位角为0°,且为水平极化. 四面体网格规模约为3亿,在天河-2超级计算机上启动400个进程,共计9 600CPU核完成自适应计算. 图10是舰船模型在频率1 GHz时的双站RCS.图9 舰船模型Fig.9 The ship model图10 频率1 GHz时舰船的双站RCSFig.10 Bistatic RCS of theship(frequency=1 GHz)2.2.3 物理绕射理论目标的电尺寸越大,其表面散射场的局部效应越明显,可利用高频方法的局部性原理来求解其散射场. JEMS中提供了PTD方法,通过考虑边缘的绕射电流达到对PO方法的修正,以提高其计算精度. 另外,采用深度缓冲器(z-buffer)算法判断遮挡,区分物体表面的照射和非照射区域,从而实现对超电大尺寸金属和多层涂覆目标的RCS计算. 算例5和算例6是采用JEMS中PTD对金属舰船模型以及涂覆介质材料的舰船模型的电磁散射特性分析.算例5 舰船模型的电磁特性分析. 仍然考虑算例4中的舰船模型. 入射平面波的频率为0.3 GHz,且采用垂直极化,当入射方向的俯仰角为90°,方位角从0°扫描到360°时,JEMS中PTD计算的舰船模型的单站RCS与商业软件CST中的SBR方法的结果如图11所示,二者吻合得较好.图11 频率0.3 GHz时舰船的单站RCSFig.11 Monostatic RCS of theship(frequency=0.3 GHz)算例6 涂覆舰船模型的电磁特性分析. 仍采用算例4中的舰船模型,表面共涂覆三层介质,表1中给出其相对介电常数、相对磁导率,以及厚度等参数. 入射平面波频率为3 GHz,入射方向的俯仰角为90°,方位角从0°扫描到360°. 图12是CST软件的PO方法与JEMS中PTD方法的计算结果对比.表1 涂层介质材料的参数Tab.1 Material parameters for dielectric coats层号相对介电常数相对磁导率涂层厚度/mm 11514.412-j12.3531.02 2151-j5.2421.77 34.254-j2.3311.96图12 频率3 GHz时涂覆舰船的单站RCSFig.12 Monostatic RCS of the coated ship(frequency=3 GHz)2.2.4 全波与高频混合方法最近,针对含金属/介质混合局部结构的电大尺寸问题的RCS分析,JEMS还研发了迭代型全波与高频混合方法FEBI-MLFMA-PO,充分利用FEBI处理复杂结构和材质的能力,以及PO方法处理电大平滑目标的高效性. 通过MLFMA实现对全波算法部分的加速,并采用自适应交叉近似方法提高全波与高频区域相互作用子矩阵的计算效率. 全波与高频区域的耦合子矩阵为稠密阵,采用自适应交叉近似方法可有效降低计算复杂度和内存需求,该算法主要包括求一行或一列的最大值、计算矩阵元素以及每步的误差．在JEMS中,将整个计算区域划分成多个块,求一行或一列中最大值转化为并行求出每一块中最大值,再通过比较块的最大值找出一行或一列的最大值;矩阵元素则是在每一块上并行计算;每步的误差则是先通过每块上计算所属部分的值,而后通过归约计算得到总的每步误差. 在保证一定精度的前提下,有效减少了未知量数目,降低了计算复杂度. 算例7是使用JEMS中FEBI-MLFMA-PO方法分析观察室内含介质体的舰船电磁散射特性.算例7 观察室内含介质体的舰船电磁散射特性分析. 模型如图13,观察室内介质体的相对介电常数为1.5,尺寸3 m×2.5 m×2.0 m．入射平面波频率为50 MHz,入射方向的俯仰角为45°,方位角为0°,且为水平极化. 网格剖分40 109个四面体,以及9 956个三角形(如果全部使用FEBI,则网格剖分含40 109个四面体,以及58 778个三角形),有效减少了未知量数目. 图14给出了利用JEMS中的FEBI-MLFMA-PO,商业软件FEKO中的全波方法MLFMA和混合方法MoM-PO三种方法的计算结果比较．可以看出,在前向和后向附近,与FEKO的MoM-PO混合方法相比,JEMS 中的FEBI-MLFMA-PO的结果与FEKO全波方法MLFMA的结果吻合更好.图13 观察室内含介质体的舰船模型Fig.13 Ship model with a cabin having dielectric object图14 观察室内含介质体的舰船的双站RCSFig.14 Bistatic RCS of the ship witha cabin having dielectric object2.2.5 时域有限差分方法此外,考虑到一些工程问题中对宽带RCS的计算需求,JEMS中的平台级时域全波电磁模拟软件JEMS-TD提供FDTD方法计算宽带RCS的功能. 应用FDTD计算瞬态近场,再由时域近远场外推公式得到特定频率的远场信息,为提高计算效率和精度,特别开发了混合阶和非均匀网格技术. 算例8给出JEMS中FDTD计算的整机模型的RCS.算例8 整机电磁散射特性分析. 整机尺寸为35 m×38 m×12 m,机身为全金属半硬壳式结构,包括四段机身结构、有机玻璃机头罩、起落架及发动机等结构. 入射波频率为1 GHz,沿机头正入射,且采用垂直极化. 利用FDTD计算该飞机模型的水平面和垂直面的双站RCS,六面体网格剖分规模约300亿,使用10 800个CPU核,计算结果如图15～16,并与CST中SBR进行了对比.图15 水平面上飞机的双站RCSFig.15 Bistatic RCS of airplane on horizontal plane图16 垂直面上飞机的双站RCSFig.16 Bistatic RCS of airplane on vertical plane3 结论本文从工程应用中目标电磁散射特性分析遇到的许多难题引出发展基于高性能计算的电磁数值方法的重要性. 首先回顾了用于RCS计算的三类方法,通过分析每种数值方法的利弊,阐明了它们具有不同的适用范围.由于没有一种数值方法能够同时解决所有问题,为从精度和效率两方面满足不同应用需求,需通过发展不同算法供实际计算选择. 本文着重介绍了以这种思路为指导的基于并行支撑框架JASMIN/JAUMIN/JCOGIN的高性能计算软件系统JEMS. JEMS本身的功能很多,这里只介绍其中针对不同类型问题的雷达散射截面计算的数值方法的研究进展,并通过一些相关算例展示出JEMS具有分析多种类型目标电磁散射特性方面的能力以及其在大规模并行计算方面的优势. 实际上,JEMS的研发团队持续通过算法模块形式,将国内外计算电磁学的最新成果融入到软件系统当中,期待通过不断丰富算法功能、优化算法效率为国内重大电磁工程问题提供基于高性能计算的复杂电磁系统分析、优化及评估解决方案.参考文献【相关文献】[1] 黄培康, 殷红成, 许小剑. 雷达目标特性[M]. 北京: 电子工业出版社, 2005.[2] 庄钊文, 袁乃昌, 莫锦军, 等. 军用目标雷达散射截面预估与测量[M]. 北京: 科学出版社, 2007.[3] 保铮, 邢孟道, 王彤. 雷达成像技术[M]. 北京: 电子工业出版社, 2005.[4] 阮颖铮. 雷达散射截面与隐身技术[M]. 北京: 国防工业出版社, 1998.[5] 聂在平, 方大纲. 目标与环境电磁散射特性建模——理论、方法与实现[M]. 北京: 国防工业出版社, 2009.[6] 桑建华. 飞行器隐身技术[M]. 北京:航空工业出版社, 2013.[7] 艾俊强, 周莉, 杨青真. S弯隐身喷管[M]. 北京: 国防工业出版社, 2017.[8] SONG J M, LU C C, CHEW W C, et al. Fast illinois solver code (FISC) [J]. IEEE antennas and propagation magazine, 1998, 40(3): 27-34.[9] PENG Z, LIM K H, LEE J F. Non-conformal domain decomposition method for solving large multi-scale electromagnetic scattering problem[J]. Proceedings of the IEEE, 2013, 101(12): 298-319.[10] 胡俊, 聂在平, 王军, 等. 三维电大尺寸目标电磁散射求解的多层快速多极子方法[J]. 电波科学学报, 2004, 19(5): 509-514.HU J, NIE Z P, WANG J, et al. Multilevel fast multipole algorithm for solving scattering from 3-D electrically large object[J]. Chinese journal of radio science, 2004, 19(5): 509-514. (in Chinese)[11] 潘小敏, 盛新庆. 电特大复杂目标电磁特性的高效精确并行计算[J]. 电波科学学报, 2008, 23(5): 888-891.PAN X M, SHENG X Q. Efficient and accurate parallel computation of electromagnetic scattering by extremely large targets[J]. Chinese journal of radio science, 2008, 23(5): 888-891.(in Chinese)[12] ZHANG Y, ZHAO X W, DONORO D G, et al. Parallelized hybrid method with higher-order MoM and Po for analysis of phased array antennas on electrically large platforms[J]. IEEE transactions on antennas and propagation, 2010, 58(2): 4110-4115.[13] MO Z Y, ZHANG A Q, CAO X L, et al. JASMIN: A software infrastructure for large scale parallel adaptive structured mesh application[J]. Frontiers of computer science in China, 2010, 4(4): 480-488.[14] KLINE M, KAY I. Electromagnetic theory and geometrical optics[M]. New York: Wiley Inter-science, 1965.[15] KELLER J B. A geometrical theory of diffraction[M]. New York: Mc Graw-hill Book Co.,。

Parallel FEM Simulation of Crack Propagation–Challenges,Status,and PerspectivesList of Authors(in alphabetical order):Bruce Carter,Chuin-Shan Chen,L.Paul Chew,Nikos Chrisochoides,Guang R.Gao,Gerd Heber,Antony R.Ingraffea,Roland Krause,Chris Myers,Demian Nave,Keshav Pingali,Paul Stodghill,Stephen Vavasis,Paul A.WawrzynekCornell Fracture Group,Rhodes Hall,Cornell University,Ithaca,NY14853CS Department,Upson Hall,Cornell University,Ithaca,NY14853CS Department,University of Notre Dame,Notre Dame,IN46556EECIS Department,University of Delaware,Newark,DE197161IntroductionUnderstanding how fractures develop in materials is crucial to many disciplines,e.g.,aeronautical engineering,mate-rial sciences,and geophysics.Fast and accurate computer simulation of crack propagation in realistic3D structures would be a valuable tool for engineers and scientists exploring the fracture process in materials.In this paper,we will describe a next generation crack propagation simulation software that aims to make this potential a reality.Within the scope of this paper,it is sufficient to think about crack propagation as a dynamic process of creating new surfaces within a solid.During the simulation,crack growth causes changes in the geometry and,sometimes,in the topology of the model.Roughly speaking,with the tools in place before the start of this project,a typical fracture analysis at a resolution of degrees of freedom,using boundary elements,would take about100hours on a state-of-the-art single processor workstation.The goal of this project is it to create a parallel environment which allows the same analysis to be done,usingfinite elements,in1hour at a resolution of degrees of freedom.In order to attain this level of performance,our system will have two features that are not found in current fracture analysis systems:Parallelism–Current trends in computer hardware suggest that in the near future,high-end engineering worksta-tions will be8-or16-way SMP“nodes”,and departmental computational servers will be built by combining a number of these nodes using a high-performance network switch.Furthermore,the performance of each processor in these nodes will continue to grow.This will happen not only because of faster clock speeds,but also becausefiner-grain parallelism will be exploited via multi-way(or superscalar)execution and multi-threading.Adaptivity–Cracks are(hopefully)very small compared with the dimension of the structure,and their growth is very dynamic in nature.Because of this,it is impossible to know a priori howfine a discretization is required to accurately predict crack growth.While it is possible to over-refine the discretization,this is undesirable,as it tends to dramatically increase the required computational resources.A better approach is to adaptively choose the discretization refinement.Initially,a coarse discretization is used,and,if this induces a large error for certain regions of the model,then the discretization is refined in those regions.The dynamic nature of crack growth and the need to do adaptive refinement make crack propagation simulation a highly irregular application.Exploiting parallelism and adaptivity presents us with three major research challenges, developing algorithms for parallel mesh generation for unstructured3D meshes with automatic element sizecontrol and provably good element quality,implementing fast and robust parallel sparse solvers,anddetermining efficient schemes for automatic,hybrid h-p refinement.To tackle the challenges of developing this system,we have assembled a multi-disciplinary and multi-institutional team that draws upon a wide-ranging pool of talent and the resources of3universities.2System OverviewFigure1gives an overview of a typical simulation.During pre-processing,a solid model is created,problem specificFigure1:Simulation loop.boundary conditions(displacements,tractions,etc.)are imposed,andflaws(cracks)are introduced.In the next step, a volume mesh is created,and(linear elasticity)equations for the displacements are formulated and solved.An error estimator determines whether the desired accuracy has been reached,or further iterations,after subsequent adaptation, are necessary.Finally,the results are fed back into a fracture analysis tool for post-processing and crack propagation.Figure1presents the simulation loop of our system in itsfinal and most advanced form.Currently,we have se-quential and parallel implementations of the outer simulation loop(i.e.,not the inner refinement loop)running with the following restrictions:currently,the parallel mesher can handle only polygonal(non-curved)boundaries,which can be handled by the sequential meshers though(see section3).We have not yet implemented unstructured h-refinement and adaptive p-refinement,although the parallel formulator can handle arbitrary p-order elements.3Geometric Modeling and Mesh GenerationThe solid modeler used in the project is called OSM.OSM,as well as the main pre-and post-processing tool, FRANC3D,is freely available from the Cornell Fracture Group’s website[12].FRANC3D-a workstation based FRacture ANalysis Code for simulating arbitrary non-planar3D crack growth-has been under development since 1987,with hydraulic fracture and crack growth in aerospace structures as the primary application targets since its inception.While there are a few3D fracture simulators available and a number of other software packages that can model cracks in3D structures,these are severely limited by the crack geometries that they can represent(typically planar elliptical or semi-elliptical only).FRANC3D differs by providing a mechanism for representing the geometry and topology of3D structures with arbitrary non-planar cracks,along with functions for1)discretizing or meshing the structure,2)attaching boundary conditions at the geometry level and allowing the mesh to inherit these values,and3) modifying the geometry to allow crack growth but with only local re-meshing required to complete the model.The simulation process is controlled by the user via a graphic user-interface,which includes windows for the display of the3D structure and a menu/dialogue-box system for interacting with the program.The creation of volume meshes for crack growth studies is quite challenging.The geometries tend to be compli-cated because of internal boundaries(cracks).The simulation requires smaller elements near each crack front in order to accurately model high stresses and curved geometry.On the other hand,larger elements might be sufficient awayfrom the crack front.There is a considerable difference between these two scales of element sizes,which amounts to three orders of magnitude in real life applications.A mesh generator must provide automatic element size control and give certain quality guarantees for elements.The mesh generators we studied so far are QMG by Steve Vavasis[15], JMESH by Joaquim Neto[11],and DMESH by Paul Chew[15].These meshers represent three different approaches: octree-algorithm based(QMG),advancing front(JMESH),and Delaunay mesh(DMESH).QMG and DMESH come with quality guarantees for elements in terms of aspect ratio.All these mesh generators are sequential and give us insight into the generation of large“engineering quality”meshes.We decided to pursue the Delaunay mesh based approachfirst for a parallel implementation,which is described in[5].Departing from traditional approaches,we simultaneously do mesh generation and partitioning in parallel.This not only eliminates most of the overhead of the traditional approach,it is almost a necessary condition to do crack growth simulations at this scale,where it is not always possible or too expensive to keep up with the geometry changes by doing structured h-refinement.The implementation is a parallelization of the so-called Bowyer-Watson(see the references in[5])algorithm:given an initial Delaunay triangulation,we add a new point to the mesh,determine the simplex containing this point and the point’s cavity(the union of simplices with non-empty circumspheres),and,fi-nally,retriangulate this cavity.One of the challenges for a parallel implementation is that this cavity might extend across several submeshes(and processors).What looks like a problem,turns out to be the key element in unifying mesh generation and partitioning:the newly created elements,together with an adequate cost function,are the best candidates to do the“partitioning on thefly”.We compared our results with Chaco and MeTis in terms of equidistri-bution of elements,relative quality of mesh separators,data migration,I/O,and total performance.Table1shows a runtime comparison between ParMeTis with PartGeomKway(PPGK)and,our implementation,called SMGP,on16 processors of an IBM SP2for meshes of up to2000K elements.The numbers behind SMGP refer to different cost functions used in driving the partitioning[5].Mesh Size PPGK SMGP0SMGP1SMGP2SMGP3200K9042424242500K215658764621000K4399716091942000K1232133310110135Table1:Total run time in seconds on16processors.4Equation Solving and PreconditioningWe chose PETSc[2,14]as the basis for our equation solver subsystem.PETSc provides a number of Krylov space solvers,and a number of widely-used preconditioners.We have augmented the basic library with third party packages, including BlockSolve95[8]and the Barnard’s SPAI[3].In addition,we have implemented a parallel version of the Global Extraction Element-By-Element(GEBE)preconditioner[7](which is unrelated to the EBE preconditioner of Winget and Hughes),and added it to the collection using PETSc’s extension mechanisms.The central idea of GEBE is to extract subblocks of the global stiffness matrix associated with elements and invert them,which is highly parallel.The ICC preconditioner is frequently used in practice,and is considered to be a good preconditioner for many elasticity problems.However,we were concerned that it would not scale well to the large number of processors required for ourfinal system.We believed that GEBE would provide a more scalable implementation,and we hoped that it would converge nearly as well as ICC.In order to test our hypothesis,we ran several experiments on the Cornell Theory Center SP-2.The preliminary performance results for the gear2and tee2models are shown in Tables2and3,respectively.(gear2is a model of a power transmission gear with a crack in one of its teeth.tee2is a model of a T steel profile).For each model,we ran the Conjugant Gradient solver with both BlockSolve95’s ICC preconditioner and our own parallel implementation of GEBE on8to64processors.(The iteration counts in Tables2and3correspond to a reduction of the residualerror,which is completely academic at this point.)The experimental results confirm our hypothesis: GEBE converges nearly as quickly as ICC for the problems that we tested.Our naive GEBE implementation scales much better than BlockSolve95’s sophisticated ICC implementation.Prec.Nodes Prec Time Per Iters.TypeTime(s)Iteration (s)ICC817.080.2416GEBE89.430.19487ICC1615.470.27422GEBE16 6.710.11486ICC328.510.32539GEBE32 3.730.08485ICC6411.000.28417GEBE 64 4.740.07485Table 2:Gear2(79,656unknowns)Prec.Nodes Prec Time Per Iters.Type Time(s)Iteration (s)ICC 3230.000.292109GEBE 3235.700.212421ICC 6423.600.292317GEBE 647.600.122418Table 3:Tee2(319,994unknowns)5AdaptivityUnderstanding the cost and impact of the different adaptivity options is the central point in our current activities.We are in the process of integrating the structured (hierarchical)h-re finement into the parallel testbed,and the final version of this paper will contain more results on that.Our implementation follows the approach of Biswas and Oliker [4]and currently handles tetrahedra,while allowing enough flexibility for an extension to non-tetrahedral element types.Error Estimation and Adaptive Strategies.For relatively simple,two-dimensional problems,stress intensity factors can be computed to an accuracy suf ficient for engineering purposes with little mesh re finement by proper use of sin-gularly enriched elements.There are many situations though when functionals other than stress intensity factors are of interest or when the singularity of the solution is not known a priori .In any case the engineer should be able to evaluate whether the data of interest have converged to some level of accuracy considered appropriate for the compu-tation.It is generally suf ficient to show,that the data of interest are converging sequences with respect to increasing degrees of freedom.Adaptive finite element methods are the most ef ficient way to achieve this goal and at the same time they are able to provide estimates of the remaining discretization error.We de fine the error of the finite elementsolutionasand a possible measure for the discretization error is the energynorm,Following an idea of Babu ˇs ka and Miller [1]the error estimator introduced by Kelly et.al.[9,6]is derived by inserting the finite element solution into the original differential equation system and calculating a norm of the residual using interpolation estimates.An error indicator computable from local results of one element of the finite element solution is then derived and the corresponding error estimator is computed by summing the contribution of the error indicators over the entire domain.The error indicator is computed with a contribution from the interior residual of the element and a contribution of the stress jumps on the faces of an element.Details on the computation of the error estimator from the finite element solution can be found in [10].Control of a Mesh Generator.For a sequence of adaptively re fined and quasi optimal meshes,the rate of conver-gence is independent of the smoothness of the exact solution.A mesh is called quasi optimal if the error associated with each element is nearly the same.The goal of an adaptive finite element algorithm is to generate a sequence of quasi optimal meshes by equilibrating the estimated error until a prescribed accuracy criterion is reached.Starting from an initial mesh,error indicators and the error estimator are computed in the post-processing step of the solutionphase.The idea is then to compute for each element the new elementsizefrom the estimated error,the present elementsize and the expected rate of convergence.6Future WorkThe main focus of our future work will be on improving the performance of the existing system.The Cornell Fracture Group is continuously extending our test-suite with new real world problems.We are considering introducing special elements at the crack tip,and non-tetrahedral elements (hexes,prisms,pyramids)elsewhere.The linear solver,after proving its robustness,will be wrapped into a Newton-type solver for nonlinear problems.Among the newly introduced test problems are some that can be made arbitrarily ill-conditioned (long thin plate or tube models with cracks in them)in order to push the iterative methods to their limits.We are exploring new precon-ditioners (e.g.,support tree preconditioning),and multigrid,as well as sparse direct solvers,to make our environmentmore effective and robust.We have not done yet any specific performance tuning,like locality optimization.This is not only highly platform dependent,but also has to be put in perspective to the forthcoming runtime optimizations,like dynamic load balancing. We are following with interest the growing importance of latency tolerant architectures in the form of multithreading and exploring for which parts of the project multithreaded architectures are the most beneficial.Finally,there is a port of our code base to the new256node NT cluster at the Cornell Theory Center underway. 7ConclusionsAt present,our project can claim two major contributions.Thefirst is our parallel mesher/partition,which is thefirst practical implementation of its kind with quality guarantees.This technology makes it possible,for thefirst time,to fully automatically solve problems using unstructured h-refinement in a parallel setting.The second major contribution is to show that GEBE outperforms ICC,at least for our problem class.We have shown that,not only does GEBE converge almost as quickly as ICC,it is much more scalable in a parallel setting than ICC.We believe that GEBE,not ICC,is the yardstick against which other parallel preconditioners should be measured.Andfinally,ourfirst experimants indicate that we should be able to meet our project’s performance goals.We are confident that,as we run our system on larger and faster machines,as we further optimize each of the subsystems,and as we incorporate adaptive h-and p-refinement,we will reach our performance goals.References[1]I.Babuˇs ka and ler,“A-posteriori error estimates and adaptive techniques for thefinite element method,”Technical Report Tech.Note BN–968,University of Maryland,Inst.for Physics,Sci.and Tech.,1981.[2]S.Balay,W.D.Gropp,L.Curfman McInnes,and B.F.Smith,“Efficient management of parallelism in object-oriented numerical software libaries”,In E.Arge,A.M.Bruaset,and ngtangen,editors,Modern Software Tools in Scientific Computing,Birkhauser Press,1997.[3]S.T.Barnard and R.Clay,“A portable MPI implementation of the SPAI preconditioner in ISIS++”,Eighth SIAMConference for Parallel Processing for Scientific Computing,March1997.[4]R.Biswas and L.Oliker,“A new procedure for dynamic adaption of three-dimensional unstructured grids”,AppliedNumerical Mathematics,13:437–452,1994.[5]N.Chrisochoides and D.Nave,“Simultaneous mesh generation and partitioning for Delaunay meshes”,In8thInt’l.Meshing Roundtable,1999.[6]J.P.de S.R.Gago,D.W.Kelly,O.C.Zienkiewicz,and I.Babuˇs ka,“A posteriori error analysis and adaptive processesin thefinite element method:Part II–Adaptive mesh refinement”,International Journal for Numerical Methods in Engineering,19:1621–1656,1983.[7]I.Hladik,M.B.Reed,and G.Swoboda,“Robust preconditioners for linear elasticity FEM analyses”,InternationalJournal for Numerical Methods in Engineering,40:2109–2127,1997.[8]M.T.Jones and P.E.Plassmann,“Blocksolve95users manual:Scalable library software for the parallel solution ofsparse linear systems”,Technical Report ANL-95/48,Argonne National Laboratory,December1995.[9]D.W.Kelly,J.P.de S.R.Gago,O.C.Zienkiewicz,and I.Babuˇs ka,“A posteriori error analysis and adaptive processesin thefinite element method:Part I–Error analysis”,International Journal for Numerical Methods in Engineering, 19:1593–1619,1983.[10]R.Krause,“Multiscale Computations with a Combined-and-Version of the Finite Element Method”,PhDthesis,Universit¨a t Dortmund,1996.[11]o et al.,“An Algorithm for Three-Dimensional Mesh Generation for Arbitrary Regions with Cracks”,submitted for publication.[12]/.[13]/People/chew/chew.html.[14]/petsc/index.html.[15]/vavasis/vavasis.html.。

MSc in Mechanical Engineering Design MSc in Structural Engineering LECTURER: Dr. K. DAVEY(P/C10)Week LectureThursday(11.00am)SB/C53LectureFriday(2.00pm)Mill/B19Tut/Example/Seminar/Lecture ClassFriday(3.00pm)Mill/B192nd Sem. Lab.Wed(9am)Friday(11am)GB/B7DeadlineforReports1 DiscreteSystems DiscreteSystems DiscreteSystems2 Discrete Systems. Discrete Systems. Tutorials/Example I.Meshing I.Deadline 3 Discrete Systems Discrete Systems Tutorials/Example IIStart4 Discrete Systems. Discrete Systems. DiscreteSystems.5 Continuous Systems Continuous Systems Tutorials/Example II. Mini Project6 Continuous Systems Continuous Systems Tutorials/Example7 Continuous Systems Continuous Systems Special elements8 Special elements Special elements Tutorials/Example III.Composite IIDeadline *9 Special elements Special elements Tutorials/Example10 Vibration Analysis Vibration Analysis Vibration Analysis III Deadline11 Vibration Analysis Vibration Analysis Tutorials/Example12 VibrationAnalysis Tutorials/Example Tutorials/Example13 Examination Period Examination Period14 Examination Period Examination Period15 Examination Period Examination Period*Week 9 is after the Easter vacation Assignment I submission (Box in GB by 3pm on the next workingday following the lab.) Assignment II and III submissions (Box in GB by 3pm on Wed.)CONTENTS OF LECTURE COURSEPrinciple of virtual work; minimum potential energy.Discrete spring systems, stiffness matrices, properties.Discretisation of a continuous system.Elements, shape functions; integration (Gauss-Legendre).Assembly of element equations and application of boundary conditions.Beams, rods and shafts.Variational calculus; Hamilton’s principleMass matrices (lumped and consistent)Modal shapes and time-steppingLarge deformation and special elements.ASSESSMENT: May examination (70%); Short Lab – Holed Plate (5%); Long Lab – Compositebeam (10%); Mini Project – Notched component (15%).COURSE BOOKSBuchanan, G R (1995), Schaum’s Outline Series: Finite Element Analysis, McGraw-Hill.Hughes, T J R (2000), The Finite Element Method, Dover.Astley, R. J., (1992), Finite Elements in Solids and Structures: An Introduction, Chapman &HallZienkiewicz, O.C. and Morgan, K., (2000), Finite Elements and Approximation, DoverZienkiewicz, O C and Taylor, R L, (2000), The Finite Element Method: Solid Mechanics,Butterworth-Heinemann.IntroductionThe finite element method (FEM) is a numerical technique that can be applied to solve a range of physical problems. The method involves the discretisation of the body (domain) of interest into subregions, which are known as elements. This enables a continuum problem to be described by a finite system of equations. In the field of solid mechanics the FEM is undoubtedly the solver of choice and its use has revolutionised design and analysis approaches. Many commercial FE codes are available for many types of analyses such as stress analysis, fluid flow, electromagnetism, etc. In fact if a physical phenomena can be described by differential or integral equations, then the FE approach can be used. Many universities, research centres and commercial software houses are involved in writing software. The differences between using and creating code are outlined below:(A) To create FE software1. Confirm nature of physical problem: solid mechanics; fluid dynamics; electromagnetic; heat transfer; 1-D, 2-D, 3-D; Linear; non-linear; etc.2. Describe mathematically: governing equations; loading conditions.3. Derive element equations: convert governing equations into algebraic form; select trial functions; prepare integrals for numerical evaluation.4. Assembly and solve: assemble system of equations; application of loads; solution of equations.5. Compute:6. Process output: select type of data; generate related data; display meaningfully and attractively.(B) To use FE software1. Define a specific problem: geometry; physical properties; loads.2. Input data to program: geometry of domain, mesh generation; physical properties; loads-interior and boundary.3. Compute:4. Process output: select type of data; generate related data; display meaningfully and attractively.DISCRETE SYSTEMSSTATICSThe finite element involves the transformation of a continuous system (infinite degrees of freedom) into a discrete system (finite degrees of freedom). It is instructive therefore to examine the behaviour of simple discrete systems and associated variational methods as this provides real insight and understanding into the more complicated systems arising from the finite element method.Work and Strain energyFLuxConsider a metal bar of uniform cross section, A , fixed at one end (unrestrained laterally) and subjected to an axial force, F , at the other.Small deflection theory is assumed to apply unless otherwise stated.The work done, W , by the applied force F is .a ()∫′′=uau d u F WIt is worth mentioning at this early stage that it is not always possible to express work in this manner for various reasons associated with reversibility and irreversibility. (To be discussed later)The work done, W , by the internal forces, denoted strain energy , is se22200se ku 21u L EA 2121EAL d EAL d AL W ==ε=ε′ε′=ε′σ=∫∫εεwhere ε=u L and stiffness k EA L=.The principle of virtual workThe principle of virtual work states that the variation in strain energy is equal to the variation in the work done by applied forces , i.e.()u F u u d u F du d W u ku u ku 21du d ku 21W u0a 22se δ=δ⎟⎟⎠⎞⎜⎜⎝⎛′′=δ=δ=δ⎟⎠⎞⎜⎝⎛=⎟⎠⎞⎜⎝⎛δ=δ∫()0u F ku =δ−⇒Note that use has been made of the relationship δf dfduu =δ where f is an arbitrary functional of u . In general displacement u is a function of position (x say) and it is understood that ()x u δ means a change in ()u x with xfixed. Appreciate that varies with from zero to ()'u F 'u ()u F F = in the above integral.Bearing in mind that δ is an arbitrary variation; then this equation is satisfied if and only if F , which is as expected. Before going on to apply the principle of virtual work to a continuous system it is worth investigating discrete systems further. This is because the finite element formulation involves the transformation of a continuous system into a discrete one. u ku =Spring systemsConsider a single spring with stiffness independent of deflection. Then, 2F21u1F1u2k()()⎟⎟⎠⎞⎜⎜⎝⎛⎥⎦⎤⎢⎣⎡−−=−=2121212se u u k k k k u u 21u u k 21W()()()⎟⎟⎠⎞⎜⎜⎝⎛⎥⎦⎤⎢⎣⎡−−δδ=δ−δ−=δ21211212se u u k k k k u u u u u u k W()⎟⎟⎠⎞⎜⎜⎝⎛δδ=δ+δ=δ21212211a F F u u u F u F W , where ()111u F F = and ()222u F F =.Note here that use has been made of the relationship δ∂∂δ∂∂δf f u u f u u =+1122, where f is an arbitrary functional of and . Observe that in this case is a functional of 1u u 2W se u u u 2121=−, so()()(121212*********se se u u u u k u u ku 21du d u du dW W δ−δ−=−δ⎟⎠⎞⎜⎝⎛=δ=δ).The principle of virtual work provides,()()()()0F u u k u F u u k u 0W W 21221121a se =−−δ+−−−δ⇒=δ−δand since δ and δ are arbitrary we have. u 1u 2F ku ku 11=−2u 2 F ku k 21=−+represented in matrix form,u F K u u k k k k F F 2121=⎟⎟⎠⎞⎜⎜⎝⎛⎥⎦⎤⎢⎣⎡−−=⎟⎟⎠⎞⎜⎜⎝⎛=where K is known as the stiffness matrix . Note that this matrix is singular (det K k k =−=220) andsymmetric (K K T=). The symmetry is a result of the fact that a unit deflection at node 1 results in a force at node 2 which is the same in magnitude at node 1 if node 2 is moved by the same amount.Could also have arrived at equation above via()⎟⎟⎠⎞⎜⎜⎝⎛⎥⎦⎤⎢⎣⎡−−=⎟⎟⎠⎞⎜⎜⎝⎛⇒=⎟⎟⎠⎞⎜⎜⎝⎛⎟⎟⎠⎞⎜⎜⎝⎛⎥⎦⎤⎢⎣⎡−−−⎟⎟⎠⎞⎜⎜⎝⎛δδ=δ−δ2121211121a se u u k k k k F F 0u u k k k k F F u u W WBoundary conditionsWith the finite element method the application of displacement constraint boundary conditions is performed after the equations are assembled. It is an interest to examine the implications of applying and not applying the displacement boundary constraints prior to applying the principle of virtual work. Consider then the single spring element above but fixed at node 1, i.e. 0u 1=. Ignoring the constraint initially gives()212se u u k 21W −=, ()()1212se u u u u k W δ−δ−=δ and 2211a u F u F W δ+δ=δ.The principle of virtual work gives 2211ku ku ku F −=−= and 2212ku ku ku F =+−=, on applicationof the constraint. Note that is the force required at node 1 to prevent the node moving and is the reaction force.21ku F −=21ku F =−Applying the constraint straightaway gives 22se ku 21W =, 22se u ku W δ=δ and 22a u F W δ=δ. The principle of virtual work gives with no information about the reaction force at node 1.22ku F =Exam Standard Question:The spring-mass system depicted in the Figure consists of three massless springs, which are attached to fixed boundaries by means of pin-joints at nodes 1, 3 and 5. The springs are connected to a rigid bar by means of pin-joints at nodes 2 and 4. The rigid bar is free to rotate about pivot A. Nodes 2 and 4 are distances and below pivot A, respectively. Each spring has the same stiffness k. Node 2 is subjected to an external horizontal force F 2/l 4/l 2. All deflections can be assumed to be small.(i) Write expressions for the extension of each spring in terms of the displacement of node 2.(ii) In terms of the degrees of freedom at node 2, write expressions for the total strain energy W of the spring-mass system. In addition, specify the variation in work done se a W δ resulting from the application of the force.2F (iii) Use Use the principle of virtual work to find a relationship between the magnitude of and the horizontal components of displacement at node 2.2F (iv) Use the principle of virtual work to show that the net vertical force imposed by the springs on the rigid-bar at node 2 is zero.Solution:(i) Directional vectors for springs are: 2112e 21e 23e +=, 2132e 21e 23e +−= and 145e e =. Extensions for bottom springs are: 221212u 23u e =⋅=δ, 223232u 23u e −=⋅=δ.Note that 2u u 24=, so 2u245−=δ.(ii)()2222222245232212se ku 87u 212323k 21k 21W =⎟⎟⎠⎞⎜⎜⎝⎛⎟⎠⎞⎜⎝⎛+⎟⎟⎠⎞⎜⎜⎝⎛−+⎟⎟⎠⎞⎜⎜⎝⎛=δ+δ+δ=, 222u F W δ=δ(iii) 2222a 22se ku 47F u F W u ku 47W =⇒δ=δ=δ=δ(iv) Need additional displacement degree of freedom at node 2. Let 22122e v e u u += and note that2221212v 21u 23u e +=⋅=δ and 2223232v 21u 23u e +−=⋅=δ.()⎟⎟⎠⎞⎜⎜⎝⎛⎟⎟⎠⎞⎜⎜⎝⎛+−+⎟⎟⎠⎞⎜⎜⎝⎛+=δ+δ=222222232212se v 21u 23v 21u 23k 21k 21W ⎟⎟⎠⎞⎜⎜⎝⎛⎟⎟⎠⎞⎜⎜⎝⎛δ+δ−⎟⎟⎠⎞⎜⎜⎝⎛+−+⎟⎟⎠⎞⎜⎜⎝⎛δ+δ⎟⎟⎠⎞⎜⎜⎝⎛+=δ22222222se v 21u 23v 21u 23v 21u 23v 21u 23k W Setting and gives0v 2=0u 2=δ2vert 222222se v F v 0v 21u 23v 21u 23k W δ=δ=⎟⎟⎠⎞⎜⎜⎝⎛⎟⎠⎞⎜⎝⎛δ⎟⎟⎠⎞⎜⎜⎝⎛−+⎟⎠⎞⎜⎝⎛δ⎟⎟⎠⎞⎜⎜⎝⎛=δ hence . 0F vert 2=Method of Minimum PotentialConsider the expression,()()F u u u TT 21212121c se K 21F F u u u u k k k k u u 21W W P −=⎟⎟⎠⎞⎜⎜⎝⎛−⎟⎟⎠⎞⎜⎜⎝⎛⎥⎦⎤⎢⎣⎡−−=−=where W F and can be considered as a work term with independent of . u F u c =+1122F i u iThe approach of minimising P is known as the method of minimum potential .Note that,()()u F 0F -u u =F u u u +u u K K K K 21W W P T T T T c se =⇒=δδ−δδ=δ−δ=δwhere use has been made of the fact that δδu u =u u T TK K as a result of K 's symmetry.It is useful at this stage to consider the minimisation of an arbitrary functional ()u P where()()3T T O H 21P P u u u u u δ+δδ+∇δ=δand the gradient ∇=P P u i i ∂∂, and the Hessian matrix coefficients H P u u ij i j=∂∂∂2.A stationary point requires that ∇=, i.e.P 0∂∂Pu i=0.Moreover, a minimum point requires that δδu u TH >0 for all δu ≠0 and matrices that possess this property are known as positive definite .Setting P W W K se c T=−=−12u u u F T provides ∇=−=P K u F 0 and H K =.It is a simple matter to check that with u 10= (to prevent rigid body movement) that K is positive definite and this is a property commonly associated with FE stiffness matrices.Exam Standard Question:The spring system depicted in the Figure consists of four massless unstretched springs, which are attached to fixed boundaries by means of pin-joints at nodes 1 to 4. The springs are connected to a slider at node 5. Theslider is constrained to move in a frictionless channel whose axis is to the horizontal. Each spring has the same stiffness k. The slider is subjected to an external force F 0453 whose direction is along the axis of the frictionless channel.(i)The deflection of node 5 can be represented by the vector 25155v u e e u +=, where and areunit orthogonal vectors which are shown in the Figure. Write the components of deflection and in terms of , where is the magnitude of , i.e. e 1e 25u 5v 5U 5U 5u 25U 5u =. Show that the extensions of eachspring, in terms of , are: 5U ()22/31U 515+=δ, ()22/31U 525−=δ, and2/U 54535−=δ−=δ.(ii) In terms of k and write expressions for the total strain energy W of the spring-mass system. Inaddition, specify the variation in work done 5U se a W δ resulting from the application of the force . 5F (iii) Use the principal of virtual work to find a relationship between the magnitude of and thedisplacement at node 5.5F 5U (iv) Use the principal of virtual work to determine an expression for the force imposed by the frictionless channel on the slider.(v)Form a potential energy function for the spring system. Assume here that nodes 1, 3 and 4 are fixed and node 5 is restricted to move in the channel. Use this function to determine the reaction force at node 2.Solution:(i) Directional vectors for springs and channel are: ()2115e e 321e +=, ()2125e e 321e +−=, 135e e −=, 45e e = and (21c 5e e 21e +=). Deflection c 555e U u =, so 2U v u 555==. Extensions springs are: ()3122U u e 551515+=⋅=δ, ()3122U u e 552525−=⋅=δ, 2Uu e 553535−=⋅=δand 2Uu e 554545=⋅=δ(ii)()()()252522245235225215se kU U 83131k 8121k 21W =⎟⎠⎞⎜⎝⎛+−++=δ+δ+δ+δ=, 55a U F W δ=δ(iii)5555a 55se kU 2F U F W U kU 2W =⇒δ=δ=δ=δ(iv) Need additional displacement degree of freedom at node 3. A unit vector perpendicular to the channel is(21p 5e e 21e +−=) and let p 55c 555e V e U u += and note that()()3122V3122U u e 5551515−++=⋅=δ and ()()3122V3122U u e 5552525++−=⋅=δ, 2V 2U u e 5553535+−=⋅=δ and 2V 2U u e 5554545−=⋅=δ()()()()()()()()⎟⎠⎞⎜⎝⎛−+++−+−++=δ+δ+δ+δ=255255255245235225215se V U 831V 31U 31V 31U k 8121k 21W ()()()()()555se V 0V 831313131kU 81W δ=δ−+−+−+=δ, where variation is onlyconsidered and is set to zero. Principle of virtual work .5V δ3V 0F V F V 0W p 55p 55se =⇒δ=δ=δ(v)()3122U u e 551515+=⋅=δ, ()()5552525V 3122Uu u e −−=−⋅=δ, where 2522e V u =. ()()()223333223232233245235225215V F U F U V 3122U 3122U k 21V F U F k 21P −−⎟⎟⎠⎞⎜⎜⎝⎛+⎥⎦⎤⎢⎣⎡−−+⎥⎦⎤⎢⎣⎡+=−−δ+δ+δ+δ=and ()0F V 3122U k V P 2232=−⎥⎦⎤⎢⎣⎡−−−=∂∂, which on setting 0V 2= gives ()⎥⎦⎤⎢⎣⎡−−=3122U k F 32.The reaction is .2F −System AssemblyConsider the following three-spring system 2F 21u 1F 1u 2kF 3F 4u 3u 4k 1k2334()()()234322322121se u u k 21u u k 21u u k 21W −+−+−=,()()()()()()343432323212121se u u u u k u u u u k u u u u k W δ−δ−+δ−δ−+δ−δ−=δ,44332211a u F u F u F u F W δ+δ+δ+δ=δ,and δδ implies that,W W se a −=0u F K u u u u k k 0k k k k 00k k k k 00k k F FF F 43213333222211114321=⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎝⎛⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡−−+−−+−−=⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎝⎛=where again it is apparent that K is symmetric but also it is banded, i.e. the non-zero coefficients are located around the principal diagonal. This is a property commonly associated with assembled FE stiffness matrices and depends on node connectivity. Note also that the summation of coefficients in individual rows or columns gives zero. The matrix is singular and 0K det =.Note that element stiffness matrices are: , and where on examination of K it is apparent how these are assembled to form K .⎥⎦⎤⎢⎣⎡−−1111k k k k ⎥⎦⎤⎢⎣⎡−−2222k k k k ⎥⎦⎤⎢⎣⎡−−3333k k k kIf a boundary constraint is imposed then row one is removed to give:0u 1=u F K u u u k k 0k k k k 0k k k F F F 432333322221432=⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡−−+−−+=⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛=. If however a boundary constraint (say) is imposed then row one is again removed but a somewhatdifferent answer is obtained: 1u 1=u F K u u u k k 0k k k k 0k k k F F k F 4323333222214312=⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡−−+−−+=⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛+=)Direct FormulationIt is possible to formulate the stiffness matrix directly by moving one node and keeping the others fixed and noting the reactions.The above system can be solved for u , once possible rigid body motion is prevented, by setting u (say) to give 10=⇒=⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡−−+−−+=⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛=u F K u u u k k 0k k k k 0k k k F F F 432333322221432⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡−−+−−+=⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛−4321333322221432F F F k k 0k k k k 0k k k u u uThe inverse stiffness matrix, K −1, is known as the flexibility matrix and, for this example at least, can be assembled directly by noting the system response to prescribed forces.In practice K −1is never calculated and the system K u F = is solved using a modern numerical linear system solver.It is a simple matter to confirm thatu u K 21u u u u k k 0k k k k 00k k k k 00k k u u u u 21W T 4321333322221111T4321se =⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎝⎛⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡−−+−−+−−⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎝⎛= with F u T4321T4321a F F F F u u u u W δ=⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎝⎛⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎝⎛δδδδ=δThus,()u F F u u K 0K W W Ta se =⇒=−δ=δ−δExample:k1F 2u23k2u3F321With use a direct method to find the assembled stiffness and flexibility matrices.0u 1=Solution:The equations of interest are of the form: 3232222u k u k F += and 3332323u k u k F +=.Consider and equilibrium at nodes 2 and 3. At node 2, 0u 3=()2212u k k F += and at node 3,.223u k F −=Consider and equilibrium at nodes 2 and 3. At node 2, 0u 2=322u k F −= and at node 3, . 323u k F =Thus: , , 2122k k k +=223k k −=232k k −= and 233k k =.For flexibility the equations of interest are of the form: 3232222F c F c u += and . 3332323F c F c u +=Consider and equilibrium at nodes 2 and 3. At node 2, 0F 3=122k F u = and at node 3,1223k F u u ==.Consider and equilibrium at nodes 2 and 3. At node 2, 0F 2=122k F u = and at node 3,()2133k 1k 1F u +=.Thus: 122k 1c =, 123k 1c =, 132k 1c = and 2133k 1k 1c +=.Can check that ⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡−−+1001k 1k 1k 1k 1k 1k k k k k 2111122221 as required,It should be noted that the direct determination requires boundary constraints to be applied to ensure that the flexibility matrix exists, which requires the stiffness to be non-singular. However, the stiffness matrix always exists, so boundary conditions need not be applied prior to constructing the stiffness matrix with the direct approach.Large deformation theory for spring elementsThus far small deflection theory has been applied where the strains are measured using the Cauchy strainxu11∂∂=ε. A conjugate stress can be obtained by differentiating with respect the expression for strain energy density (energy per unit volume) 11ε211E 21ε=ω, i.e. 111111E ε=ε∂ω∂=σ, where E is Young’s Modulusand is the Cauchy stress (sometimes referred to as the Euler stress). 11σIn the case of large deformation theory we will restrict our attention to hyperelastic materials which are materials that possess an expression for strain energy density Ω (say) that is analytical in strain.The strain used in large deformation theory is Green’s strain (see Appendix II) which for a uniformly loadeduniaxial bar is 211x u 21x u E ⎟⎠⎞⎜⎝⎛∂∂+∂∂=.An expression for strain energy density (energy per unit volume) 211EE 21=Ω and the derived stress is 111111EE E S =∂Ω∂=, where E is Young’s Modulus and is known as the 211S nd Piola-Kirchoff stress . 2F21u1F1u2kBar subject to longitudinal deformationConsider a bar of length L and cross sectional area A represented by a spring element and subject to nodal forces and . 1F 2FThe strain energy is∫∫∫∫⎥⎥⎦⎤⎢⎢⎣⎡⎟⎠⎞⎜⎝⎛∂∂+∂∂==Ω=Ω=212121x x 22x x 211x x V se dx x u 21x u EA 21dx E EA 21dx A dV WConsider further a linear displacement field of the form ()21u L x u L x L x u ⎟⎠⎞⎜⎝⎛+⎟⎠⎞⎜⎝⎛−= and note thatL u u xu 12−=∂∂. ()()221212x x 221212se u u L 21u u L EA 21dx L u u 21L u u EA 21W 21⎥⎦⎤⎢⎣⎡−+−=⎥⎥⎦⎤⎢⎢⎣⎡⎟⎠⎞⎜⎝⎛−+−=∫ ()()()⎥⎦⎤⎢⎣⎡−+−+−=4122312212se u u L 41u u L 1u u k 21W()()()(12312221212se u u u u L 21u u L 23u u k W δ−δ⎥⎦⎤⎢⎣⎡−+−+−=δ) and 2211a u F u F W δ+δ=δ.The principle of virtual work gives()()()⎥⎦⎤⎢⎣⎡−+−+−−=3122212121u u L 21u u L 23u u k F and()()(⎥⎦⎤⎢⎣⎡−+−+−=3122212122u u L 21u u L 23u u k F ), represented in matrix form as()()()()()⎟⎟⎠⎞⎜⎜⎝⎛⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡−+−−−−−−−+−+⎥⎦⎤⎢⎣⎡−−=⎟⎟⎠⎞⎜⎜⎝⎛21121212121221u u L 3u u 1L 3u u 1L 3u u 1L 3u u 1L 2u u k 3k k k k F Fwhich is of the form[]u F G L K K += where is called the geometrical stiffness matrix and is the usual linear stiffnessmatrix. G K L KA common approximation used, depending on the magnitude of L /u u 12−, is⎟⎟⎠⎞⎜⎜⎝⎛⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡−−+⎥⎦⎤⎢⎣⎡−−=⎟⎟⎠⎞⎜⎜⎝⎛2121u u 1111L 2P 3k k k k F F where ()12u u k P −=.The fact that is non-linear (even in its approximate form) means that iterative solution procedures are required to be employed to determine the unknown displacements. G KNote that the approximate form is arrived at using the following strain energy expression()()⎥⎦⎤⎢⎣⎡−+−=312212se u u L 1u u k 21WExample:The strain energies for the springs in the above system (fixed at node 1) are k 1 F 2u 23k 2u3F321⎥⎦⎤⎢⎣⎡+=1322211seL u u k 21W and ()()⎥⎦⎤⎢⎣⎡−+−=323222322se u u L 1u u k 21WUse the principle of virtual work to obtain the assembled linear and geometrical stiffness matrices.()()()3322a 2322322322122212se1sese u F u F W u u u u L 23u u k u L 2u 3u k W W W δ+δ=δ=δ−δ⎥⎦⎤⎢⎣⎡−+−+δ⎥⎦⎤⎢⎣⎡+=δ+δ=δThus ()(⎥⎦⎤⎢⎣⎡−+−−⎥⎦⎤⎢⎣⎡+=2232232122212u u L 23u u k L 2u 3u k F ) and ()()⎥⎦⎤⎢⎣⎡−+−=22322323u u L 23u u k F⎟⎟⎠⎞⎜⎜⎝⎛⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡αα−α−α+α+⎥⎦⎤⎢⎣⎡−−+=⎟⎟⎠⎞⎜⎜⎝⎛32222212222132u u k k k k k F F where 1211L 2u k 3=α and ()23222u u L 2k 3−=α.Note that the element stiffness matrices are[][]111k K α+= and ⎥⎦⎤⎢⎣⎡αα−α−α+⎥⎦⎤⎢⎣⎡−−=222222222k k k k Kand it is evident how these should be assembled to form the assembled linear and geometrical stiffness matrices.2v21u 1v1u 2kxBar subject to longitudinal and lateral deflectionConsider a bar of length L and cross sectional area A represented by a spring element and subject to longitudinal and lateral displacements u and v, respectively.The normal strain is 2211x v 21x u 21x u E ⎟⎠⎞⎜⎝⎛∂∂+⎟⎠⎞⎜⎝⎛∂∂+∂∂= and the associated strain energy∫∫∫∫⎥⎥⎦⎤⎢⎢⎣⎡⎟⎠⎞⎜⎝⎛∂∂+⎟⎠⎞⎜⎝⎛∂∂+∂∂==Ω=Ω=212121x x 22x x 211x x V se dx x v 21x u 21x u EA 21dx E EA 21dx A dV W ∫⎥⎥⎦⎤⎢⎢⎣⎡⎟⎠⎞⎜⎝⎛∂∂⎟⎠⎞⎜⎝⎛∂∂+⎟⎠⎞⎜⎝⎛∂∂+⎟⎠⎞⎜⎝⎛∂∂≈21x x 232se dx x v x u x u x u EA 21WConsider further a linear displacement field of the form ()21u L x u L x L x u ⎟⎠⎞⎜⎝⎛+⎟⎠⎞⎜⎝⎛−= and()21v L x v L x L x v ⎟⎠⎞⎜⎝⎛+⎟⎠⎞⎜⎝⎛−=, and note thatL u u x u 12−=∂∂ and L v v x v 12−=∂∂. ()()()()⎥⎦⎤⎢⎣⎡−−+−+−=L v v u u L u u u u L EA 21W 21212312212se()()()()()()()1212121221221212se v v L v v u u k u u L 2v v L 2u u 3u u k W δ−δ⎦⎤⎢⎣⎡−−+δ−δ⎥⎦⎤⎢⎣⎡−+−+−=δ2v 22h 21v 11h 1a v F u F v F u F W δ+δ+δ+δ=δ and the principle of virtual work gives()()()⎥⎦⎤⎢⎣⎡−+−+−−=L 2v v L 2u u 3u u k F 21221212h1and ()()⎥⎦⎤⎢⎣⎡−−−=L v v u u k F 1212v1 ()()()⎥⎦⎤⎢⎣⎡−+−+−=L 2v v L 2u u 3u u k F 21221212h2and ()()⎥⎦⎤⎢⎣⎡−−=L v v u u k F 1212v2()()⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎝⎛⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎝⎛⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡−−−−−+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡−−−+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡−−=⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎝⎛22111212v 2h 2v 1h 1v u v u 101005.105.1101005.105.1Lu u k 1010000010100000L2v v k 0000010100000101k F F F FExam Standard Question:The spring system depicted in the Figure consists of two massless springs of equal length , which are attached to fixed boundaries by means of pin-joints at nodes 1 and 2. The springs are connected to a slider atnode 3. The slider is constrained to move in a frictionless channel whose axis is 45 to the horizontal. Each spring has the same stiffness . The slider is subjected to an external force F 1L =0L /EA k =3 whose direction is along the axis of the frictionless channel.FigureAssume the springs have strain density ⎥⎥⎦⎤⎢⎢⎣⎡⎟⎠⎞⎜⎝⎛∂∂⎟⎠⎞⎜⎝⎛∂∂+⎟⎠⎞⎜⎝⎛∂∂+⎟⎠⎞⎜⎝⎛∂∂=Ω232x v x u x u x u E 21.(i) Write expressions for the longitudinal and lateral displacements for each spring at node 3 in terms of thedisplacement along the channel at node 3.(ii) In terms of displacement along the channel at node 3, write expressions for the total strain energy W of thespring-mass system. In addition, specify the variation in work done se a W δ resulting from the application of the force .3F (iii) Use the principle of virtual work to find a relationship between the magnitude of and the displacementalong the channel at node 3. 3FSolution:(i) Directional vectors for springs and channel are: ()2113e e 321e +=and ()2123e e 321e +−= and (21c 3e e 21e +=). Perpendicular vectors are: ()2113e 3e 21e +−=⊥and ()2123e 3e 21e +=⊥Deflection c 333e U u =, so 2U v u 333==.Longitudinal displacement: ()3122U u e 331313+=⋅=δ, ()3122U u e 332323−=⋅=δ.Lateral displacement: ()3122U u e 331313+−=⋅=δ⊥⊥, ()3122U u e 332323+=⋅=δ⊥⊥(ii) The strain energy density for element 1 is ⎥⎥⎦⎤⎢⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛δ⎟⎠⎞⎜⎝⎛δ+⎟⎠⎞⎜⎝⎛δ+⎟⎠⎞⎜⎝⎛δ=Ω⊥21313313213L L L L E 21 The strain energy density for element 2 is ⎥⎥⎦⎤⎢⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛δ⎟⎠⎞⎜⎝⎛δ+⎟⎠⎞⎜⎝⎛δ+⎟⎠⎞⎜⎝⎛δ=Ω⊥22323323223L L L L E 21 The total strain energy with substitution of 1L = gives()()()()[]()()()()[][]3322312232332322321313313213se U Uk 21k 21k 21W α+α=δδ+δ+δ+δδ+δ+δ=⊥⊥where and are constants determined on collecting up terms on substitution of and .1α2α231313,,δδδ⊥⊥δ2333a U F W δ=δ.(iii) The principle of virtual work gives⎥⎦⎤⎢⎣⎡α+α=⇒δ=δ=δ⎥⎦⎤⎢⎣⎡α+α=δ32133332323231se U 23kU F U F W U U 23U k WPin-jointed structuresThe example above is a pin-jointed structure. A reasonable good approximation reported in the literature for strain energy density, commonly used with pin-jointed structures, is⎥⎥⎦⎤⎢⎢⎣⎡⎟⎠⎞⎜⎝⎛∂∂⎟⎠⎞⎜⎝⎛∂∂+⎟⎠⎞⎜⎝⎛∂∂=Ω22x v x u x u E 21This arises from strain-energy approximation 211x v 21x u E ⎟⎠⎞⎜⎝⎛∂∂+∂∂=. Can be used when 22x v x u ⎟⎠⎞⎜⎝⎛∂∂<<⎟⎠⎞⎜⎝⎛∂∂.。

双向板等效均布荷载计算公式英文回答：Calculation of Equivalent Uniform Load for Two-Way Slabs.The design of two-way slabs requires the calculation of the equivalent uniform load (EUL) that represents the actual distributed loading on the slab. The EUL is used in the design of the slab reinforcement and can be calculated using one of several methods.Method 1: Strip Method.The strip method is a simplified approach that divides the slab into a series of parallel strips and assumes that the load is uniformly distributed along each strip. The EUL for each strip is then calculated as follows:w = (L/2) q.where:w is the EUL (in lb/ft or kN/m)。

L is the length of the strip (in ft or m)。

q is the uniformly distributed load (in lb/ft² orkN/m²)。

The EUL for the entire slab is then calculated by summing the EULs for all the strips.Method 2: Finite Element Method.The finite element method (FEM) is a more accurate method that uses computer software to analyze the behavior of the slab under various loading conditions. The FEM can account for the effects of boundary conditions, concentrated loads, and other factors that may affect the distribution of the load on the slab.The FEM is a complex method that requires specialized software and expertise to use. However, it can provide more accurate results than the strip method, especially for slabs with irregular shapes or complex loading conditions.Choice of Method.The choice of method for calculating the EUL depends on the complexity of the slab and the desired level of accuracy. The strip method is a simple and straightforward approach that is suitable for most regular slabs with uniformly distributed loads. The FEM is a more accurate method that should be used for slabs with irregular shapes or complex loading conditions.中文回答：双向板等效均布荷载计算公式。

自动化英语专业英语词汇表文章摘要：本文介绍了自动化英语专业的一些常用的英语词汇，包括自动化技术、控制理论、系统工程、人工智能、模糊逻辑等方面的专业术语。

本文按照字母顺序，将这些词汇分为26个表格，每个表格包含了以相应字母开头的词汇及其中文释义。

本文旨在帮助自动化专业的学习者和从业者掌握和使用这些专业英语词汇，提高他们的英语水平和专业素养。

A英文中文acceleration transducer加速度传感器acceptance testing验收测试accessibility可及性accumulated error累积误差AC-DC-AC frequency converter交-直-交变频器AC (alternating current) electric drive交流电子传动active attitude stabilization主动姿态稳定actuator驱动器，执行机构adaline线性适应元adaptation layer适应层adaptive telemeter system适应遥测系统adjoint operator伴随算子admissible error容许误差aggregation matrix集结矩阵AHP (analytic hierarchy process)层次分析法amplifying element放大环节analog-digital conversion模数转换annunciator信号器antenna pointing control天线指向控制anti-integral windup抗积分饱卷aperiodic decomposition非周期分解a posteriori estimate后验估计approximate reasoning近似推理a priori estimate先验估计articulated robot关节型机器人assignment problem配置问题，分配问题associative memory model联想记忆模型associatron联想机asymptotic stability渐进稳定性attained pose drift实际位姿漂移B英文中文attitude acquisition姿态捕获AOCS (attritude and orbit control system)姿态轨道控制系统attitude angular velocity姿态角速度attitude disturbance姿态扰动attitude maneuver姿态机动attractor吸引子augment ability可扩充性augmented system增广系统automatic manual station自动-手动操作器automaton自动机autonomous system自治系统backlash characteristics间隙特性base coordinate system基座坐标系Bayes classifier贝叶斯分类器bearing alignment方位对准bellows pressure gauge波纹管压力表benefit-cost analysis收益成本分析bilinear system双线性系统biocybernetics生物控制论biological feedback system生物反馈系统C英文中文calibration校准，定标canonical form标准形式canonical realization标准实现capacity coefficient容量系数cascade control级联控制causal system因果系统cell单元，元胞cellular automaton元胞自动机central processing unit (CPU)中央处理器certainty factor确信因子characteristic equation特征方程characteristic function特征函数characteristic polynomial特征多项式characteristic root特征根英文中文charge-coupled device (CCD)电荷耦合器件chaotic system混沌系统check valve单向阀，止回阀chattering phenomenon颤振现象closed-loop control system闭环控制系统closed-loop gain闭环增益cluster analysis聚类分析coefficient of variation变异系数cogging torque齿槽转矩，卡齿转矩cognitive map认知图，认知地图coherency matrix相干矩阵collocation method配点法，配置法combinatorial optimization problem组合优化问题common mode rejection ratio (CMRR)共模抑制比，共模抑制率commutation circuit换相电路，换向电路commutator motor换向电动机D英文中文damping coefficient阻尼系数damping ratio阻尼比data acquisition system (DAS)数据采集系统data fusion数据融合dead zone死区decision analysis决策分析decision feedback equalizer (DFE)决策反馈均衡器decision making决策，决策制定decision support system (DSS)决策支持系统decision table决策表decision tree决策树decentralized control system分散控制系统decoupling control解耦控制defuzzification去模糊化，反模糊化delay element延时环节，滞后环节delta robot德尔塔机器人demodulation解调，检波density function密度函数，概率密度函数derivative action微分作用，微分动作design matrix设计矩阵E英文中文eigenvalue特征值，本征值eigenvector特征向量，本征向量elastic element弹性环节electric drive电子传动electric potential电势electro-hydraulic servo system电液伺服系统electro-mechanical coupling system电机耦合系统electro-pneumatic servo system电气伺服系统electronic governor电子调速器encoder编码器，编码装置end effector末端执行器，末端效应器entropy熵equivalent circuit等效电路error analysis误差分析error bound误差界，误差限error signal误差信号estimation theory估计理论Euclidean distance欧几里得距离，欧氏距离Euler angle欧拉角Euler equation欧拉方程F英文中文factor analysis因子分析factorization method因子法，因式分解法feedback反馈，反馈作用feedback control反馈控制feedback linearization反馈线性化feedforward前馈，前馈作用feedforward control前馈控制field effect transistor (FET)场效应晶体管filter滤波器，滤波环节finite automaton有限自动机finite difference method有限差分法finite element method (FEM)有限元法finite impulse response (FIR) filter有限冲激响应滤波器first-order system一阶系统fixed-point iteration method不动点迭代法flag register标志寄存器flip-flop circuit触发器电路floating-point number浮点数flow chart流程图，流程表fluid power system流体动力系统G英文中文gain增益gain margin增益裕度Galerkin method伽辽金法game theory博弈论Gauss elimination method高斯消元法Gauss-Jordan method高斯-约当法Gauss-Markov process高斯-马尔可夫过程Gauss-Seidel iteration method高斯-赛德尔迭代法genetic algorithm (GA)遗传算法gradient method梯度法，梯度下降法graph theory图论gravity gradient stabilization重力梯度稳定gray code格雷码，反向码gray level灰度，灰阶grid search method网格搜索法ground station地面站，地面控制站guidance system制导系统，导航系统gyroscope陀螺仪，陀螺仪器H英文中文H∞ control H无穷控制Hamiltonian function哈密顿函数harmonic analysis谐波分析harmonic oscillator谐振子，谐振环节Hartley transform哈特利变换Hebb learning rule赫布学习规则Heisenberg uncertainty principle海森堡不确定性原理hidden layer隐层，隐含层hidden Markov model (HMM)隐马尔可夫模型hierarchical control system分层控制系统high-pass filter高通滤波器Hilbert transform希尔伯特变换Hopfield network霍普菲尔德网络hysteresis滞后，迟滞，磁滞I英文中文identification识别，辨识identity matrix单位矩阵，恒等矩阵image processing图像处理impulse response冲激响应impulse response function冲激响应函数inadmissible control不可接受控制incremental encoder增量式编码器indefinite integral不定积分index of controllability可控性指标index of observability可观测性指标induction motor感应电动机inertial navigation system (INS)惯性导航系统inference engine推理引擎，推理机inference rule推理规则infinite impulse response (IIR) filter无限冲激响应滤波器information entropy信息熵information theory信息论input-output linearization输入输出线性化input-output model输入输出模型input-output stability输入输出稳定性J英文中文Jacobian matrix雅可比矩阵jerk加加速度，冲击joint coordinate system关节坐标系joint space关节空间Joule's law焦耳定律jump resonance跳跃共振K英文中文Kalman filter卡尔曼滤波器Karhunen-Loeve transform卡尔胡南-洛维变换kernel function核函数，核心函数kinematic chain运动链，运动链条kinematic equation运动方程，运动学方程kinematic pair运动副，运动对kinematics运动学kinetic energy动能L英文中文Lagrange equation拉格朗日方程Lagrange multiplier拉格朗日乘子Laplace transform拉普拉斯变换Laplacian operator拉普拉斯算子laser激光，激光器latent root潜根，隐根latent vector潜向量，隐向量learning rate学习率，学习速度least squares method最小二乘法Lebesgue integral勒贝格积分Legendre polynomial勒让德多项式Lennard-Jones potential莱纳德-琼斯势level set method水平集方法Liapunov equation李雅普诺夫方程Liapunov function李雅普诺夫函数Liapunov stability李雅普诺夫稳定性limit cycle极限环，极限圈linear programming线性规划linear quadratic regulator (LQR)线性二次型调节器linear system线性系统M英文中文machine learning机器学习machine vision机器视觉magnetic circuit磁路，磁电路英文中文magnetic flux磁通量magnetic levitation磁悬浮magnetization curve磁化曲线magnetoresistance磁阻，磁阻效应manipulability可操作性，可操纵性manipulator操纵器，机械手Markov chain马尔可夫链Markov decision process (MDP)马尔可夫决策过程Markov property马尔可夫性质mass matrix质量矩阵master-slave control system主从控制系统matrix inversion lemma矩阵求逆引理maximum likelihood estimation (MLE)最大似然估计mean square error (MSE)均方误差measurement noise测量噪声，观测噪声mechanical impedance机械阻抗membership function隶属函数N英文中文natural frequency固有频率，自然频率natural language processing (NLP)自然语言处理navigation导航，航行negative feedback负反馈，负反馈作用neural network神经网络neuron神经元，神经细胞Newton method牛顿法，牛顿迭代法Newton-Raphson method牛顿-拉夫逊法noise噪声，噪音nonlinear programming非线性规划nonlinear system非线性系统norm范数，模，标准normal distribution正态分布，高斯分布notch filter凹槽滤波器，陷波滤波器null space零空间，核空间O英文中文observability可观测性英文中文observer观测器，观察器optimal control最优控制optimal estimation最优估计optimal filter最优滤波器optimization优化，最优化orthogonal matrix正交矩阵oscillation振荡，振动output feedback输出反馈output regulation输出调节P英文中文parallel connection并联，并联连接parameter estimation参数估计parity bit奇偶校验位partial differential equation (PDE)偏微分方程passive attitude stabilization被动姿态稳定pattern recognition模式识别PD (proportional-derivative) control比例-微分控制peak value峰值，峰值幅度perceptron感知器，感知机performance index性能指标，性能函数period周期，周期时间periodic signal周期信号phase angle相角，相位角phase margin相位裕度phase plane analysis相平面分析phase portrait相轨迹，相图像PID (proportional-integral-derivative) control比例-积分-微分控制piezoelectric effect压电效应pitch angle俯仰角pixel像素，像元Q英文中文quadratic programming二次规划quantization量化，量子化quantum computer量子计算机quantum control量子控制英文中文queueing theory排队论quiescent point静态工作点，静止点R英文中文radial basis function (RBF) network径向基函数网络radiation pressure辐射压random variable随机变量random walk随机游走range范围，区间，距离rank秩，等级rate of change变化率，变化速率rational function有理函数Rayleigh quotient瑞利商real-time control system实时控制系统recursive algorithm递归算法recursive estimation递归估计reference input参考输入，期望输入reference model参考模型，期望模型reinforcement learning强化学习relay control system继电器控制系统reliability可靠性，可信度remote control system遥控系统，远程控制系统residual error残差误差，残余误差resonance frequency共振频率S英文中文sampling采样，取样sampling frequency采样频率sampling theorem采样定理saturation饱和，饱和度scalar product标量积，点积scaling factor缩放因子，比例系数Schmitt trigger施密特触发器Schur complement舒尔补second-order system二阶系统self-learning自学习，自我学习self-organizing map (SOM)自组织映射sensitivity灵敏度，敏感性sensitivity analysis灵敏度分析，敏感性分析sensor传感器，感应器sensor fusion传感器融合servo amplifier伺服放大器servo motor伺服电机，伺服马达servo valve伺服阀，伺服阀门set point设定值，给定值settling time定常时间，稳定时间T英文中文tabu search禁忌搜索，禁忌表搜索Taylor series泰勒级数，泰勒展开式teleoperation遥操作，远程操作temperature sensor温度传感器terminal终端，端子testability可测试性，可检测性thermal noise热噪声，热噪音thermocouple热电偶，热偶threshold阈值，门槛time constant时间常数time delay时延，延时time domain时域time-invariant system时不变系统time-optimal control时间最优控制time series analysis时间序列分析toggle switch拨动开关，切换开关tolerance analysis公差分析torque sensor扭矩传感器transfer function传递函数，迁移函数transient response瞬态响应U英文中文uncertainty不确定性，不确定度underdamped system欠阻尼系统undershoot低于量，低于值unit impulse function单位冲激函数unit step function单位阶跃函数unstable equilibrium point不稳定平衡点unsupervised learning无监督学习upper bound上界，上限utility function效用函数，效益函数V英文中文variable structure control变结构控制variance方差，变异vector product向量积，叉积velocity sensor速度传感器verification验证，校验virtual reality虚拟现实viscosity粘度，黏度vision sensor视觉传感器voltage电压，电位差voltage-controlled oscillator (VCO)电压控制振荡器W英文中文wavelet transform小波变换weighting function加权函数Wiener filter维纳滤波器Wiener process维纳过程work envelope工作空间，工作范围worst-case analysis最坏情况分析X英文中文XOR (exclusive OR) gate异或门，异或逻辑门Y英文中文yaw angle偏航角Z英文中文Z transform Z变换zero-order hold (ZOH)零阶保持器zero-order system零阶系统zero-pole cancellation零极点抵消。

基于MPI+FreeFem++的有限元并行计算摘要：有限元方法是一种灵活而高效的数值求解偏微分方程的计算方法，是工程分析和计算中不可缺少的重要工具之一。

在计算机技术的快速发展使得并行机的价格日益下降的今天，并行有限元计算方法受到了学术界和工程界的普遍关注。

讨论了基于MPI+FreeFem++的有限元并行计算环境的构建，阐述了在该环境下有限元并行程序的编写、编译及运行等过程，并通过具体编程实例，说明了MPI+FreeFem++环境下的有限元并行编程的简单和高效。

关键词：有限元方法；并行计算；MPI；FreeFem++0 引言有限元方法是20世纪50年代伴随电子计算机的诞生，在计算数学和计算工程领域里诞生的一种高效而灵活的计算方法，它将古典变分法与分片多项式插值相结合，易于处理复杂的边值问题，具有有限差分法无可比拟的优越性，广泛应用于求解热传导、电磁场、流体力学等相关问题，已成为当今工程分析和计算中不可缺少的最重要的工具之一。

有限元方法的“化整为零、积零为整”的基本思想与并行处理技术的基本原则“分而治之”基本一致，因而具有高度的内在并行性。

在计算机技术快速发展使得并行机的价格日益下降的今天，有限元并行计算引起了学术界和工程界的普遍关注，吸引了众多科研与工程技术人员。

但要实现有限元并行编程，并不是一件容易的事，特别是对于复杂区域问题，若从网格生成、任务的划分、单元刚度矩阵的计算、总刚度矩阵的组装，到有限元方程组的求解以及后处理，都需要程序员自己编写代码的话，将是一件十分繁琐的事情。

本文探讨了构建基于MPI+FreeFem++的有限元并行计算环境，在该环境下，程序员可避免冗长代码的编写，进而轻松、快速、高效地实现复杂问题的有限元并行计算。

1 FreeFem++简介FreeFem++ 是一款免费的偏微分方程有限元计算软件，它集成网格生成器、线性方程组的求解器、后处理及计算结果可视化于一体，能快速而高效地实现复杂区域问题的有限元数值计算。

有限元仿真的英语Finite Element SimulationThe field of engineering has seen a remarkable evolution in recent decades, with the advent of advanced computational tools and techniques that have revolutionized the way we approach design, analysis, and problem-solving. One such powerful tool is the finite element method (FEM), a numerical technique that has become an indispensable part of the modern engineer's toolkit.The finite element method is a powerful computational tool that allows for the simulation and analysis of complex physical systems, ranging from structural mechanics and fluid dynamics to heat transfer and electromagnetic phenomena. At its core, the finite element method involves discretizing a continuous domain into a finite number of smaller, interconnected elements, each with its own set of properties and governing equations. By solving these equations numerically, the finite element method can provide detailed insights into the behavior of the system, enabling engineers to make informed decisions and optimize their designs.One of the key advantages of the finite element method is its abilityto handle complex geometries and boundary conditions. Traditional analytical methods often struggle with intricate shapes and boundary conditions, but the finite element method can easily accommodate these complexities by breaking down the domain into smaller, manageable elements. This flexibility allows engineers to model real-world systems with a high degree of accuracy, leading to more reliable and efficient designs.Another important aspect of the finite element method is its versatility. The technique can be applied to a wide range of engineering disciplines, from structural analysis and fluid dynamics to heat transfer and electromagnetic field simulations. This versatility has made the finite element method an indispensable tool in the arsenal of modern engineers, allowing them to tackle a diverse array of problems with a single computational framework.The power of the finite element method lies in its ability to provide detailed, quantitative insights into the behavior of complex systems. By discretizing the domain and solving the governing equations numerically, the finite element method can generate comprehensive data on stresses, strains, temperatures, fluid flow patterns, and other critical parameters. This information is invaluable for engineers, as it allows them to identify potential failure points, optimize designs, and make informed decisions that lead to more reliable and efficient products.The implementation of the finite element method, however, is not without its challenges. The process of discretizing the domain, selecting appropriate element types, and defining boundary conditions can be complex and time-consuming. Additionally, the accuracy of the finite element analysis is heavily dependent on the quality of the input data, the selection of appropriate material models, and the proper interpretation of the results.To address these challenges, researchers and software developers have invested significant effort in improving the finite element method and developing user-friendly software tools. Modern finite element analysis (FEA) software packages, such as ANSYS, ABAQUS, and COMSOL, provide intuitive graphical user interfaces, advanced meshing algorithms, and powerful post-processing capabilities, making the finite element method more accessible to engineers of all levels of expertise.Furthermore, the ongoing advancements in computational power and parallel processing have enabled the finite element method to tackle increasingly complex problems, pushing the boundaries of what was previously possible. High-performance computing (HPC) clusters and cloud-based computing resources have made it possible to perform large-scale, multi-physics simulations, allowing engineers to gain deeper insights into the behavior of their designs.As the engineering field continues to evolve, the finite element method is poised to play an even more pivotal role in the design, analysis, and optimization of complex systems. With its ability to handle a wide range of physical phenomena, the finite element method has become an indispensable tool in the modern engineer's toolkit, enabling them to push the boundaries of innovation and create products that are more reliable, efficient, and sustainable.In conclusion, the finite element method is a powerful computational tool that has transformed the field of engineering. By discretizing complex domains and solving the governing equations numerically, the finite element method provides engineers with detailed insights into the behavior of their designs, allowing them to make informed decisions and optimize their products. As the field of engineering continues to evolve, the finite element method will undoubtedly remain a crucial component of the modern engineer's arsenal, driving innovation and shaping the future of technological advancement.。