一种确定精确格林函数的方法及装置与流程

本发明属于声学
技术领域：
，更具体的涉及一种确定精确格林函数的方法及装置。
背景技术：
：在声学领域中，格林函数反映了声场的基本结构，是单位强度点声源所辐射的声场。如果声源和格林函数已知，就可获取声场的解。比如，对于空间中一个位于y点的强度为q、圆频率为ω的单极子点声源q(y,ω)，在空间中观察点x点处产生的声压p(x,ω)可以通过下列公式确定：p(x,ω)＝q(y,ω)·G(x,y,ω)其中，G(x,y,ω)为关联声源y和观察点x的格林函数。在实际应用中，对于无界空间，可直接采用具有解析解的自由空间格林函数G0(x,y,ω)进行声场的计算；然而对于有界空间，因声波与边界相互作用的原因，点声源产生的声场远较无界空间复杂。仅当边界形状及声学性质相对简单的少数情形，才有可能获得格林函数的解析解。如果声波传播空间存在复杂的边界条件，比如边界几何形状不规则且有阻抗特性，则需要采用数值方法求解满足相应几何边界条件和声学边界条件的格林函数，该格林函数称为精确格林函数。若边界的几何尺寸远小于声波的波长，则该边界是紧致的，边界对声波的散射作用可以忽略，声波相当于在自由空间内传播，可以采用自由空间格林函数求解声场；反之，若边界的几何尺寸大于或者接近于声波的波长，则边界是非紧致的，声源激发的噪声在向空间辐射的同时，会在边界发生散射，使得噪声的大小和分布与自由空间下的声场存在很大的不同。另一方面，声波传播对不同性质的边界需要满足特定的声学边界条件，当边界为刚性边界时，入射的声波会在边界上发生全反射，称为声学硬边界；如果边界不是刚性的，比如边界材料为阻尼材料，此时的边界属于阻抗边界，入射的声波不仅会被边界反射，还可能被边界吸收。一般的数学物理问题，其边界条件可以概括为Dirichlet条件、Neumann条件和Robin条件三种，声学硬边界和阻抗边界分别对应于Neumann边界条件和Robin边界条件。近年来，大量学者对精确格林函数开展了研究工作，比如，有人提出一种带高阶正交多项式的谱配置边界元方法来计算精确格林函数有人提出提出一种基于边界元方法的精确格林函数数值计算方法；还有人利用光学中的惠更斯原理，提出了一种求解精确格林函数的级数叠代方法。综上所述，现有技术中，精确格林函数的计算方法都假定边界为声学硬边界，忽略了边界阻抗特性对声传播的影响，从而导致边界阻抗特性在声学分析中的应用受到限定的问题。技术实现要素：本发明实施例提供一种确定精确格林函数的方法及装置，用以解决现有技术中，精确格林函数的计算方法都假定边界为声学硬边界，忽略了边界阻抗特性对声传播的影响，从而导致边界阻抗特性在声学分析中的应用受到限定的问题。本发明实施例提供一种确定精确格林函数的方法，包括：建立边界元模型，将所述边界元模型S离散成M个网格单元，依次确定所述M个网格单元的中心点，单位外法线方向以及面积；在所述M个网格单元中确定任意两个中心点为zm和zn的网格单元，根据所述zm和所述zn点所在散射边界的声学特性以及公式(1)，确定所述zm和所述zn点处的声学参数β(zm)和β(zn)，确定关联所述zm点和所述zn点的自由空间格林函数G0(zm,zn,ω)以及所述自由空间格林函数G0(zm,zn,ω)沿所述zm点外法线方向的偏导数G0(zm,zn,ω)/n(zm)；在所述边界元模型S外确定观察点x，根据所述zn点和所述观察点x，确定关联所述观察点x和所述zn点之间的自由空间格林函数G0(x,zn,ω)；根据所述G0(zm,zn,ω)、所述G0(zm,zn,ω)/n(zm)和所述G0(x,zn,ω)，通过公式(2)确定关联所述观察点x和所述zn点之间的格林函数G(x,zn,ω)；在所述边界元模型S外确定声源点y，根据所述观察点x和所述声源点y，确定关联所述观察点x和所述声源点y之间的自由空间格林函数G0(x,y,ω)；根据所述zn点和所述声源点y，确定关联所述zn点和所述声源点y之间的自由空间格林函数G0(zn,y,ω)及所述自由空间格林函数G0(zn,y,ω)沿zn点外法线方向的偏导数G0(zn,y,ω)/n(zn)；根据所述G(x,zn,ω)、所述G0(zn,y,ω)和所述G0(zn,y,ω)/n(zn)通过公式(3)确定散射格林函数GS(x,y,ω)；根据所述G0(x,y,ω)和所述GS(x,y,ω)，通过公式(4)确定精确格林函数G(x,y,ω)；其中，公式(1)如下所示：β(zm)＝αjρ0ω/Z(zm,ω)，β(zn)＝αjρ0ω/Z(zn,ω)公式(2)如下所示：12G(x,zn,ω)=G0(x,zn,ω)+∫S′G(x,zm,ω)∂G0(zm,zn,ω)∂n(zm)+β(zm)G0(zm,zn,ω)ds(zm)]]>公式(3)如下所示：GS(x,y,ω)=∫SG(x,zn,ω)∂G0(zn,y,ω)∂n(zn)+β(zn)G0(zn,y,ω)ds(zn)]]>公式(4)如下所示：G(x,y,ω)＝G0(x,y,ω)+GS(x,y,ω)其中，ω表示圆频率，j为虚数单位，ρ0为声传播介质的密度，G0是频域自由空间格林函数；x和y分别为观察点和声源点的位置矢量，zm和zn是边界上网格单元的中心点；n(zm)和n(zn)分别表示zm和zn处边界的外法线方向，s(zm)和s(zn)分别表示zm和zn点所在的网格单元面积，S是边界元模型，S'是除去zm点所在网格单元的散射边界；Z(zm)和Z(zn)分别表示zm和zn点所在的网格单元声学阻抗，α为常数，对于声学硬边界α＝0，对于声学阻抗边界α＝1。本发明实施例还提供一种确定精确格林函数的装置，包括：建立单元，用于建立边界元模型，将所述边界元模型S离散成M个网格单元，依次确定所述M个网格单元的中心点，单位外法线方向以及面积；第一确定单元，用于在所述M个网格单元中确定任意两个中心点为zm和zn的网格单元，根据所述zm和所述zn点所在散射边界的声学特性以及公式(1)，确定所述zm和所述zn点处的声学参数β(zm)和β(zn)，确定关联所述zm点和所述zn点的自由空间格林函数G0(zm,zn,ω)以及所述自由空间格林函数G0(zm,zn,ω)沿所述zm点外法线方向的偏导数G0(zm,zn,ω)/n(zm)；第二确定单元，用于在所述边界元模型S外确定观察点x，根据所述zn点和所述观察点x，确定关联所述观察点x和所述zn点之间的自由空间格林函数G0(x,zn,ω)；第三确定单元，用于根据所述G0(zm,zn,ω)、所述G0(zm,zn,ω)/n(zm)和所述G0(x,zn,ω)，通过公式(2)确定关联所述观察点x和所述zn点之间的格林函数G(x,zn,ω)；第四确定单元，用于在所述边界元模型S外确定声源点y，根据所述观察点x和所述声源点y，确定关联所述观察点x和所述声源点y之间的自由空间格林函数G0(x,y,ω)；根据所述zn点和所述声源点y，确定关联所述zn点和所述声源点y之间的自由空间格林函数G0(zn,y,ω)及所述自由空间格林函数G0(zn,y,ω)沿zn点外法线方向的偏导数G0(zn,y,ω)/n(zn)；第五确定单元，用于根据所述G(x,zn,ω)、所述G0(zn,y,ω)和所述G0(zn,y,ω)/n(zn)通过公式(3)确定散射格林函数GS(x,y,ω)；第六确定单元，用于根根据所述G0(x,y,ω)和所述GS(x,y,ω)，通过公式(4)确定精确格林函数G(x,y,ω)；其中，公式(1)如下所示：β(zm)＝αjρ0ω/Z(zm,ω)，β(zn)＝αjρ0ω/Z(zn,ω)公式(2)如下所示：12G(x,zn,ω)=G0(x,zn,ω)+∫S′G(x,zm,ω)∂G0(zm,zn,ω)∂n(zm)+β(zm)G0(zm,zn,ω)ds(zm)]]>公式(3)如下所示：GS(x,y,ω)=∫SG(x,zn,ω)∂G0(zn,y,ω)∂n(zn)+β(zn)G0(zn,y,ω)ds(zn)]]>公式(4)如下所示：G(x,y,ω)＝G0(x,y,ω)+GS(x,y,ω)其中，ω表示圆频率，j为虚数单位，ρ0为声传播介质的密度，G0是频域自由空间格林函数；x和y分别为观察点和声源点的位置矢量，zm和zn是边界上网格单元的中心点；n(zm)和n(zn)分别表示zm和zn处边界的外法线方向，s(zm)和s(zn)分别表示zm和zn点所在的网格单元面积，S是边界元模型，S'是除去zm点所在网格单元的散射边界；Z(zm)和Z(zn)分别表示zm和zn点所在的网格单元声学阻抗，α为常数，对于声学硬边界α＝0，对于声学阻抗边界α＝1。本发明实施例中，提供了一种确定精确格林函数的方法及装置，包括：建立边界元模型，将所述边界元模型S离散成M个网格单元，依次确定所述M个网格单元的中心点，单位外法线方向以及面积；在所述M个网格单元中确定任意两个中心点为zm和zn的网格单元，根据所述zm和所述zn点所在散射边界的声学特性以及β(zm)＝αjρ0ω/Z(zm,ω)、β(zn)＝αjρ0ω/Z(zn,ω)，确定所述zm和所述zn点处的声学参数β(zm)和β(zn)，确定关联所述zm点和所述zn点的自由空间格林函数G0(zm,zn,ω)以及所述自由空间格林函数G0(zm,zn,ω)沿所述zm点外法线方向的偏导数G0(zm,zn,ω)/n(zm)；在所述边界元模型S外确定观察点x，根据所述zn点和所述观察点x，确定关联所述观察点x和所述zn点之间的自由空间格林函数G0(x,zn,ω)；根据所述G0(zm,zn,ω)、所述G0(zm,zn,ω)/n(zm)和所述G0(x,zn,ω)，通过确定关联所述观察点x和所述zn点之间的格林函数G(x,zn,ω)；在所述边界元模型S外确定声源点y，根据所述观察点x和所述声源点y，确定关联所述观察点x和所述声源点y之间的自由空间格林函数G0(x,y,ω)；根据所述zn点和所述声源点y，确定关联所述zn点和所述声源点y之间的自由空间格林函数G0(zn,y,ω)及所述格林函数沿zn点外法线方向的偏导数G0(zn,y,ω)/n(zn)；根据所述G(x,zn,ω)、所述G0(zn,y,ω)和所述G0(zn,y,ω)/n(zn)通过确定散射格林函数GS(x,y,ω)；根据所述G0(x,y,ω)和所述GS(x,y,ω)，通过G(x,y,ω)＝G0(x,y,ω)+GS(x,y,ω)确定精确格林函数G(x,y,ω)；其中，ω表示圆频率，j为虚数单位，ρ0为声传播介质的密度，G0是频域自由空间格林函数；x和y分别为观察点和声源点的位置矢量，zm和zn是边界上网格单元的中心点；n(zm)和n(zn)分别表示zm和zn处边界的外法线方向，s(zm)和s(zn)分别表示zm和zn点所在的网格单元面积，S是边界元模型，S'是除去zm点所在网格单元的散射边界；Z(zm)和Z(zn)分别表示zm和zn点所在的网格单元声学阻抗，α为常数，对于声学硬边界α＝0，对于声学阻抗边界α＝1。本发明实施例所提供的精确格林函数计算方法，采用了Robin边界条件，不仅能分析声学非紧致硬界的散射，还能考虑非紧致阻抗边界对声传播的影响，是一种应用范围更广的精确格林函数计算方法。附图说明为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。图1为本发明实施例提供的一种确定精确格林函数的方法流程示意图；图2为本发明实施例一提供的二维非紧致圆柱声散射示意图；图3为本发明实施例一提供的圆柱声学硬边界条件下精确格林函数数值解、解析解与自由空间格林函数的对比图；图4为本发明实施例一提供的圆柱声学阻抗边界条件下精确格林函数数值解、解析解与声学硬边界条件下精确格林函数数值解的对比图；图5为本发明实施例提供的一种确定精确格林函数的装置结构示意图。图中符号说明如下：1-非紧致圆柱；2-观察点；3-声源点；4-圆柱表面点；5-自由空间格林函数曲线；6-圆柱声学硬边界条件下精确格林函数解析解曲线；7-圆柱声学硬边界条件下精确格林函数数值解曲线；8-圆柱声学硬边界条件下精确格林函数数值解曲线；9-圆柱声学阻抗边界条件下精确格林函数解析解曲线；10-圆柱声学阻抗边界条件下精确格林函数数值解曲线。具体实施方式下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。图1示例性的示出了本发明实施例提供的一种确定精确格林函数的方法流程示意图，该方法至少可以应用于声学分析中。如图1所示，为本发明实施例提供的一种确定精确格林函数的方法，包括以下步骤：步骤101，建立边界元模型，将所述边界元模型S离散成M个网格单元，依次确定所述M个网格单元的中心点，单位外法线方向以及面积；步骤102，在所述M个网格单元中确定任意两个中心点为zm和zn的网格单元，根据所述zm和所述zn点所在散射边界的声学特性以及公式(1)，确定所述zm和所述zn点处的声学参数β(zm)和β(zn)，确定关联所述zm点和所述zn点的自由空间格林函数G0(zm,zn,ω)以及所述自由空间格林函数G0(zm,zn,ω)沿所述zm点外法线方向的偏导数G0(zm,zn,ω)/n(zm)；步骤103，在所述边界元模型S外确定观察点x，根据所述zn点和所述观察点x，确定关联所述观察点x和所述zn点之间的自由空间格林函数G0(x,zn,ω)；步骤104，根据所述G0(zm,zn,ω)、所述G0(zm,zn,ω)/n(zm)和所述G0(x,zn,ω)，通过公式(2)确定关联所述观察点x和所述zn点之间的格林函数G(x,zn,ω)；步骤105，在所述边界元模型S外确定声源点y，根据所述观察点x和所述声源点y，确定关联所述观察点x和所述声源点y之间的自由空间格林函数G0(x,y,ω)；根据所述zn点和所述声源点y，确定关联所述zn点和所述声源点y之间的自由空间格林函数G0(zn,y,ω)及所述格林函数沿zn点外法线方向的偏导数G0(zn,y,ω)/n(zn)；步骤106，根据所述G(x,zn,ω)、所述G0(zn,y,ω)和所述G0(zn,y,ω)/n(zn)通过公式(3)确定散射格林函数GS(x,y,ω)；步骤107，根据所述G0(x,y,ω)和所述GS(x,y,ω)，通过公式(4)确定精确格林函数G(x,y,ω)；其中，公式(1)如下所示：β(zm)＝αjρ0ω/Z(zm,ω)，β(zn)＝αjρ0ω/Z(zn,ω)公式(2)如下所示：12G(x,zn,ω)=G0(x,zn,ω)+∫S′G(x,zm,ω)∂G0(zm,zn,ω)∂n(zm)+β(zm)G0(zm,zn,ω)ds(zm)]]>公式(3)如下所示：GS(x,y,ω)=∫SG(x,zn,ω)∂G0(zn,y,ω)∂n(zn)+β(zn)G0(zn,y,ω)ds(zn)]]>公式(4)如下所示：G(x,y,ω)＝G0(x,y,ω)+GS(x,y,ω)其中，ω表示圆频率，j为虚数单位，ρ0为声传播介质的密度，G0是频域自由空间格林函数；x和y分别为观察点声源点的位置矢量，zm和zn是边界上网格单元的中心点；n(zm)和n(zn)分别表示zm和zn处边界的外法线方向，s(zm)和s(zn)分别表示zm和zn点所在的网格单元面积，S是边界元模型，S'是除去zm点所在网格单元的散射边界；Z(zm)和Z(zn)分别表示zm和zn点所在的网格单元声学阻抗，β(zm)和β(zn)分别表示zm和zn点所在的网格单元的声学边界条件参数，α为常数，对于声学硬边界α＝0，对于声学阻抗边界α＝1。以下结合上述步骤，具体介绍本发明实施例提供一种确定精确格林函数的方法：需要说明的是，在实际应用中，若精确格林函数的计算方法都是假定边界为声学硬边界，而忽略了边界阻抗特性这个因素，在本发明实施例中，在计算精确格林函数之前，需要先确认参数的边界特性。在步骤102中，在M个网格单元中确定任意两个中心点为zm和zn的网格单元，根据zm和zn点所在散射边界的声学特性以及公式(1)，可以确定zm和zn点处的声学参数β(zm)和β(zn)，具体地，如下公式(1)所示：β(zm)＝αjρ0ω/Z(zm,ω)β(zn)＝αjρ0ω/Z(zn,ω)(1)其中，ω表示圆频率，Z(ω)为声学边界的阻抗，j为虚数单位，ρ0为声传播介质的密度，α为常数，对于声学硬边界α＝0，对于声学阻抗边界α＝1。在本发明实施例中，由于边界模型可以为二维空间，也可以是三维空间，即在确定zm点和zn点的自由空间格林函数G0(zm,zn,ω)，需要区分二维空间和三维空间。具体地，当边界模型为二维空间时，zm点和zn点的自由空间格林函数G0(zm,zn,ω)可以通过下列公式(5)确定：G0(zm,zn,ω)=j4H0(1)(k|zm-zn|)---(5)]]>其中，G0(zm,zn,ω)为zm点和zn点的自由空间格林函数，为0阶第一类Hankel函数，j为虚数单位，k＝ω/c0为声学波数，zm和zn分别是编号为m和n的网格单元的中心点。当边界模型为三维空间时，zm点和zn点的自由空间格林函数G0(zm,zn,ω)可以通过下列公式(6)确定：G0(zm,zn,ω)=14π|zm-zn|ejk|zm-zn|---(6)]]>进一步地，关联zm点和zn点的自由空间格林函数G0(zm,zn,ω)沿zm点外法线方向的偏导数G0(zm,zn,ω)/n(zm)的确定，也需要区分二维空间和三维空间，具体地：当所述边界元模型为二维空间时，关联zm点和zn点的自由空间格林函数G0(zm,zn,ω)沿zm点外法线方向的偏导数G0(zm,zn,ω)/n(zm)可以通过公式(7)确定：∂G0(zm,zn,ω)∂n(zm)=-jk4H1(1)(k|zm-zn|)|zm-zn|(zm-zn)·n(zm)---(7)]]>其中，G0(zm,zn,ω)/n(zm)为关联zm点和zn点的自由空间格林函数G0(zm,zn,ω)沿zm点外法线方向的偏导数，为1阶第一类Hankel函数，n(zm)表示zm处边界的外法线方向。当所述边界元模型为三维空间时，关联zm点和zn点的自由空间格林函数G0(zm,zn,ω)沿zm点外法线方向的偏导数G0(zm,zn,ω)/n(zm)可以通过公式(8)确定：∂G0(zm,zn,ω)∂n(zm)=14πejk|zm-zn|(jk|zm-zn|-1)|zm-zn|3(zm-zn)·n(zm)---(8)]]>在步骤103中，在边界元模型S外确定观察点x，根据zn点和观察点x，确定关联观察点x和zn点之间的自由空间格林函数G0(x,zn,ω)。需要说明的是，在确定关联观察点x和zn点之间的自由空间格林函数G0(x,zn,ω)时，也需要区分二维空间和三维空间。具体地，当边界元模型S为二维空间时，关联观察点x和zn点之间的自由空间格林函数G0(x,zn,ω)可以通过下列公式(9)确定：G0(x,zn,ω)=j4H0(1)(k|x-zn|)---(9)]]>其中，G0(x,zn,ω)为关联观察点x和zn点之间的自由空间格林函数。当所述边界元模型S为三维空间时，关联观察点x和zn点之间的自由空间格林函数G0(x,zn,ω)可以通过下列公式(10)确定：G0(x,zn,ω)=14π|x-zn|ejk|x-zn|---(10)]]>在步骤104中，根据所有zn点和G0(x,zn,ω)，通过下列公式(11)确定矩阵G0(x,z,ω)：G0(x,z,ω)=G0(x,z1,ω)G0(x,z2,ω)...G0(x,zM,ω)---(11)]]>根据G0(zm,zn,ω)和G0(zm,zn,ω)/n(zm)，通过下列公式(12)确定H矩阵的子项：Hmn=∫S′∂G0(zm,zn,ω)∂n(zm)+βG0(zm,zn,ω)ds(zm),m≠n0,m=n---(12)]]>需要说明的是，对于二维问题，可以将公式(5)、(7)确定的G0(zm,zn,ω)代入上述公式(12)内，确定Hmn；而对于三维问题，可以将公式(6)、(8)确定的G0(zm,zn,ω)代入上述公式(12)内，确定Hmn。进一步地，根据H矩阵的子项和G0(x,z,ω)矩阵，通过下列公式(13)确定代数离散方程组：|E2+H|G(x,z,ω)=G0(x,z,ω)---(13)]]>其中，E是单位对角矩阵，Hmn为矩阵H下标为m,n的子项，矩阵G(x,z,ω)为对应于所述所有zn点需要求解的未知矩阵，具体如以下公式(14)所示：G(x,z,ω)=G(x,z1,ω)G(x,z2,ω)...G(x,zM,ω)---(14)]]>进一步地，根据上述公式(13)确定的代数离散方程组和G(x,z,ω)，通过公式(2)，可以确定关联观察点x和zn点之间的格林函数G(x,zn,ω)：12G(x,zn,ω)=G0(x,zn,ω)+∫S′G(x,zm,ω)∂G0(zm,zn,ω)∂n(zm)+β(zm)G0(zm,zn,ω)ds(zm)---(2)]]>其中，G(x,zn,ω)为关联观察点x和zn点之间的格林函数。在步骤105中，在边界元模型S外确定声源点y，根据观察点x和声源点y，确定关联观察点x和声源点y之间的自由空间格林函数G0(x,y,ω)时，需要区分二维空间和三维空间。具体地，当所述边界元模型S为二维空间时，通过下列公式(15)确定关联观察点x和声源点y之间的自由空间格林函数G0(x,y,ω)：G0(x,y,ω)=j4H0(1)(k|x-y|)---(15)]]>其中，G0(x,y,ω)为关联观察点x和声源点y之间的自由空间格林函数。当所述边界元模型S为三维空间时，通过下列公式(16)确定关联观察点x和声源点y之间的自由空间格林函数G0(x,y,ω)：G0(x,y,ω)=14π|x-y|ejk|x-y|---(16)]]>进一步地，根据zn点和声源点y，确定关联zn点和声源点y之间的自由空间格林函数G0(zn,y,ω)时，需要区分二维空间和三维空间。具体地，当所述边界元模型S为二维空间时，通过下列公式(17)确定关联zn点和声源点y之间的自由空间格林函数G0(zn,y,ω)：G0(zn,y,ω)=j4H0(1)(k|zn-y|)---(17)]]>其中，G0(zn,y,ω)为关联zn点和声源点y之间的自由空间格林函数。当所述边界元模型S为三维空间时，通过下列公式(18)确定关联zn点和声源点y之间的自由空间格林函数G0(zn,y,ω)：G0(zn,y,ω)=14π|zn-y|ejk|z-y|---(18)]]>进一步地，根据zn点和声源点y，确定关联zn点和声源点y之间的自由空间格林函数G0(zn,y,ω)沿zn点外法线方向的偏导数G0(zn,y,ω)/n(zn)，也需要区分二维空间和三维空间。具体地，当所述边界元模型S为二维空间时，通过下列公式(19)关联zn点和声源点y之间的自由空间格林函数G0(zn,y,ω)沿zn点外法线方向的偏导数G0(zn,y,ω)/n(zn)：∂G0(zn,y,ω)∂n(zn)=-jk4H1(1)(k|zn-y|)|zn-y|(zn-y)·n(z)---(19)]]>其中，G0(zn,y,ω)/n(zn)为关联zn点和声源点y之间的自由空间格林函数G0(zn,y,ω)沿zn点外法线方向的偏导数，n(zn)表示zn处边界的外法线方向。当所述边界元模型S为三维空间时，通过下列公式(20)关联zn点和声源点y之间的自由空间格林函数G0(zn,y,ω)沿zn点外法线方向的偏导数G0(zn,y,ω)/n(zn)：∂G0(zn,y,ω)∂n(zn)=14πejk|z-y|(jk|zn-y|-1)|zn-y|3(zn-y)·n(zn)---(20)]]>在步骤106中，根据上述步骤104公式(2)确定的G(x,zn,ω)、步骤105确定所述G0(zn,y,ω)和所述G0(zn,y,ω)/n(zn)，可以通过公式(3)确定散射格林函数GS(x,y,ω)。具体地，对于二维问题，可以将公式(2)确定的G(x,zn,ω)，公式(17)确定的G0(zn,y,ω)和公式(19)确定的G0(zn,y,ω)/n(zn)代入公式(3)内，确定散射格林函数GS(x,y,ω)；对于三维问题，可以将公式(2)确定的G(x,zn,ω)，公式(18)确定的G0(zn,y,ω)和公式(20)确定的G0(zn,y,ω)/n(zn)代入公式(3)内，确定散射格林函数GS(x,y,ω)。具体地，公式(3)如下所示：GS(x,y,ω)=∫SG(x,zn,ω)∂G0(zn,y,ω)∂n(zn)+β(zn)G0(zn,y,ω)ds(zn)---(3)]]>其中，GS(x,y,ω)为关联观察点x和声源点y之间的散射格林函数。在步骤107中，根据上述G0(x,y,ω)和GS(x,y,ω)，通过公式(4)确定精确格林函数G(x,y,ω)。具体的，对于二维问题，可以将公式(15)确定的G0(x,y,ω)和公式(3)确定的GS(x,y,ω)，代入公式(4)中，确定精确格林函数G(x,y,ω)；对于三维问题，可以将公式(16)确定的G0(x,y,ω)和公式(3)确定的GS(x,y,ω)，代入公式(4)中，确定精确格林函数G(x,y,ω)。具体地，公式(4)如下所示：G(x,y,ω)＝G0(x,y,ω)+GS(x,y,ω)(4)其中，G(x,y,ω)为为关联观察点x和声源点y之间的精确格林函数。综上所述，本发明实施例提供的一种确定精确格林函数的方法，采用了Robin边界条件，不仅能分析声学非紧致硬界的散射，还能考虑非紧致阻抗边界对声传播的影响，是一种应用范围更广的精确格林函数计算方法。为了清楚介绍本发明实施例提供的一种确定精确格林函数的方法，下列实施例一中，以一种二维非紧致圆柱外空间为例，介绍确定精确格林函数的算法：实施例一如图2所示，在空间中有一直径D＝100mm的二维圆柱，在圆柱附近y处有一单位强度的单极子点声源，声源的圆频率为ω，观察点x与圆柱中心的距离为r＝12800mm。对阻抗圆柱，其表面阻抗参数可以通过下列公式(21)确定：Z(ω)＝0.2ρ0c0+j(-13.48ρ0c0/ω+0.0739ρ0c0ω)(21)其中，Z(ω)为声学边界的阻抗，ρ0＝1.215kg/m3，c0＝340m/s。具体地，本发明实施例中，针对二维圆柱，确定精确格林函数的算法主要包括下列步骤：步骤201，划分网格，建立边界元模型，其中，边界元模型需要将结构表面进行离散，在本发明实施例一中，二维圆柱结构简单，采用自编程序完成网格划分，将180个离散结点(m＝1,2,…,180代表编号为m的离散点)均匀布置在圆柱表面，从而将二维圆柱表面离散为180个线网格单元。步骤202，提取各边界网格单元的中心点坐标zm、单位外法线方向n(zm)以及面积Δs(zm)，m＝1,2,…,180代表编号为m的网格单元。对于本发明实施例一，坐标原点位于圆柱中心，提取的过程包括：a.提取离散结点的坐标：从已建立好的边界元模型中，提取出各离散结点矢量的坐标，m＝1,2,…,180代表编号为m的离散点，和分别为第m个离散点在x和y轴方向的坐标值。b.确定各线网格单元的中心点坐标：第m个线网格单元由离散点和决定，其中心点矢量的坐标通过下列公式(22)确定：x1m=(x0m+x0m+1)/2,y1m=(y0m+y0m+1)/2---(22)]]>c.确定各线网格单元的单位外法线方向矢量：第m个线网格单元单位外法线方向矢量的坐标通过下列公式(23)确定：xnm=x1m/(x1m)2+(y1m)2,ynm=y1n/(x1m)2+(y1m)2---(23)]]>d.确定各线网格单元的面积：对于第m个线网格单元，其面积等于该线网格单元的长度，长度可以通过下列公式(24)确定：Δs(zm)=(x0m+1-x0m)2+(y0m+1-y0m)2---(24)]]>步骤203，根据边界的声学特性，通过公式(1)，确定参数β：β(zm)＝αjρ0ω/Z(zm,ω)β(zn)＝αjρ0ω/Z(zn,ω)(1)其中，ω表示圆频率，Z(ω)为声学边界的阻抗，j为虚数单位，ρ0为声传播介质的密度。需要说明的是，对于声学硬边界，β＝0；对声学阻抗边界时，β＝jρ0ω/Z(ω)。步骤204，对观察点x和边界元模型上的点z，采用边界元方法求解下述积分方程用于计算关联x和z的格林函数G(x,z,ω)：12G(x,zn,ω)=G0(x,zn,ω)+∫S′G(x,zm,ω)∂G0(zm,zn,ω)∂n(zm)+β(zm)G0(zm,zn,ω)ds(zm)---(2)]]>公式(2)中，x为观察点的位置矢量，z和是边界上网格单元的中心点，ω为圆频率，G0是频域自由空间格林函数，Φ(z)是边界上z点处的立体角函数，n(z)表示z点处边界的外法线方向，表示点所在的网格单元面积，S'是除去点部分的边界面，π是圆周率；β是声学边界条件参数。其中，公式(2)的具体计算过程如下：a.将公式(2)离散成一系列代数线性方程组，具体如公式(13)所示：|E2+H|G(x,z,ω)=G0(x,z,ω)---(13)]]>其中，E是单位对角矩阵，可以用下列公式(25)确定：E=10...01180×180---(25)]]>进一步地，公式(13)的H矩阵的子项可以通过下列公式(12)确定：Hmn=∫S′∂G0(zm,zn,ω)∂n(zm)+βG0(zm,zn,ω)ds(zm),m≠n0,m=n---(12)]]>公式(12)中，下标m,n＝1,2,…,180为网格单元编号，zm和zn分别是编号为m和n的网格单元的中心点，S'是除去zm点部分的边界面，Δs(zm)为编号为m的网格单元的面积。矩阵G(x,z,ω)和G0(x,z,ω)分别由下列公式(14)和公式(11)确定：G(x,z,ω)=G(x,z1,ω)G(x,z2,ω)...G(x,z180,ω)---(14)]]>G0(x,z,ω)=G0(x,z1,ω)G0(x,z2,ω)...G0(x,z180,ω)---(11)]]>b.对于二维空间的声传播问题，可以采用公式公式(5)和公式(7)确G0(zm,zn,ω)和G0(zm,zn,ω)/n(zm)，同时，采用公式(9)确定G0(x,zn,ω)：G0(zm,zn,ω)=j4H0(1)(k|zm-zn|)---(5)]]>∂G0(zm,zn,ω)∂n(zm)=-jk4H1(1)(k|zm-zn|)|zm-zn|(zm-zn)·n(zm)---(7)]]>G0(x,zn,ω)=j4H0(1)(k|x-zn|)---(9)]]>其中，为0阶第一类Hankel函数，为1阶第一类Hankel函数，j为虚数单位，k＝ω/c0为声学波数，ω为圆频率，c0为声波在介质中传播的速度。c.将根据公式(5)和公式(7)确定的G0(zm,zn,ω)和G0(zm,zn,ω)/n(zm)代入公式(12)后，采用Gauss-Legendre方法积分求解H矩阵。d.将公式(9)、矩阵E和矩阵H，采用Gauss方法求解代数线性方程组(13)，得到格林函数G(x,z,ω)。步骤205，对观察点x和源点y，采用边界元方法求解下述积分方程用于计算关联x和y的精确格林函数G(x,y,ω)：G(x,y,ω)＝G0(x,y,ω)+GS(x,y,ω)(4)其中，x和y分别为观察点和源点的位置矢量，z是边界上网格单元的中心点，ω为圆频率，G0是频域自由空间格林函数，n(z)表示z点处边界的外法线方向，s(z)表示z点所在的网格单元面积，S是非紧致边界；β为声学边界条件参数。精确格林函数G(x,y,ω)的具体计算过程如下：a.确定G0(x,y,ω)、G0(z,y,ω)和G0(z,y,ω)/n(z)的解。对于二维空间的声传播问题，可以采用公式(17)和公式(19)确定G0(zn,y,ω)和G0(zn,y,ω)/n(zn)，同时，采用公式(15)确定G0(x,y,ω)：G0(x,y,ω)=j4H0(1)(k|x-y|)---(15)]]>G0(zn,y,ω)=j4H0(1)(k|zn-y|)---(17)]]>∂G0(zn,y,ω)∂n(zn)=-jk4H1(1)(k|zn-y|)|zn-y|(zn-y)·n(z)---(19)]]>其中，为0阶第一类Hankel函数，为1阶第一类Hankel函数，j为虚数单位，k＝ω/c0为声学波数，ω为圆频率，c0为声波在介质中传播的速度。b.将上述公式(21)，公式(23)结合下列公式(3)，可以确定散射格林函数GS(x,y,ω)：GS(x,y,ω)=∫SG(x,zn,ω)∂G0(zn,y,ω)∂n(zn)+β(zn)G0(zn,y,ω)ds(zn)---(3)]]>c.将公式(19)确定的自由空间格林函数G0(x,y,ω)与公式(3)确定的散射格林函数GS(x,y,ω)通过公式(4)，可以确定精确格林函数G(x,y,ω)：G(x,y,ω)＝G0(x,y,ω)+GS(x,y,ω)(4)在实际应用中，图3为本发明实施例一提供的圆柱声学硬边界条件下精确格林函数数值解、解析解与自由空间格林函数的对比图；图4为本发明实施例一提供的圆柱声学阻抗边界条件下精确格林函数数值解、解析解与声学硬边界条件下精确格林函数数值解的对比图。当声学波数k＝20时，声波波长与圆柱直径的比约为3，圆柱不满足声学紧致条件。此时，对于声学硬边界圆柱，采用本发明实施例所提供的一种确定精确格林函数数值解与精确格林函数解析解以及自由空间格林函数的对比如图3所示，其中，图3中曲线5为自由空间格林函数结果，曲线6为精确格林函数解析解结果，曲线7为精确格林函数数值解结果。采用本发明实施例提供的确定精确格林函数与自由空间格林函数存在明显的差异，是因为本发明实施例考虑了非紧致圆柱边界的声散射。对于声学阻抗边界圆柱，采用本发明方实施例得到的精确格林函数数值解与精确格林函数解析解以及声学硬边界圆柱精确格林函数的对比如图4所示，图4中曲线8为声学硬边界圆柱精确格林函数数值解结果，曲线9为声学阻抗边界圆柱精确格林函数解析解结果，曲线10为声学阻抗边界圆柱精确格林函数数值解结果。采用本发明实施例得到的声学阻抗边界圆柱精确格林函数与声学硬边界圆柱精确格林函数存在明显的差异，是因为本发明实施例考虑了圆柱阻抗边界的声吸收。从图3和图4可以看出，无论圆柱表面是声学硬边界条件还是声学阻抗边界条件，采用本发明实施例提供精确格林函数数值解均与解析解完全吻合，从而可以证明了本发明实施例所提供的确定精确格林函数的方法的正确性。基于同一发明构思，本发明实施例提供了一种确定精确格林函数的装置，由于该装置解决技术问题的原理与一种确定精确格林函数的方法相似，因此该装置的实施可以参见方法的实施，重复之处不再赘述。图5为本发明实施例提供的一种确定精确格林函数的装置结构示意图，如图5所示，本发明实施例所提供的一种确定精确格林函数的装置主要包括：建立单元501，第一确定单元502，第二确定单元503，第三确定单元504，第四确定单元505，第五确定单元506和第六确定单元507。建立单元501，用于建立边界元模型，将所述边界S离散成M个网格单元，依次确定所述M个网格单元的中心点，单位外法线方向以及面积；第一确定单元502，用于在所述M个网格单元中确定任意两个中心点为zm和zn的网格单元，根据所述zm和所述zn点所在散射边界的声学特性以及公式(1)，确定所述zm和所述zn点处的声学参数β(zm)和β(zn)，确定关联所述zm点和所述zn点的自由空间格林函数G0(zm,zn,ω)以及所述自由空间格林函数G0(zm,zn,ω)沿所述zm点外法线方向的偏导数G0(zm,zn,ω)/n(zm)；第二确定单元503，用于在所述边界元模型S外确定观察点x，根据所述zn点和所述观察点x，确定关联所述观察点x和所述zn点之间的自由空间格林函数G0(x,zn,ω)；第三确定单元504，用于根据所述G0(zm,zn,ω)、所述G0(zm,zn,ω)/n(zm)和所述G0(x,z,ω)，通过公式(2)确定关联所述观察点x和所述zn点之间的格林函数G(x,zn,ω)；第四确定单元505，用于在所述边界S外确定声源点y，根据所述观察点x和所述声源点y，确定关联所述观察点x和所述声源点y之间的自由空间格林函数G0(x,y,ω)；根据所述zn点和所述声源点y，确定关联所述zn点和所述声源点y之间的自由空间格林函数G0(zn,y,ω)及所述格林函数沿zn点外法线方向的偏导数G0(zn,y,ω)/n(zn)；第五确定单元506，用于根据所述G(x,zn,ω)、所述G0(zn,y,ω)和所述G0(zn,y,ω)/n(zn)通过公式(3)确定散射格林函数GS(x,y,ω)；第六确定单元507，用于根根据所述G0(x,y,ω)和所述GS(x,y,ω)，通过公式(4)确定精确格林函数G(x,y,ω)；其中，公式(1)如下所示：β(zm)＝αjρ0ω/Z(zm,ω)，β(zn)＝αjρ0ω/Z(zn,ω)公式(2)如下所示：12G(x,zn,ω)=G0(x,zn,ω)+∫S′G(x,zm,ω)∂G0(zm,zn,ω)∂n(zm)+β(zm)G0(zm,zn,ω)ds(zm)]]>公式(3)如下所示：GS(x,y,ω)=∫SG(x,zn,ω)∂G0(zn,y,ω)∂n(zn)+β(zn)G0(zn,y,ω)ds(zn)]]>公式(4)如下所示：G(x,y,ω)＝G0(x,y,ω)+GS(x,y,ω)其中，ω表示圆频率，j为虚数单位，ρ0为声传播介质的密度，G0是频域自由空间格林函数；x和y分别为观察点和源点的位置矢量，zm和zn是边界上网格单元的中心点；n(zm)和n(zn)分别表示zm和zn处边界的外法线方向，s(zm)和s(zn)分别表示zm和zn点所在的网格单元面积，S是散射边界，S'是除去zm点所在网格单元的散射边界；Z(zm)和Z(zn)分别表示zm和zn点所在的网格单元声学阻抗，对于声学硬边界α＝0，对于声学阻抗边界α＝1。优选地，所述第一确定单元502具体用于：当所述边界元模型为二维空间时，通过下列公式确定G0(zm,zn,ω)：G0(zm,zn,ω)=j4H0(1)(k|zm-zn|)]]>当所述边界元模型为三维空间时，通过下列公式确定G0(zm,zn,ω)：G0(zm,zn,ω)=14π|zm-zn|ejk|zm-zn|]]>当所述边界元模型为二维空间时，通过下列公式确定G0(zm,zn,ω)/n(zm)：∂G0(zm,zn,ω)∂n(zm)=-jk4H1(1)(k|zm-zn|)|zm-zn|(zm-zn)·n(zm)]]>当所述边界元模型为三维空间时，通过下列公式确定G0(zm,zn,ω)/n(zm)：∂G0(zm,zn,ω)∂n(zm)=14πejk|zm-zn|(jk|zm-zn|-1)|zm-zn|3(zm-zn)·n(zm)]]>所述第二确定单元503具有用于：当所述边界元模型为二维空间时，通过下列公式确定G0(x,zn,ω)：G0(x,zn,ω)=j4H0(1)(k|x-zn|)]]>当所述边界元模型为三维空间时，通过下列公式确定G0(x,zn,ω)：G0(x,zn,ω)=14π|x-zn|ejk|x-zn|]]>其中，下标m,n＝1,2,…,M为网格单元编号，zm和zn分别是编号为m和n的网格单元的中心点，x和y分别为观察点和源点的位置矢量，S'是除去zm点部分的边界面，n(zm)表示zm点处边界的外法线方向，为0阶第一类Hankel函数，j为虚数单位，k＝ω/c0为声学波数，ω为圆频率，c0为声波在介质中传播的速度。优选地，所述第三确定单元504还用于：根据所述所有zn点和所述G0(x,zn,ω)，通过下列公式确定矩阵G0(x,z,ω)G0(x,z,ω)=G0(x,z1,ω)G0(x,z2,ω)...G0(x,zM,ω)]]>根据所述G0(zm,zn,ω)和所述G0(zm,zn,ω)/n(zm)，通过下列公式确定H矩阵的子项：Hmn=∫S′∂G0(zm,zn,ω)∂n(zm)+βG0(zm,zn,ω)ds(zm),m≠n0,m=n]]>根据所述H矩阵的子项和所述G0(x,z,ω)矩阵，通过下列公式确定代数离散方程组：|E2+H|G(x,z,ω)=G0(x,z,ω)]]>其中，E是单位对角矩阵，Hmn为矩阵H下标为m,n的子项，矩阵G(x,z,ω)为对应于所述所有zn点需要求解的未知矩阵G(x,z,ω)=G(x,z1,ω)G(x,z2,ω)...G(x,zM,ω)]]>优选地，所述第四确定单元505具体用于：当所述边界元模型为二维空间时，通过下列公式确定G0(zn,y,ω)：G0(zn,y,ω)=j4H0(1)(k|zn-y|)]]>当所述边界元模型为三维空间时，通过下列公式确定G0(zn,y,ω)：G0(zn,y,ω)=14π|zn-y|ejk|z-y|]]>确定关联所述zn点和所述声源点y之间的自由空间格林函数G0(zn,y,ω)沿zn点外法线方向的偏导数G0(zn,y,ω)/n(zn)，包括：当所述边界元模型为二维空间时，通过下列公式确定G0(zn,y,ω)/n(zn)：∂G0(zn,y,ω)∂n(zn)=-jk4H1(1)(k|zn-y|)|zn-y|(zn-y)·n(z)]]>当所述边界元模型为三维空间时，通过下列公式确定G0(zn,y,ω)/n(zn)：∂G0(zn,y,ω)∂n(zn)=14πejk|z-y|(jk|zn-y|-1)|zn-y|3(zn-y)·n(zn)]]>确定关联所述观察点x和所述声源点y之间的自由空间格林函数G0(x,y,ω)，包括：当所述边界元模型为二维空间时，通过下列公式确定G0(x,y,ω)：G0(x,y,ω)=j4H0(1)(k|x-y|)]]>当所述边界元模型为三维空间时，通过下列公式确定G0(x,y,ω)：G0(x,y,ω)=14π|x-y|ejk|x-y|]]>其中，x和y分别为观察点和源点的位置矢量，n(zn)表示zn点处边界的外法线方向，j为虚数单位，为0阶第一类Hankel函数，ω为圆频率，c0为声波在介质中传播的速度，k＝ω/c0为声学波数。应当理解，以上一种确定精确格林函数的装置包括的单元仅为根据该设备装置实现的功能进行的逻辑划分，实际应用中，可以进行上述单元的叠加或拆分。并且该实施例提供的一种确定精确格林函数的装置所实现的功能与上述实施例提供的一种确定精确格林函数的方法一一对应，对于该装置所实现的更为详细的处理流程，在上述方法实施例一中已做详细描述，此处不再详细描述。本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。当前第1页1 2 3