CN114357907B - 一种适用于拉格朗日型粒子类数值模拟的并行方法 - Google Patents

一种适用于拉格朗日型粒子类数值模拟的并行方法 Download PDF

Info

Publication number
CN114357907B
CN114357907B CN202210014743.9A CN202210014743A CN114357907B CN 114357907 B CN114357907 B CN 114357907B CN 202210014743 A CN202210014743 A CN 202210014743A CN 114357907 B CN114357907 B CN 114357907B
Authority
CN
China
Prior art keywords
sub
calculation
particle
region
calculation region
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.)
Active
Application number
CN202210014743.9A
Other languages
English (en)
Other versions
CN114357907A (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.)
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Original Assignee
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
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 Computational Aerodynamics Institute of China Aerodynamics Research and Development Center filed Critical Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority to CN202210014743.9A priority Critical patent/CN114357907B/zh
Publication of CN114357907A publication Critical patent/CN114357907A/zh
Application granted granted Critical
Publication of CN114357907B publication Critical patent/CN114357907B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • Y—GENERAL 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
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明适用于流体力学技术领域，提供了一种适用于拉格朗日型粒子类数值模拟的并行方法，包括获取计算区域内的粒子，将计算区域划分为多个子计算区域，每个子计算区域内的粒子数量相同；每个子计算区域对应的线程获取当前子计算区域中的粒子信息、计算参数信息，以及线程属性信息；设置每个子计算区域的联结区；获取当前时刻与各子计算区域相邻子计算区域的联结区内粒子信息，并对下一时刻各子计算区域中的粒子信息进行更新；当下一时刻的各子计算区域粒子越过当前各子计算区域的边界时，将越界的粒子发送至对应的子计算区域，该方法提高了粒子类数值模拟算法的计算效率和内存供给量。

Description

一种适用于拉格朗日型粒子类数值模拟的并行方法
技术领域
本发明涉及流体力学技术领域，尤其是涉及一种适用于拉格朗日型粒子类数值模拟的并行方法。
背景技术
对于理工类的科学研究，主要有理论、实验和数值模拟这几种科研方法。理论科研方法遇到的难点主要包括理论模型建立的困难和理论验证的困难，实验科研方法遇到的难点主要包括技术实现难度大、耗费成本高。数值模拟通过在计算机上执行针对特定问题的数学方程和公式，求得计算结果并进行相应分析，具有高效低成本、获取数据全面的优点。目前主要的数值模拟方法可分为基于网格的数值模拟方法和基于粒子的数值模拟方法两类。有限差分方法、有限体积方法、有限元方法等基于网格的数值模拟方法发展时间更早，发展得更为成熟，已经被大量运用于解决结构力学、流体力学、爆炸力学等领域。
光滑粒子流体动力学方法、再生核质点方法等拉格朗日型粒子类方法由于无网格限制，无需因为物质运动而对网格进行复杂操作防止网格破碎，且能自然地捕捉界面，所以在模拟碰撞冲击、造浪等界面变形问题时有较大优势。但由于粒子类方法发展时间尚短，且目前更多聚焦的是方法本身精度、稳定性、边界条件等方面的发展，在大规模工程计算方面还较弱。以光滑粒子流体动力学方法方法为例，其在每个计算步内都需要找出每个粒子周边的相邻粒子，从而实现粒子配对计算，这导致其计算效率较低，目前采用的串行计算更是只能计算二维问题或是简单的三维问题，对于上百万的粒子数，由于单线程内存限制甚至可能导致没法计算或计算时间过长无法接受。近期用的一种较多的基于粒子分解方式的并行算法，其是一套支持跨平台共享内存方式的多线程并发的编程框架，优点是对源代码改动小，实施方便，但其缺点是所有线程共享内存空间，受硬件制约较大，且主要针对的是程序循环代码处，故也无法处理过大规模的计算问题，一般来说千万量级粒子问题已经很难处理。
发明内容
本发明的目的是提供一种适用于拉格朗日型粒子类数值模拟的并行方法，来解决现有技术中存在的上述技术问题，包括如下步骤：
步骤100：获取计算区域内的粒子，将计算区域划分为多个子计算区域，每个子计算区域内的粒子数量相同；每个子计算区域对应的线程获取当前子计算区域中的粒子信息、计算参数信息，以及线程属性信息；
步骤200：设置每个子计算区域的联结区；
步骤300：获取当前时刻与各子计算区域相邻子计算区域的联结区内粒子信息，并对下一时刻各子计算区域中的粒子信息进行更新；
步骤400：当下一时刻的各子计算区域粒子越过当前各子计算区域的边界时，将更新后下一时刻各子计算区域中越界的粒子发送至对应的子计算区域中。
进一步地，步骤200中：每个子计算区域的联结区的宽度B不小于阈值宽度R。
进一步地，宽度阈值R的计算方法如下：
Figure 955868DEST_PATH_IMAGE001
其中，k为预设常数，
Figure 545112DEST_PATH_IMAGE002
表示两个相邻子计算区域的联结区的宽度值，n表示第n个子计算区域，m表示第m个子计算区域，第n个子计算区域与第m个子计算区域相邻，
Figure 59270DEST_PATH_IMAGE003
，n=0~N-1，m=0~N-1，N表示子计算区域的总数。
进一步地，
Figure 126583DEST_PATH_IMAGE002
的计算方法如下：
Figure 675376DEST_PATH_IMAGE004
其中，
Figure 169943DEST_PATH_IMAGE005
表示第n个子计算区域中所有粒子中的最大光滑长度，
Figure 171397DEST_PATH_IMAGE006
表示第n个子计算区域中的第i个粒子的光滑长度，i表示第n个子计算区域中的第i个粒子；
Figure 42401DEST_PATH_IMAGE007
表示第个子计算区域中所有粒子中的最大光滑长度，
Figure 914542DEST_PATH_IMAGE008
表示第m个子计算区域中的第j个粒子的光滑长度，j表示第个子计算区域中的第j个粒子。
进一步地，当计算区域中的粒子直径只随空间改变，不随时间改变时，
Figure 580010DEST_PATH_IMAGE006
和
Figure 68760DEST_PATH_IMAGE008
的计算方法如下：
Figure 477875DEST_PATH_IMAGE009
Figure 938944DEST_PATH_IMAGE010
其中，
Figure 837630DEST_PATH_IMAGE011
表示特定系数，d i 表示第i个粒子的直径，d j 表示第j个粒子的直径。
进一步地，
Figure 16938DEST_PATH_IMAGE012
和
Figure 229745DEST_PATH_IMAGE013
的计算方法如下：
Figure 607637DEST_PATH_IMAGE014
或
Figure 614907DEST_PATH_IMAGE015
其中，
Figure 281512DEST_PATH_IMAGE016
表示第n个子计算区域中第i个粒子的面积，
Figure 94747DEST_PATH_IMAGE017
表示第n个子计算区域中第i个粒子的体积；
Figure 264828DEST_PATH_IMAGE018
或
Figure 770896DEST_PATH_IMAGE019
其中，
Figure 393638DEST_PATH_IMAGE020
表示第m个子计算区域中第j个粒子的面积，
Figure 744985DEST_PATH_IMAGE021
表示第m个子计算区域中第j个粒子的体积。
进一步地，当计算区域中的粒子直径随时间变化时，
Figure 831890DEST_PATH_IMAGE006
和
Figure 446542DEST_PATH_IMAGE008
随时间的更新方法如下：
Figure 353318DEST_PATH_IMAGE022
Figure 508356DEST_PATH_IMAGE023
其中，t为计算时长，D为整个模拟对象的维度，D=1，2，3，
Figure 449767DEST_PATH_IMAGE024
表示第n个子计算区域中的第i个粒子的密度，
Figure 235320DEST_PATH_IMAGE025
表示第m个子计算区域中的第j个粒子的密度。
进一步地，步骤S100中计算区域划分为多个子计算区域的步骤如下：
获取计算区域所在的坐标系x-y-z；
分别遍历坐标系中x、y、z方向上粒子，在每个1/X、1/Y、1/Z粒子数所对应的x轴位置进行切割，其中X为x方向上预设的计算区域切割数量，Y为y方向上预设的计算区域切割数量，Z为z方向上预设的计算区域切割数量。
进一步地，步骤S300中对当前线程中的粒子信息进行更新的方法如下：
基于粒子类数值算法计算当前时刻下各子计算区域粒子信息的变化率；
根据当前时刻下各子计算区域粒子信息的变化率，更新获得下一时刻的各子计算区域粒子信息。
进一步地，步骤400中，下一时刻的各子计算区域粒子越过当前各子计算区域联结区的边界的判断方法如下：
获取下一时刻各计算区域更新后粒子的坐标信息；
获取当前时刻各子计算区域联结区的边界的坐标信息；
当粒子的坐标信息值大于边界的坐标信息时，则该粒子越过当前联结区的边界。
本发明相对于现有技术至少具有如下技术效果：
（1）本发明提供的一种适用于拉格朗日型粒子类数值模拟的并行方法，通过将计算区域划分为多个子计算区域，每个子计算区域配置有对应的线程，各线程同时运行进行粒子的数值计算，并且在每个子计算区域中设置了联结区，该方法可以解决大规模，如亿量级粒子规模问题计算，提高了粒子类数值模拟算法的计算效率和内存供给量，且根据粒子类数值模拟方法的特殊性，通过设置的联结区提高了粒子信息在计算过程中传递的完整度，很大程度上提高了粒子类数值方法的工程应用能力，如提升了粒子类算法在处理航空、航天、航海、预防自然灾害等领域的工程应用能力。
（2）本发明提供的一种适用于拉格朗日型粒子类数值模拟的并行方法中，在每个子计算区域中设置了联结区，并且设置了每个联结区的宽度，该宽度可以根据粒子的光滑长度进行动态调整，且每个联结区的宽度也是结合核函数计算性质及各子计算区域包含粒子信息、相邻子计算区域信息进行设置，保证了各粒子插值计算所需周边粒子信息不被子区域边界截断的同时尽量减少了信息交互消耗，提高了数值模拟的准确性。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对本发明实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面所描述的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明中适用于拉格朗日型粒子类数值模拟的并行方法的流程图；
图2是本发明中的将计算区域划分为多个子计算区域的示意图；
图3是本发明中的设置每个子计算区域联结区的示意图；
图4是本发明中的相邻子计算区域联结区粒子信息获取的示意图。
具体实施方式
以下的说明提供了许多不同的实施例、或是例子，用来实施本发明的不同特征。以下特定例子所描述的元件和排列方式，仅用来精简的表达本发明，其仅作为例子，而并非用以限制本发明。
在下文中将参考附图对本发明的各方面进行更充分的描述。然而，本发明可以具体化成许多不同形式且不应解释为局限于贯穿本发明所呈现的任何特定结构或功能。相反地，提供这些方面将使得本发明周全且完整，并且本发明将给本领域技术人员充分地传达本发明的范围。基于本文所教导的内容，本领域的技术人员应意识到，无论是单独还是结合本发明的任何其它方面实现本文所公开的任何方面，本发明的范围旨在涵盖本文中所公开的任何方面。例如，可以使用本文所提出任意数量的装置或者执行方法来实现。另外，除了本文所提出本发明的多个方面之外，本发明的范围更旨在涵盖使用其它结构、功能或结构和功能来实现的装置或方法。应可理解，其可通过权利要求的一或多个元件具体化本文所公开的任何方面。
在此使用的术语仅仅是为了描述具体实施例，而并非意在限制本公开。在此使用的术语“包括”、“包含”等表明了所述特征、步骤、操作和/或部件的存在，但是并不排除存在或添加一个或多个其他特征、步骤、操作或部件。
在此使用的所有术语(包括技术和科学术语)具有本领域技术人员通常所理解的含义，除非另外定义。应注意，这里使用的术语应解释为具有与本说明书的上下文相一致的含义，而不应以理想化或过于刻板的方式来解释。
如图1所示，本发明的目的是提供一种适用于拉格朗日型粒子类数值模拟的并行方法，来解决现有技术中存在的上述技术问题，包括如下步骤：
步骤S100：获取计算区域内的粒子，将计算区域划分为多个子计算区域，每个子计算区域内的粒子数量相同；每个子计算区域对应的线程获取当前子计算区域中的粒子信息、计算参数信息，以及线程属性信息；
步骤S200：设置每个子计算区域的联结区；
步骤S300：获取当前时刻与各子计算区域相邻子计算区域的联结区内粒子信息，并对下一时刻各子计算区域中的粒子信息进行更新；
步骤S400：当下一时刻的各子计算区域粒子越过当前各子计算区域的边界时，将更新后下一时刻各子计算区域中越界的粒子发送至对应的子计算区域中。
上述方案中，在步骤S100前，预先在需要计算的区域模拟生成初始粒子，用于后续的数值模拟计算。计算机获取需要计算区域内的粒子，并将需要计算区域划分为多个子计算区域，具体需要将计算区域划分为多少个子计算区域，可以结合粒子的总数量、计算机中线程的数量、计算的时长等方面综合考虑进行选择。
对于划分完成之后计算区域中的子计算区域进行排序，每个子计算区域的序号按照0~(N-1)进行编号，每个子计算区域对应分配有线程，该子计算区域中的编号与分配的线程的编号相同，如子计算区域的编号为1，1号子计算区域对应分配的线程编号也为1，这样设置避免了设置多种编号，增加模拟方法的复杂性，两者共用一个编号，使得子计算区域与线程一一对应，方便后续在数值模拟的过程中两者之间能够快速进行配对。
当给每个子计算区域编号完成、线程分配完成后，每个子计算区域对应分配的线程读取该子计算区域中的粒子的信息，粒子的信息主要包括：位移、属性（气、液、固等属性）、速度、密度、光滑长度等信息，可以根据实际需求来获取粒子的相关信息。除此之外，位于每个子计算区域中的线程还需要读取该子计算区域中计算时间步长、计算总时间步数、模块设置等计算参数信息，还包括线程属性信息，线程对应的子计算区域中的边界信息，如x min ，x max ，y min ，y max ，z min ，z max ，子计算区域的以及与其相邻的子计算区域的编号信息、
Figure 894972DEST_PATH_IMAGE026
信息等，其中，线程获取相邻线程的
Figure 588121DEST_PATH_IMAGE026
信息是用于第一个时间步时设置联结区宽度。线程获取子计算区域以及与其相邻的子计算区域的编号信息，是为了方便在后续的计算过程中，可以及时快速的找到与之相邻的子计算区域及该区域对应分配的线程，两个线程之间进行信息的交换更新等操作，需要说明的是，线程不需要获取与该子计算区域不相邻的其他子计算区域的信息，避免信息混乱，增加数据量。
如图2所示，粒子处于二维空间中，将计算区域中的粒子按照二维坐标系x-y中x方向将粒子均分分割，沿y方向将粒子均分分割，分为12个子计算区域，每个子计算区域进行编号，每个子计算区域分配有对应的线程，每个子计算区域中的粒子数量相等，均为4个。
拉格朗日型粒子类数值模拟方法大多数基于核函数积分进行计算，如光滑粒子流体动力学方法：
Figure 321722DEST_PATH_IMAGE027
其中
Figure 74915DEST_PATH_IMAGE028
表示函数变量，x，
Figure 221862DEST_PATH_IMAGE029
表示不同位置的坐标信息，
Figure 921965DEST_PATH_IMAGE030
为包括x的积分体积，W表示核函数，h表示粒子光滑长度，每个线程在对应计算区域边界处的粒子需要依赖相邻子计算区域中边界处的粒子进行积分插值计算，由此可以看出，每个线程需要对对应的子计算区域中的粒子进行计算时，都需要获取相邻子计算区域划分时的边界处的粒子的信息，因此，划分之后的子计算区域需要设置划分边界处的联结区，来保证子计算区域边界处粒子的计算精度，如图3所示，在二维空间的水平方向上，第n个子计算区域在设置联结区的时候需要设置两个，与n-1靠近的左联结区，与n+1靠近的右联结区，说明：二维空间中本来需要4个联结区，水平方向2个，垂直方向2个，这里为简单示意，只画了水平方向的2个，三维空间需要6个联结区。
在对划分之后的子计算区域进行联结区设置之后，为了子计算区域的粒子的数值计算，需要将当前时刻与其相邻的子计算区域联结区的粒子信息进行获取，并对下一时刻各子计算区域中的粒子信息进行更新，如当前时刻中，各个子计算区域对应线程对与其相邻的子计算区域中处于联结区中的粒子的信息进行获取，如图4所示，水平方向需要获取第n-1个子计算区域中的右联结区中粒子的信息，以及第n+1个子计算区域中的左联结区中的粒子的信息，两个相邻联结区中的粒子只是进行信息的获取，为了进行后续的计算，并不是将相邻的联结区的粒子进行交换。在获取到相邻子计算区域的粒子信息并计算当前时刻各子计算区域中粒子的信息的变化率，并根据当前的变化率来获得下一时刻各子计算区域中粒子的信息，然后判断下一时刻各子计算区域中粒子是否越过了当前各子计算区域的边界，当没有越过时，继续执行后面的步骤，当越过时，将对应的越过的粒子发送至越界后所对应的子计算区域中，再进行后续的计算。
然后再判断是否到达计算机中读取的计算时长，若没有到达，则重复步骤S200-S400，直到到达计算时长，结束计算。
需要说明的是，本实施例中的线程指的是计算机中的核处理器，也可以是其他种类的线程，在此不做限制。
需要说明的是，本实施例提供的方法还适用于三维空间的拉格朗日型粒子类数值模拟方法，本实施例中仅是为了方便理解，采用二维空间进行举例。
因此，本发明提供的一种适用于拉格朗日型粒子类数值模拟的并行方法，通过将计算区域划分为多个子计算区域，每个子计算区域中配置有对应的线程，各线程同时运行进行粒子的数值计算，并且在每个子计算区域中设置了联结区，该方法可以解决大规模，如亿量级粒子规模问题计算，提高了粒子类数值模拟算法的计算效率和内存供给量，且根据粒子类数值模拟方法的特殊性，通过设置的联结区提高了粒子信息在计算过程中传递的完整度，很大程度上提高了粒子类数值方法的工程应用能力，如提升了粒子类算法在处理航空、航天、航海、预防自然灾害等领域的工程应用能力。
进一步地，步骤S200中：每个子计算区域的联结区的宽度B不小于阈值宽度R。
上述方案中，在针对每个子计算区域在设置联结区时，为了保证粒子插值子计算区域不被联结区边界截断，联结区的宽度B不能小于阈值宽度R，所述阈值宽度和核函数支持域半径相关，连接区的宽度指的是沿x方向的长度，如图3所示，
本实施例中以光滑粒子流体动力学方法（SPH）来进行示例：
进一步地，宽度阈值R的计算方法如下：
Figure 572389DEST_PATH_IMAGE001
其中，k为预设常数，
Figure 496483DEST_PATH_IMAGE002
表示两个相邻子计算区域的联结区的宽度值，n表示第n个子计算区域，m表示第m个子计算区域，第n个子计算区域与第m个子计算区域相邻，
Figure 68410DEST_PATH_IMAGE003
，n=0~N-1，m=0~N-1，N表示子计算区域的总数。
上述方案中，其中，k的取值和选择的核函数相关，如当采用三次样条核函数时，k=2。
进一步地，
Figure 103362DEST_PATH_IMAGE002
的计算方法如下：
Figure 139451DEST_PATH_IMAGE004
其中，
Figure 996348DEST_PATH_IMAGE005
表示第n个子计算区域中所有粒子中的最大光滑长度，
Figure 117887DEST_PATH_IMAGE006
表示第n个子计算区域中的第i个粒子的光滑长度，i表示第n个子计算区域中的第i个粒子；
Figure 956530DEST_PATH_IMAGE007
表示第个子计算区域中所有粒子中的最大光滑长度，
Figure 519230DEST_PATH_IMAGE008
表示第m个子计算区域中的第j个粒子的光滑长度，j表示第个子计算区域中的第j个粒子。
上述方案中，为了使得第n个子计算区域中粒子的光滑长度更接近实际时，使得计算结果的精度更高，
Figure 785126DEST_PATH_IMAGE002
取的是第n个子计算区域中所有粒子中最大的光滑长度和相邻的第m个子计算区域中所有粒子中最大的光滑长度的均值，选取每个子计算区域中粒子光滑长度的最大值的原因是使得联结区的宽度不小于该子区域内所有粒子核函数所需的支持域大小，充分保证各粒子插值计算时所需的周边粒子信息不被联结区边界截断。
进一步地，当计算区域中的粒子直径只随空间改变，不随时间改变时，
Figure 128383DEST_PATH_IMAGE006
和
Figure 973979DEST_PATH_IMAGE008
的计算方法如下：
Figure 719081DEST_PATH_IMAGE009
Figure 155879DEST_PATH_IMAGE010
其中，
Figure 924115DEST_PATH_IMAGE011
表示特定系数，d i 表示第i个粒子的直径，d j 表示第j个粒子的直径。
上述方案中，
Figure 370139DEST_PATH_IMAGE011
的取值可选地为1.1~1.3之间，在大多数情况下，每个粒子的直径是不随时间进行变化，在此种情况下，粒子的光滑长度可以通过上述的方法计算获得，不同子计算区域中的每个粒子的光滑长度不同。
进一步地，
Figure 704169DEST_PATH_IMAGE012
和
Figure 515130DEST_PATH_IMAGE013
的计算方法如下：
Figure 832979DEST_PATH_IMAGE014
或
Figure 754798DEST_PATH_IMAGE015
其中，
Figure 474493DEST_PATH_IMAGE016
表示第n个子计算区域中第i个粒子的面积，
Figure 253093DEST_PATH_IMAGE017
表示第n个子计算区域中第i个粒子的体积；
Figure 995921DEST_PATH_IMAGE018
或
Figure 783748DEST_PATH_IMAGE019
其中，
Figure 295632DEST_PATH_IMAGE020
表示第m个子计算区域中第j个粒子的面积，
Figure 245134DEST_PATH_IMAGE021
表示第m个子计算区域中第j个粒子的体积。
上述方案中，在二维情况下，粒子的直径可以通过粒子所占的面积求平方根进行计算，在三维情况下，粒子的直径可以通过粒子所占的体积求立方根进行计算。
进一步地，当计算区域中的粒子直径随时间变化时，
Figure 271996DEST_PATH_IMAGE006
和
Figure 801197DEST_PATH_IMAGE008
随时间的更新方法如下：
Figure 964325DEST_PATH_IMAGE022
Figure 350307DEST_PATH_IMAGE023
其中，t为计算时长，D为整个模拟对象的维度，D=1，2，3，
Figure 67727DEST_PATH_IMAGE024
表示第n个子计算区域中的第i个粒子的密度，
Figure 197357DEST_PATH_IMAGE025
表示第m个子计算区域中的第j个粒子的密度。
Figure 214992DEST_PATH_IMAGE006
和
Figure 975138DEST_PATH_IMAGE008
可以根据更新率来获得当前时刻对应的粒子直径。
上述方案中，在大多数情况下，每个粒子的直径是不随时间进行变化，但是对于某些粒子在受到局部剧烈膨胀或者压缩问题时，粒子的密度会随之发生剧烈变化，导致粒子光滑长度不光在空间分布上不均匀，也会随着时间进行变化，此时可以通过将粒子的光滑长度随时间的变化率和密度变化率联合起来获得，更好地保持全局的一致性。
进一步地，步骤S100中计算区域划分为多个子计算区域的步骤如下：
获取计算区域所在的坐标系x-y-z;
分别遍历坐标系中x、y、z方向上粒子，在每个1/X、1/Y、1/Z粒子数所对应的x轴位置进行切割，其中X为x方向上预设的计算区域切割数量，Y为y方向上预设的计算区域切割数量，Z为z方向上预设的计算区域切割数量。
上述方案中，在对计算区域进行划分时，现有技术中，通常是采用简单的空间等分的方式，这样分割后，每个子计算区域中的粒子数量可能相差较大，导致每个子计算区域中的线程需要计算的粒子数量不同，会出现有的子计算区域数值计算快，有的慢，导致不能同时计算完，共同进行下一步数值模拟中的计算，而本实施例中，采用的是将粒子等分的方式进行划分的，先获取x方向上粒子需要切割的数量X，然后遍历x方向上的粒子，在每个1/X粒子数量对应的x轴位置进行划分，然后获取y方向上的粒子需要切割的数量Y，遍历y方向的粒子，在每个1/Y粒子数量对应的y轴位置进行划分，当粒子处于三维空间中时，同理对Z轴进行划分。
进一步地，步骤S300中对当前线程中的粒子信息进行更新的方法如下：
基于粒子类数值算法计算当前时刻下各子计算区域粒子信息的变化率；
根据当前时刻下各子计算区域粒子信息的变化率，更新获得下一时刻的各子计算区域粒子信息。
上述方案中，各线程基于粒子类数值算法计算获得t时刻下粒子的密度、速度、能量等参数的变化率，计算方法如下所示：
Figure 976592DEST_PATH_IMAGE031
其中，
Figure 644333DEST_PATH_IMAGE032
为密度，v为速度，e为内能，
Figure 516474DEST_PATH_IMAGE033
为总应力张量，希腊字母α和β代表坐标方向。
然后基于上述的参数变化率更新
Figure 447521DEST_PATH_IMAGE034
时刻下的粒子的位移、密度、速度、能量等参数信息，从相邻线程联结区所接收粒子虽然参与所有计算流程，但接收得到的这些粒子参数信息不在此线程内更新，计算结束前这些粒子信息从线程中去除。
拉格朗日型粒子类数值模拟计算时，粒子的功能除了表示物质粒子信息，还用于场函数的近似，以光滑粒子流体动力学方法为例：
Figure 936271DEST_PATH_IMAGE035
其中，
Figure 345387DEST_PATH_IMAGE036
表示对粒子i支持域内的M个粒子的参数进行叠加求和操作，v ij 表示粒子i和粒子j的速度差。
Figure 337614DEST_PATH_IMAGE037
表示粒子j的质量，W表示核函数。
进一步地，步骤S400中，下一时刻的各子计算区域粒子越过当前各子计算区域边界的判断方法如下：
获取下一时刻各子计算区域更新后粒子的坐标信息；
获取当前时刻各子计算区域边界的坐标信息；
当粒子的坐标信息值大于边界的坐标信息时，则该粒子越过当前子计算区域的边界。
上述方案中，根据
Figure 501879DEST_PATH_IMAGE034
时刻更新得到的各子计算区域更新得到的粒子i坐标信息，对比对应子计算区域边界的坐标，判断更新后的粒子是否越过该子计算区域的边界，判断越界的方法如下：
Figure 150029DEST_PATH_IMAGE038
上式中x i ，y i ，z i ，代表第i个粒子在x， y， z方向上的坐标，x min ，x max ，y min ，y max ，z min ， z max 等代表当前子计算区域边界坐标。
将越过子计算区域的边界处的粒子信息打包发送到粒子现在所在子计算区域以及线程中，并和所在子计算区域中原有的粒子进行递补排序。如果粒子没有越过子计算区域的边界时，则下一时刻依然由原有线程进行计算。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种适用于拉格朗日型粒子类数值模拟的并行方法，其特征在于，包括如下步骤：
步骤S100：获取计算区域内的粒子，将计算区域划分为多个子计算区域，每个子计算区域内的粒子数量相同；每个子计算区域对应的线程获取当前子计算区域中的粒子信息、计算参数信息，以及线程属性信息；
步骤S200：设置每个子计算区域的联结区；
步骤S300：获取当前时刻与各子计算区域相邻子计算区域的联结区内粒子信息，并对下一时刻各子计算区域中的粒子信息进行更新；
步骤S400：当下一时刻的流体力学领域各子计算区域粒子越过当前各子计算区域的边界时，将更新后下一时刻各子计算区域中越界的粒子发送至对应的子计算区域中；
步骤S200中：
每个子计算区域的联结区的宽度B不小于阈值宽度R；
宽度阈值R的计算方法如下：
R＝k×rn，m
其中，k为预设常数，rn,m表示两个相邻子计算区域的联结区的宽度值，n表示第n个子计算区域，m表示第m个子计算区域，第n个子计算区域与第m个子计算区域相邻，n≠m，n＝0～N-1，m＝0～N-1，N表示子计算区域的总数；
rn,m的计算方法如下：
Figure QLYQS_1
其中，max(hi)n表示第n个子计算区域中所有粒子中的最大光滑长度，(hi)n表示第n个子计算区域中的第i个粒子的光滑长度，i表示第n个子计算区域中的第i个粒子；max(hj)m表示第m个子计算区域中所有粒子中的最大光滑长度，(hj)m表示第m个子计算区域中的第j个粒子的光滑长度，j表示第m个子计算区域中的第j个粒子。
2.如权利要求1所述的适用于拉格朗日型粒子类数值模拟的并行方法，其特征在于，当计算区域中的粒子直径只随空间改变，不随时间改变时，(hi)n和(hj)m的计算方法如下：
(hi)n＝ε(di)n
(hj)m＝ε(dj)m
其中，ε表示特定系数，di表示第i个粒子的直径，dj表示第j个粒子的直径。
3.如权利要求2所述的适用于拉格朗日型粒子类数值模拟的并行方法，其特征在于，(di)n和(dj)m的计算方法如下：
Figure QLYQS_2
或
Figure QLYQS_3
其中，(Si)n表示第n个子计算区域中第i个粒子的面积，(Vi)n表示第n个子计算区域中第i个粒子的体积；
Figure QLYQS_4
或
Figure QLYQS_5
其中，(Sj)m表示第m个子计算区域中第j个粒子的面积，(Vj)m表示第m个子计算区域中第j个粒子的体积。
4.如权利要求1所述的适用于拉格朗日型粒子类数值模拟的并行方法，其特征在于，当计算区域中的粒子直径随时间变化时，(hi)n和(hj)m随时间的更新方法如下：
Figure QLYQS_6
Figure QLYQS_7
其中，t为计算时长，D为整个模拟对象的维度，D＝1，2，3，(ρi)n表示第n个子计算区域中的第i个粒子的密度，(ρj)m表示第m个子计算区域中的第j个粒子的密度。
5.如权利要求1-4之一所述的适用于拉格朗日型粒子类数值模拟的并行方法，其特征在于，步骤S100中计算区域划分为多个子计算区域的步骤如下：
获取计算区域所在的坐标系x-y-z；
分别遍历坐标系中x、y、z方向上粒子，在每个1/X、1/Y、1/Z粒子数所对应的x轴位置进行切割，其中X为x方向上预设的计算区域切割数量，Y为y方向上预设的计算区域切割数量，Z为z方向上预设的计算区域切割数量。
6.如权利要求1-4之一所述的适用于拉格朗日型粒子类数值模拟的并行方法，其特征在于，步骤S300中对当前线程中的粒子信息进行更新的方法如下：
基于粒子类数值算法计算当前时刻下各子计算区域粒子信息的变化率；
根据当前时刻下各子计算区域粒子信息的变化率，更新获得下一时刻的各子计算区域粒子信息。
7.如权利要求1-4之一所述的适用于拉格朗日型粒子类数值模拟的并行方法，其特征在于，步骤S400中，下一时刻的各子计算区域粒子越过当前各子计算区域的边界的判断方法如下：
获取下一时刻各子计算区域更新后粒子的坐标信息；
获取当前时刻各子计算区域边界的坐标信息；
当粒子的坐标信息值大于边界的坐标信息时，则该粒子越过当前子计算区域的边界。
CN202210014743.9A 2022-01-07 2022-01-07 一种适用于拉格朗日型粒子类数值模拟的并行方法 Active CN114357907B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210014743.9A CN114357907B (zh) 2022-01-07 2022-01-07 一种适用于拉格朗日型粒子类数值模拟的并行方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210014743.9A CN114357907B (zh) 2022-01-07 2022-01-07 一种适用于拉格朗日型粒子类数值模拟的并行方法

Publications (2)

Publication Number Publication Date
CN114357907A CN114357907A (zh) 2022-04-15
CN114357907B true CN114357907B (zh) 2023-03-21

Family

ID=81106805

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210014743.9A Active CN114357907B (zh) 2022-01-07 2022-01-07 一种适用于拉格朗日型粒子类数值模拟的并行方法

Country Status (1)

Country Link
CN (1) CN114357907B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111241742A (zh) * 2019-12-27 2020-06-05 西安交通大学 一种多相流计算方法
CN112989683A (zh) * 2021-04-19 2021-06-18 中国人民解放军国防科技大学 一种sph的向量化并行计算方法及装置

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2014035614A (ja) * 2012-08-08 2014-02-24 Furukawa Electric Co Ltd:The 粒子成長過程のシミュレーション方法、プログラム及びシミュレーション装置
CN107633123B (zh) * 2017-09-13 2021-05-18 浙江工业大学 一种用于光滑粒子流体动力学模拟出血及处理加速的方法
CN109948109A (zh) * 2019-01-31 2019-06-28 天津大学 含截面变化的明渠非恒定流无网格粒子模拟方法
CN110929456B (zh) * 2019-11-13 2021-07-06 西安交通大学 移动粒子法并行计算等效粒子负载均衡加速方法
CN112733415B (zh) * 2021-01-14 2022-07-01 中国海洋大学 一种薄壁弹性体边界的无网格处理方法、装置、终端设备及计算介质
CN113850032B (zh) * 2021-12-02 2022-02-08 中国空气动力研究与发展中心计算空气动力研究所 一种数值模拟计算中的负载均衡方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111241742A (zh) * 2019-12-27 2020-06-05 西安交通大学 一种多相流计算方法
CN112989683A (zh) * 2021-04-19 2021-06-18 中国人民解放军国防科技大学 一种sph的向量化并行计算方法及装置

Also Published As

Publication number Publication date
CN114357907A (zh) 2022-04-15

Similar Documents

Publication Publication Date Title
CN109063275B (zh) 基于feap的三维多晶微观结构材料模型的构建方法
Li et al. Arc–surface intersection method to calculate cutter–workpiece engagements for generic cutter in five-axis milling
CN102681972A (zh) 一种利用GPU加速格子-Boltzmann的方法
CN104933225B (zh) 实现计算流体力学大规模实时模拟的方法
WO2011122423A1 (en) System and method for optimizing machining simulation
JPH0816629A (ja) 解析用メッシュ作成方法及び装置
US20190176405A1 (en) Computer aided design with high resolution lattice structures using graphics processing units (gpu)
Chao et al. Improved hybrid bounding box collision detection algorithm
Yu et al. A robust Delaunay-AFT based parallel method for the generation of large-scale fully constrained meshes
CN114357907B (zh) 一种适用于拉格朗日型粒子类数值模拟的并行方法
Jones et al. Visualizing flow trajectories using locality-based rendering and warped curve plots
CN1395196A (zh) 过程冶金中计算流体力学的数据可视化方法
Eddy et al. Multidimensional design visualization in multiobjective optimization
Zhang et al. A GPU-based parallel slicer for 3D printing
CN117115393A (zh) 一种基于gpu的nurbs曲面并行求交方法、设备及存储介质
Gidaspov et al. A software package for simulation of unsteady flows of the reacting gas in the channel
JP5429214B2 (ja) プログラム作成支援装置、プログラム作成支援方法およびコンピュータプログラム
Huang et al. Parallel Performance and Optimization of the Lattice Boltzmann Method Software Palabos Using CUDA
US10650172B1 (en) Systems and methods for implementing iterative simulation manipulations and results display
JP4622987B2 (ja) 工具参照面データの作成装置と作成方法
Austern et al. Rationalization and Optimization of Concrete Façade Panels
JP2005078416A (ja) 解析モデル生成方法および装置ならびにプログラムおよびその記憶媒体
JPH09305651A (ja) 解析シミュレーション装置とその高速表示方法
Zhu et al. Realtime simulation of burning solids on GPU with CUDA
WO2019026357A1 (ja) 情報処理装置、情報処理方法、情報処理システム、プログラム製造方法及びプログラム

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