CN114818422A - 一种弹性波数值仿真分析方法与系统 - Google Patents

一种弹性波数值仿真分析方法与系统 Download PDF

Info

Publication number
CN114818422A
CN114818422A CN202210412326.XA CN202210412326A CN114818422A CN 114818422 A CN114818422 A CN 114818422A CN 202210412326 A CN202210412326 A CN 202210412326A CN 114818422 A CN114818422 A CN 114818422A
Authority
CN
China
Prior art keywords
grid
boundary
module
grid points
points
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202210412326.XA
Other languages
English (en)
Other versions
CN114818422B (zh
Inventor
富明慧
林美鸿
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Sun Yat Sen University
Original Assignee
Sun Yat Sen University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Sun Yat Sen University filed Critical Sun Yat Sen University
Priority to CN202210412326.XA priority Critical patent/CN114818422B/zh
Publication of CN114818422A publication Critical patent/CN114818422A/zh
Application granted granted Critical
Publication of CN114818422B publication Critical patent/CN114818422B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/10Geometric CAD
    • G06F30/17Mechanical parametric or variational design
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/10Numerical modelling
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14Force analysis or force optimisation, e.g. static or dynamic forces
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02TCLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Complex Calculations (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明还提供一种弹性波数值仿真分析方法与系统,系统包括边界定义、材料定义和绑定、网格生成和检查、求解器前处理、求解器求解、求解模型打开和保存、交错网格显示、结果显示和脚本处理模块。边界定义、材料定义和绑定、网格生成和检查、求解器前处理、求解器求解、求解模型打开和保存、交错网格显示、结果显示以上每一步操作均生成脚本并保存。本发明的一种弹性波数值仿真分析方法与系统,建立了任意形状组合结构、任意多个吸收边界条件、任意区域不均匀材料参数的通用求解器,并采用了有限差分法,能扩展弹性波在复杂结构健康检测监测的应用。

Description

一种弹性波数值仿真分析方法与系统
技术领域
本发明属于计算机仿真技术领域,涉及一种弹性波数值仿真分析方法与系统。
背景技术
复杂结构中的弹性波传播分析是解决结构健康检测监测的重要技术手段,对于包含吸收边界条件的任意形状的不均匀材料参数的结构中的弹性波传播问题,传统有限元法在建模速度和计算速度方面都面临挑战。
交错网格有限差分法作为研究弹性波传播的有效方法,由于不用构建大型刚度矩阵,在建模速度、并行求解效率方面都有较大优势。为了扩展弹性波在结构中的应用,亟待设计一种弹性波数值仿真分析方法与系统,以建立解决包括任意结构形状、任意多个吸收边界条件、任意区域不均匀材料参数的通用求解器,能恰好地任意形状组合结构中复杂的弹性波传播分析问题,比如土木工程中的方桩、管桩、桥梁桩-平台系统损伤识别评估分析,以及航空钢结构弹性波传播的损伤识别问题。
发明内容
本发明针对上述现有技术的不足,提供一种弹性波数值仿真分析方法与系统,以有效地进行建模、求解和结果显示,本发明的方法适应任意结构形状、快速建模、快速求解和结果显示。
本发明的技术解决方案是:一种弹性波数值仿真分析方法,所述仿真分析方法包括以下步骤:
材料参数定义:材料参数包括:材料名称、拉梅常数(λ、μ、κ)、弹性模量E、泊松比υ、密度ρ、横波波速Vs、纵波波速Vp等参数。材料参数实现根据弹性模量E、泊松比υ、密度ρ计算拉梅常数(λ、μ、κ)、横波波速Vs、纵波波速Vp的功能。
分析边界定义:通过边界对象实现不同类型的边界,边界通过实现判断点是否在边界内和是否在边界上的功能,边界类型包括长方体、圆柱体、自定义边界、组边界、排斥边界等类型;边界组合是多个边界对象的集合;排斥边界包含2个边界对象,实现在一个边界内和不在另外一个边界内的功能。通过边界组对象实现不均匀材料的组合计算区域定义。
计算域边界组创建:创建多个边界对象并对每个边界对象进行材料参数绑定,创建能够适应任意形状计算区域的边界组。作为网格生成的参数输入。
交错网格定义:每个网格点采用6个整数索引标记:[i,j,k],[m,n,l],其中[i,j,k]用来标记计算区域内部网格点,[m,n,l]用来标记边界网格点;对于均匀网格的网格点(x,y,z):x=x0+iΔx/2,y=y0+jΔy/2,z=z0+kΔz/2,(x0,y0,z0)起始坐标,Δx,Δy,Δz分别为网格步长;对于不均匀网格点(xi,yj,zk),分别为网格点[i,j,k]的坐标,xi,yj,zk通过预先定义输入。每个网格点[i,j,k]存储6个方向(X+:[i+1,j,k],X-:[i-1,j,k],Y+:[i,j+1,k],Y-:[i,j-1,k],Z+:[i,j,k+1],Z-:[i,j,k-1])的附近网格点坐标。
交错网格变量定义和存储位置:网格采用交错网格,交错网格变量包括9个:u,v,w,σxx,σyy,σzz,τxy,τyz,τxz;3个速度分量u,v,w分别为X,Y,Z方向的速度分量;σxx,σyy,σzz分别为X,Y,Z方向的正应力分量;τxy,τyz,τxz为剪应力分量。交错网格存储:σxx,σyy,σzz存储在同一网格点;u,v,w分别存储在σxx,σyy,σzz网格点之间;τxy存储在u,v中间;τyz存储在v,w网格点中间;τxz存储在u,w网格点中间。
吸收区域定义:吸收边界通过边界对象定义,并根据网格点坐标计算吸收边界的吸收系数并存储在网格点对象;
自由表面定义:对于任意形状的计算区域,自由表面定义共7种:X(垂直于X轴的平面),Y(垂直于Y轴的平面),Z(垂直于Z轴的平面),XY(垂直于X轴和Y轴的线),YZ(垂直于Y轴和Z轴的线),XZ(垂直于X轴和Z轴的线),XYZ(角点,即三条边界线的交点)。
自由表面外镜像网格点定义:共12种类型:边界X网格点、边界Y网格点、边界Z网格点;边界XY的靠近X轴侧网格点、边界XY的靠近Y轴侧网格点、边界YZ的靠近Y轴侧网格点、边界YZ的靠近Z轴侧网格点、边界XZ的靠近X轴侧网格点、边界XZ的靠近Z轴侧网格点;靠近边界XYZ的靠近X轴侧网格点,靠近边界XYZ的靠近Y轴侧网格点,靠近边界XYZ的靠近Z轴侧网格点。镜像网格点分为三层,分别为自由表面外边界层1,边界层2和边界层3。每个网格点的[m,n,l]存储6个方向(X+:[m+1,n,l],X-:[m-1,n,l],Y+:[m,n+1,l],Y-:[m,n-1,l],Z+:[m,n,l+1],Z-:[m,n,l-1])的附近网格点坐标,网格点的[i,j,k]和自由表面的网格点一致。
网格生成:网格生成采用均匀步长网格生成算法和不均匀步长网格算法两种。均匀网格步长生成算法从初始点(x0,y0,z0)开始,对第一个网格点赋予任意一个变量,对计算区域内的每一个[i,j,k]进行判断是否在边界组内,如果在,创建该网格点对象,并检查周围网格点对象是否已经创建,如果已经创建则建立指针和确定当前网格点存储的变量,并根据边界材料参数赋予网格点材料属性;不均匀步长网格算法根据预定义的X,Y,Z方向的种子,对任意一个网格点(xi,yj,zk),对计算区域内的每一个[i,j,k]进行判断是否在边界组内,如果在,创建该网格点对象,并检查周围网格点对象是否已经创建,如果已经创建则建立指针和确定当前网格点存储的变量,并根据边界材料参数赋予网格点材料属性。网格生成过程还需判断网格点是否在吸收区域内,如果是则根据网格点坐标计算吸收边界的吸收系数并存储在网格点对象。
自由表面外侧镜像网格点生成:镜像网格点包括边界层1,边界层2和边界层3,镜像网格点坐标通过以自由表面位置镜像中心,通过内部镜像的1、2、3层的网格点的坐标镜像确定,并复制材料参数。均匀网格和不均匀网格均采取该算法确定镜像网格点。
初始条件定义:通过边界对象和激发时程数据或激发时程函数定义初始条件,求解器求解过程中根据激发时程数据或激发时程函数的值,对边界对象内的网格点进行幅值。
求解器全局参数定义:求解时间步数m、时间步长dt、并行计算线程数量n、均匀网格步长(dx,dy,dz)和材料参数。
网格点差分格式数据库定义:网格点类型包括:9个分量(包括3个速度分量和6个应力分量),对于计算区域内部和自由表面,有(1+7)×9=72种网格点类型;对于自由表面外镜像网格点,共有6(the type count of boundary layers)×12(the type count ofimage points)×9=648种网格点,一共有72+648=720种网格点类型;另外考虑吸收边界范围内的网格点类型(1+7)×9=72种;一共792种网格点。对每一种网格点类型建立有限差分计算公式,并根据网格点位置(内部网格点、自由表面、镜像网格点第一、二和三层)、网格类型(792种网格类型)对有限差分计算公式进行标记。
网格点差分格式数据库根据计算类型分为等步长均匀材料、等步长不均匀材料、不等步长均匀材料、不等步长不均匀材料4种网格点差分格式数据库。
网格点集合生成:共9种:(1)内部和自由表面应力网格点;(2)边界层1上的应力网格点;(3)边界层2上的应力网格点;(4)边界层3上的应力网格点;(5)初始条件定义的应力网格点;(6)内部和自由表面速度网格点;(7)边界层1上的速度网格点;(8)边界层2上的速度网格点;(9)初始条件定义的速度网格点。
并行求解集合生成:对计算区域内所有网格点,根据并行线程数量n,对网格点集合内的每个网格点,根据网格点差分格式数据库,对每一个网格点存储的速度或应力分量,根据网格类型查找对应的有限差分计算公式,添加到求解集合。对每个网格点集合生成n个求解集合,一共9×n=9n个求解集合;
通用并行求解器计算:创建并行计算线程数量n线程的线程池,对m个时间步的每一个时间步:对9n个并行求解集合,分9次,对n个求解集合生成n个计算线程,提交给线程池进行计算。计算过程中对每一时间步进行结果输出。
结果输出包括计算区域网格点坐标和类型数据输出、任意设定的多个截面速度或应力分量的数据输出、任意设定的网格点的分量数据时程输出;结果输出可选择时间步进行输出。
本发明还提供一种弹性波数值仿真分析系统,包括边界定义、材料定义和绑定、网格生成和检查、求解器前处理、求解器求解、求解模型打开和保存、交错网格显示、结果显示和脚本处理模块。边界定义、材料定义和绑定、网格生成和检查、求解器前处理、求解器求解、求解模型打开和保存、交错网格显示、结果显示以上每一步操作均生成脚本并保存,其中:
边界定义模块输入边界数据,生成长方体、圆柱体、自定义边界、组边界;选取2个边界对象生成排斥边界;选取多个边界生成多个边界对象的集合。实现边界的增加、删除、修改和查询功能。
材料定义和绑定模块,实现边界的三维显示、边界绑定显示、增加、删除、修改和查询功能。
网格生成和检查模块,采用等步长均匀网格生成和不等步长网格生成算法进行网格生成。
求解器前处理模块,包括网格点集合生成和并行求解集合生成功能,生成求解过程需要的并行求解集合数据。
求解器求解模块,实现求解过程状态监控和输出、启动、暂停、重新启动等功能。
求解模型打开和保存模块,实现计算区域的网格点数据、求解过程分量的时程数据的打开保存功能。
交错网格显示模块,对网格点分量、网格连接等进行三维显示,并包括网格选取和显示功能。
结果显示模块包括速度应力分量的任意设定截面结果显示和整个计算区域结果显示。
脚本处理模块包括边界定义、材料定义和绑定、网格生成和检查、求解器前处理、求解器求解、求解模型打开和保存、交错网格显示、结果显示等模块的脚本生成、执行功能;实现脚本的新建、删除、修改保存功能。
本发明的有益效果是:
本发明的一种弹性波数值仿真分析方法与系统,建立了任意形状组合结构、任意多个吸收边界条件、任意区域不均匀材料参数的通用求解器,并采用了有限差分法,能扩展弹性波在复杂结构健康检测监测的应用。
以下将结合附图对本发明的构思、具体结构及产生的技术效果作进一步说明,以充分地了解本发明的目的、特征和效果。
附图说明
图1是本发明方法的流程示意图
具体实现方式
为了使本发明的目的、内容和流程更加清楚,结合附图对本发明具体实施步骤进行详细说明,本发明的具体实施方式不限于此。
一种弹性波数值仿真分析方法,所述仿真分析方法包括以下步骤:
材料参数定义:材料参数包括:材料名称、拉梅常数(λ、μ、κ)、弹性模量E、泊松比υ、密度ρ、横波波速Vs、纵波波速Vp等参数。材料参数实现根据弹性模量E、泊松比υ、密度ρ计算拉梅常数(λ、μ、κ)、横波波速Vs、纵波波速Vp的功能。
分析边界定义:通过边界对象实现不同类型的边界,边界通过实现判断点是否在边界内和是否在边界上的功能,边界类型包括长方体、圆柱体、自定义边界、组边界、排斥边界等类型;边界组合是多个边界对象的集合;排斥边界包含2个边界对象,实现在一个边界内和不在另外一个边界内的功能。通过边界组对象实现不均匀材料的组合计算区域定义。
计算域边界组创建:创建多个边界对象并对每个边界对象进行材料参数绑定,创建能够适应任意形状计算区域的边界组。作为网格生成的参数输入。
交错网格定义:每个网格点采用6个整数索引标记:[i,j,k],[m,n,l],其中[i,j,k]用来标记计算区域内部网格点,[m,n,l]用来标记边界网格点;对于均匀网格的网格点(x,y,z):x=x0+iΔx/2,y=y0+jΔy/2,z=z0+kΔz/2,(x0,y0,z0)起始坐标,Δx,Δy,Δz分别为网格步长;对于不均匀网格点(xi,yj,zk),分别为网格点[i,j,k]的坐标,xi,yj,zk通过预先定义输入。每个网格点[i,j,k]存储6个方向(X+:[i+1,j,k],X-:[i-1,j,k],Y+:[i,j+1,k],Y-:[i,j-1,k],Z+:[i,j,k+1],Z-:[i,j,k-1])的附近网格点坐标。
交错网格变量定义和存储位置:网格采用交错网格,交错网格变量包括9个:u,v,w,σxx,σyy,σzz,τxy,τyz,τxz;3个速度分量u,v,w分别为X,Y,Z方向的速度分量;σxx,σyy,σzz分别为X,Y,Z方向的正应力分量;τxy,τyz,τxz为剪应力分量。交错网格存储:σxx,σyy,σzz存储在同一网格点;u,v,w分别存储在σxx,σyy,σzz网格点之间;τxy存储在u,v中间;τyz存储在v,w网格点中间;τxz存储在u,w网格点中间。
吸收区域定义:吸收边界通过边界对象定义,并根据网格点坐标计算吸收边界的吸收系数并存储在网格点对象;
自由表面定义:对于任意形状的计算区域,自由表面定义共7种:X(垂直于X轴的平面),Y(垂直于Y轴的平面),Z(垂直于Z轴的平面),XY(垂直于X轴和Y轴的线),YZ(垂直于Y轴和Z轴的线),XZ(垂直于X轴和Z轴的线),XYZ(角点,即三条边界线的交点)。
自由表面外镜像网格点定义:共12种类型:边界X网格点、边界Y网格点、边界Z网格点;边界XY的靠近X轴侧网格点、边界XY的靠近Y轴侧网格点、边界YZ的靠近Y轴侧网格点、边界YZ的靠近Z轴侧网格点、边界XZ的靠近X轴侧网格点、边界XZ的靠近Z轴侧网格点;靠近边界XYZ的靠近X轴侧网格点,靠近边界XYZ的靠近Y轴侧网格点,靠近边界XYZ的靠近Z轴侧网格点。镜像网格点分为三层,分别为自由表面外边界层1,边界层2和边界层3。每个网格点的[m,n,l]存储6个方向(X+:[m+1,n,l],X-:[m-1,n,l],Y+:[m,n+1,l],Y-:[m,n-1,l],Z+:[m,n,l+1],Z-:[m,n,l-1])的附近网格点坐标,网格点的[i,j,k]和自由表面的网格点一致。
网格生成:网格生成采用均匀步长网格生成算法和不均匀步长网格算法两种。均匀网格步长生成算法从初始点(x0,y0,z0)开始,对第一个网格点赋予任意一个变量,对计算区域内的每一个[i,j,k]进行判断是否在边界组内,如果在,创建该网格点对象,并检查周围网格点对象是否已经创建,如果已经创建则建立指针和确定当前网格点存储的变量,并根据边界材料参数赋予网格点材料属性;不均匀步长网格算法根据预定义的X,Y,Z方向的种子,对任意一个网格点(xi,yj,zk),对计算区域内的每一个[i,j,k]进行判断是否在边界组内,如果在,创建该网格点对象,并检查周围网格点对象是否已经创建,如果已经创建则建立指针和确定当前网格点存储的变量,并根据边界材料参数赋予网格点材料属性。网格生成过程还需判断网格点是否在吸收区域内,如果是则根据网格点坐标计算吸收边界的吸收系数并存储在网格点对象。
自由表面外侧镜像网格点生成:镜像网格点包括边界层1,边界层2和边界层3,镜像网格点坐标通过以自由表面位置镜像中心,通过内部镜像的1、2、3层的网格点的坐标镜像确定,并复制材料参数。均匀网格和不均匀网格均采取该算法确定镜像网格点。
初始条件定义:通过边界对象和激发时程数据或激发时程函数定义初始条件,求解器求解过程中根据激发时程数据或激发时程函数的值,对边界对象内的网格点进行幅值。
求解器全局参数定义:求解时间步数m、时间步长dt、并行计算线程数量n、均匀网格步长(dx,dy,dz)和材料参数。
网格点差分格式数据库定义:网格点类型包括:9个分量(包括3个速度分量和6个应力分量),对于计算区域内部和自由表面,有(1+7)×9=72种网格点类型;对于自由表面外镜像网格点,共有6(the type count of boundary layers)×12(the type count ofimage points)×9=648种网格点,一共有72+648=720种网格点类型;另外考虑吸收边界范围内的网格点类型(1+7)×9=72种;一共792种网格点。对每一种网格点类型建立有限差分计算公式,并根据网格点位置(内部网格点、自由表面、镜像网格点第一、二和三层)、网格类型(792种网格类型)对有限差分计算公式进行标记。
网格点差分格式数据库根据计算类型分为等步长均匀材料、等步长不均匀材料、不等步长均匀材料、不等步长不均匀材料4种网格点差分格式数据库。
网格点集合生成:共9种:(1)内部和自由表面应力网格点;(2)边界层1上的应力网格点;(3)边界层2上的应力网格点;(4)边界层3上的应力网格点;(5)初始条件定义的应力网格点;(6)内部和自由表面速度网格点;(7)边界层1上的速度网格点;(8)边界层2上的速度网格点;(9)初始条件定义的速度网格点。
并行求解集合生成:对计算区域内所有网格点,根据并行线程数量n,对网格点集合内的每个网格点,根据网格点差分格式数据库,对每一个网格点存储的速度或应力分量,根据网格类型查找对应的有限差分计算公式,添加到求解集合。对每个网格点集合生成n个求解集合,一共9×n=9n个求解集合;
通用并行求解器计算:创建并行计算线程数量n线程的线程池,对m个时间步的每一个时间步:对9n个并行求解集合,分9次,对n个求解集合生成n个计算线程,提交给线程池进行计算。计算过程中对每一时间步进行结果输出。
结果输出包括计算区域网格点坐标和类型数据输出、任意设定的多个截面速度或应力分量的数据输出、任意设定的网格点的分量数据时程输出;结果输出可选择时间步进行输出。
本发明还提供一种弹性波数值仿真分析系统,包括边界定义、材料定义和绑定、网格生成和检查、求解器前处理、求解器求解、求解模型打开和保存、交错网格显示、结果显示和脚本处理模块。边界定义、材料定义和绑定、网格生成和检查、求解器前处理、求解器求解、求解模型打开和保存、交错网格显示、结果显示以上每一步操作均生成脚本并保存,其中:
边界定义模块输入边界数据,生成长方体、圆柱体、自定义边界、组边界;选取2个边界对象生成排斥边界;选取多个边界生成多个边界对象的集合。实现边界的增加、删除、修改和查询功能。
材料定义和绑定模块,实现边界的三维显示、边界绑定显示、增加、删除、修改和查询功能。
网格生成和检查模块,采用等步长均匀网格生成和不等步长网格生成算法进行网格生成。
求解器前处理模块,包括网格点集合生成和并行求解集合生成功能,生成求解过程需要的并行求解集合数据。
求解器求解模块,实现求解过程状态监控和输出、启动、暂停、重新启动等功能。
求解模型打开和保存模块,实现计算区域的网格点数据、求解过程分量的时程数据的打开保存功能。
交错网格显示模块,对网格点分量、网格连接等进行三维显示,并包括网格选取和显示功能。
结果显示模块包括速度应力分量的任意设定截面结果显示和整个计算区域结果显示。
脚本处理模块包括边界定义、材料定义和绑定、网格生成和检查、求解器前处理、求解器求解、求解模型打开和保存、交错网格显示、结果显示等模块的脚本生成、执行功能;实现脚本的新建、删除、修改保存功能。
以上详细描述了本发明的较佳具体实施例。应当理解,本领域的普通技术人员无需创造性劳动就可以根据本发明的构思做出诸多修改和变化。因此,凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案,皆应在由权利要求书所确定的保护范围内。

Claims (6)

1.一种弹性波数值仿真分析方法,其特征在于,所述仿真分析方法包括以下步骤:
S1、材料参数定义及分析边界定义,其中分析边界定义:通过边界对象实现不同类型的边界,边界通过实现判断点是否在边界内和是否在边界上的功能,边界类型包括长方体、圆柱体、自定义边界、组边界、排斥边界的类型;边界组合是多个边界对象的集合;排斥边界包含2个边界对象,实现在一个边界内和不在另外一个边界内的功能,通过边界组对象实现不均匀材料的组合计算区域定义;
S2、计算域边界组创建:创建多个边界对象并对每个边界对象进行材料参数绑定,创建能够适应任意形状计算区域的边界组,作为网格生成的参数输入;
S3、交错网格定义:每个网格点采用6个整数索引标记:[i,j,k],[m,n,l],其中[i,j,k]用来标记计算区域内部网格点,[m,n,l]用来标记边界网格点;对于均匀网格的网格点(x,y,z):x=x0+iΔx/2,y=y0+jΔy/2,z=z0+kΔz/2,(x0,y0,z0)起始坐标,Δx,Δy,Δz分别为网格步长;对于不均匀网格点(xi,yj,zk),分别为网格点[i,j,k]的坐标,xi,yj,zk通过预先定义输入,每个网格点[i,j,k]存储6个方向(X+:[i+1,j,k],X-:[i-1,j,k],Y+:[i,j+1,k],Y-:[i,j-1,k],Z+:[i,j,k+1],Z-:[i,j,k-1])的附近网格点坐标;
S4、交错网格变量定义和存储位置:网格采用交错网格,交错网格变量包括9个:u,v,w,σxx,σyy,σzz,τxy,τyz,τxz;3个速度分量u,v,w分别为X,Y,Z方向的速度分量;σxx,σyy,σzz分别为X,Y,Z方向的正应力分量;τxy,τyz,τxz为剪应力分量,交错网格存储:σxx,σyy,σzz存储在同一网格点;u,v,w分别存储在σxx,σyy,σzz网格点之间;τxy存储在u,v中间;τyz存储在v,w网格点中间;τxz存储在u,w网格点中间;
S5、吸收区域定义:吸收边界通过边界对象定义,并根据网格点坐标计算吸收边界的吸收系数并存储在网格点对象;
S6、自由表面定义;
S7、自由表面外镜像网格点定义:共12种类型:边界X网格点、边界Y网格点、边界Z网格点;边界XY的靠近X轴侧网格点、边界XY的靠近Y轴侧网格点、边界YZ的靠近Y轴侧网格点、边界YZ的靠近Z轴侧网格点、边界XZ的靠近X轴侧网格点、边界XZ的靠近Z轴侧网格点;靠近边界XYZ的靠近X轴侧网格点,靠近边界XYZ的靠近Y轴侧网格点,靠近边界XYZ的靠近Z轴侧网格点,镜像网格点分为三层,分别为自由表面外边界层1,边界层2和边界层3,每个网格点的[m,n,l]存储6个方向(X+:[m+1,n,l],X-:[m-1,n,l],Y+:[m,n+1,l],Y-:[m,n-1,l],Z+:[m,n,l+1],Z-:[m,n,l-1])的附近网格点坐标,网格点的[i,j,k]和自由表面的网格点一致;
S8、网格生成:网格生成采用均匀步长网格生成算法和不均匀步长网格算法两种,均匀网格步长生成算法从初始点(x0,y0,z0)开始,对第一个网格点赋予任意一个变量,对计算区域内的每一个[i,j,k]进行判断是否在边界组内,如果在,创建该网格点对象,并检查周围网格点对象是否已经创建,如果已经创建则建立指针和确定当前网格点存储的变量,并根据边界材料参数赋予网格点材料属性;不均匀步长网格算法根据预定义的X,Y,Z方向的种子,对任意一个网格点(xi,yj,zk),对计算区域内的每一个[i,j,k]进行判断是否在边界组内,如果在,创建该网格点对象,并检查周围网格点对象是否已经创建,如果已经创建则建立指针和确定当前网格点存储的变量,并根据边界材料参数赋予网格点材料属性,网格生成过程还需判断网格点是否在吸收区域内,如果是则根据网格点坐标计算吸收边界的吸收系数并存储在网格点对象;
S9、自由表面外侧镜像网格点生成:镜像网格点包括边界层1,边界层2和边界层3,镜像网格点坐标通过以自由表面位置镜像中心,通过内部镜像的1、2、3层的网格点的坐标镜像确定,并复制材料参数,均匀网格和不均匀网格均采取该算法确定镜像网格点;
S10、初始条件定义:通过边界对象和激发时程数据或激发时程函数定义初始条件,求解器求解过程中根据激发时程数据或激发时程函数的值,对边界对象内的网格点进行幅值;
S11、求解器全局参数定义:求解时间步数m、时间步长dt、并行计算线程数量n、均匀网格步长(dx,dy,dz)和材料参数;
S12、网格点差分格式数据库定义:网格点类型包括:9个分量,包括3个速度分量和6个应力分量,对于计算区域内部和自由表面,有(1+7)×9=72种网格点类型;对于自由表面外镜像网格点,共有6×12×9=648种网格点,一共有72+648=720种网格点类型;另外考虑吸收边界范围内的网格点类型(1+7)×9=72种;一共792种网格点,对每一种网格点类型建立有限差分计算公式,并根据网格点位置、网格类型对有限差分计算公式进行标记,其中并根据网格点位置包括内部网格点、自由表面、镜像网格点第一、二和三层,网格类型包括792种网格类型;
S13、网格点差分格式数据库根据计算类型分为等步长均匀材料、等步长不均匀材料、不等步长均匀材料、不等步长不均匀材料4种网格点差分格式数据库;
S14、网格点集合生成;
S15、并行求解集合生成:对计算区域内所有网格点,根据并行线程数量n,对网格点集合内的每个网格点,根据网格点差分格式数据库,对每一个网格点存储的速度或应力分量,根据网格类型查找对应的有限差分计算公式,添加到求解集合,对每个网格点集合生成n个求解集合,一共9×n=9n个求解集合;
S16、通用并行求解器计算:创建并行计算线程数量n线程的线程池,对m个时间步的每一个时间步:对9n个并行求解集合,分9次,对n个求解集合生成n个计算线程,提交给线程池进行计算,计算过程中对每一时间步进行结果输出;结果输出包括计算区域网格点坐标和类型数据输出、任意设定的多个截面速度或应力分量的数据输出、任意设定的网格点的分量数据时程输出;结果输出可选择时间步进行输出。
2.根据权利要求1所述的一种弹性波数值仿真分析方法,其特征在于,所述材料参数定义中,材料参数包括:材料名称、拉梅常数(λ、μ、κ)、弹性模量E、泊松比υ、密度ρ、横波波速Vs、纵波波速Vp;材料参数实现根据弹性模量E、泊松比υ、密度ρ计算拉梅常数(λ、μ、κ)、横波波速Vs、纵波波速Vp的功能。
3.根据权利要求1所述的一种弹性波数值仿真分析方法,其特征在于,所述网格点集合生成共9种,具体为:(1)内部和自由表面应力网格点;(2)边界层1上的应力网格点;(3)边界层2上的应力网格点;(4)边界层3上的应力网格点;(5)初始条件定义的应力网格点;(6)内部和自由表面速度网格点;(7)边界层1上的速度网格点;(8)边界层2上的速度网格点;(9)初始条件定义的速度网格点。
4.根据权利要求1所述的一种弹性波数值仿真分析方法,其特征在于,所述自由表面定义,对于任意形状的计算区域,自由表面定义共7种:垂直于X轴的平面,垂直于Y轴的平面,垂直于Z轴的平面,垂直于X轴和Y轴的线,垂直于Y轴和Z轴的线,垂直于X轴和Z轴的线,角点,即三条边界线的交点。
5.一种弹性波数值仿真分析系统,其特征在于,包括边界定义模块、材料定义和绑定模块、网格生成和检查模块、求解器前处理模块、求解器求解模块、求解模型打开和保存模块、交错网格显示模块、结果显示模块和脚本处理模块,其中:
所述边界定义模块输入边界数据,生成长方体、圆柱体、自定义边界、组边界;选取2个边界对象生成排斥边界;选取多个边界生成多个边界对象的集合,实现边界的增加、删除、修改和查询功能;
所述材料定义和绑定模块,实现边界的三维显示、边界绑定显示、增加、删除、修改和查询功能;
所述网格生成和检查模块,采用等步长均匀网格生成和不等步长网格生成算法进行网格生成;
所述求解器前处理模块,包括网格点集合生成和并行求解集合生成功能,生成求解过程需要的并行求解集合数据;
所述求解器求解模块,实现求解过程状态监控和输出、启动、暂停、重新启动的功能;
所述求解模型打开和保存模块,实现计算区域的网格点数据、求解过程分量的时程数据的打开保存功能;
所述交错网格显示模块,对网格点分量、网格连接等进行三维显示,并包括网格选取和显示功能;
所述结果显示模块包括速度应力分量的任意设定截面结果显示和整个计算区域结果显示;
所述脚本处理模块包括边界定义、材料定义和绑定、网格生成和检查、求解器前处理、求解器求解、求解模型打开和保存、交错网格显示、结果显示等模块的脚本生成、执行功能;实现脚本的新建、删除、修改保存功能。
6.根据权利要求5所述的一种弹性波数值仿真分析系统,其特征在于,所述边界定义模块、材料定义和绑定模块、网格生成和检查模块、求解器前处理模块、求解器求解模块、求解模型打开和保存模块、交错网格显示模块、结果显示模块和脚本处理模块每一步操作均生成脚本并保存。
CN202210412326.XA 2022-04-19 2022-04-19 一种弹性波数值仿真分析方法与系统 Active CN114818422B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210412326.XA CN114818422B (zh) 2022-04-19 2022-04-19 一种弹性波数值仿真分析方法与系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210412326.XA CN114818422B (zh) 2022-04-19 2022-04-19 一种弹性波数值仿真分析方法与系统

Publications (2)

Publication Number Publication Date
CN114818422A true CN114818422A (zh) 2022-07-29
CN114818422B CN114818422B (zh) 2024-03-22

Family

ID=82505101

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210412326.XA Active CN114818422B (zh) 2022-04-19 2022-04-19 一种弹性波数值仿真分析方法与系统

Country Status (1)

Country Link
CN (1) CN114818422B (zh)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO1999056153A1 (en) * 1998-04-27 1999-11-04 Phillips Petroleum Company Method and apparatus for cancelling reflections in wave propagation models
FR3031210A1 (fr) * 2014-12-31 2016-07-01 Landmark Graphics Corp Simulation d'onde elastique sismique pour milieu isotrope transversalement incline utilisant une grille decalee de lebedev
CN106250102A (zh) * 2015-06-12 2016-12-21 中国石油化工股份有限公司 交错网格有限差分正演模拟优化的方法
CN109725351A (zh) * 2018-12-18 2019-05-07 中国石油天然气集团有限公司 一种3d弹性波混合吸收边界条件的确定方法、装置及系统
CN110032756A (zh) * 2019-02-27 2019-07-19 上海交通大学 基于流函数分数坐标系变换的流动边界层数值分析方法
CN111859766A (zh) * 2020-07-28 2020-10-30 深圳拳石科技发展有限公司 可变计算域的拉格朗日积分点有限元数值仿真系统及方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO1999056153A1 (en) * 1998-04-27 1999-11-04 Phillips Petroleum Company Method and apparatus for cancelling reflections in wave propagation models
FR3031210A1 (fr) * 2014-12-31 2016-07-01 Landmark Graphics Corp Simulation d'onde elastique sismique pour milieu isotrope transversalement incline utilisant une grille decalee de lebedev
CN106250102A (zh) * 2015-06-12 2016-12-21 中国石油化工股份有限公司 交错网格有限差分正演模拟优化的方法
CN109725351A (zh) * 2018-12-18 2019-05-07 中国石油天然气集团有限公司 一种3d弹性波混合吸收边界条件的确定方法、装置及系统
CN110032756A (zh) * 2019-02-27 2019-07-19 上海交通大学 基于流函数分数坐标系变换的流动边界层数值分析方法
CN111859766A (zh) * 2020-07-28 2020-10-30 深圳拳石科技发展有限公司 可变计算域的拉格朗日积分点有限元数值仿真系统及方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
张明财;熊章强;张大洲;: "基于MPI的三维瑞雷面波有限差分并行模拟", 石油物探, no. 04, 25 July 2013 (2013-07-25) *
裴正林: "任意起伏地表弹性波方程交错网格高阶有限差分法数值模拟", 石油地球物理勘探, no. 06, 15 December 2004 (2004-12-15) *
陈可洋;: "完全匹配层吸收边界条件在弹性波波场分离数值模拟中的应用", 石油工业计算机应用, no. 01, 20 March 2010 (2010-03-20) *

Also Published As

Publication number Publication date
CN114818422B (zh) 2024-03-22

Similar Documents

Publication Publication Date Title
Abedi et al. An h-adaptive spacetime-discontinuous Galerkin method for linear elastodynamics
US20130151551A1 (en) Computer-implemented method of geometric feature detection
CN110765679B (zh) 一种基于有限元模型与SVM回归算法的大坝监测web展示方法
Khokhlov et al. On the class of compact grid-characteristic schemes
CN114722745B (zh) 湍流壁面距离计算方法、装置、计算机设备和存储介质
US20170017739A1 (en) Post-processing system for finite element analysis
CN102495926B (zh) 三维原始模型的检验方法及装置
CN114861500B (zh) 基于三维点云自动生成隧道结构有限元模型的方法及系统
CN105608283B (zh) 一种飞机结构强度单元合并方法与装置
CN114818422A (zh) 一种弹性波数值仿真分析方法与系统
CN112433003A (zh) 一种用于t形结构件超声波检测的三维仿真方法
Krishnan et al. BOOLE: A boundary evaluation system for boolean combinations of sculptured solids
Liu et al. Fast solving the Cauchy problems of Poisson equation in an arbitrary three-dimensional domain
Nie et al. Development of an object-oriented finite element program with adaptive mesh refinement for multi-physics applications
CN112560386B (zh) 一种大规模复杂版图电阻提取加速方法
US11113436B1 (en) Method for simulation of flow in fluid flow network having one-dimensional and multi-dimensional flow components
Graysmith et al. Automated procedures for Boolean operations on finite element meshes
EP3451191A1 (en) Computer implemented method for manipulating a numerical model of a 3d domain
CN117494539B (zh) 物面流体仿真中粒子的近邻搜索方法、装置、电子设备
Kercher Generic implementation of CAD models for nuclear simulation
US9760662B1 (en) Simulating performance, reliability, weight, and cost using a physical model
Sabino Discrete oriented polytopes with orthogonal bases for the construction of tighter Bounding Volume Hierarchies
Oanta Integration of the original software applications for mechanical engineering
CN115587509A (zh) 一种任意缺陷漏磁场的快速计算方法、设备及存储设备
Franklin et al. Computing intersection areas of overlaid 2D meshes

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant