CN113031263A - 基于小波边界元模型的二维正方晶格光子晶体带隙设计方法 - Google Patents
基于小波边界元模型的二维正方晶格光子晶体带隙设计方法 Download PDFInfo
- Publication number
- CN113031263A CN113031263A CN202110335850.7A CN202110335850A CN113031263A CN 113031263 A CN113031263 A CN 113031263A CN 202110335850 A CN202110335850 A CN 202110335850A CN 113031263 A CN113031263 A CN 113031263A
- Authority
- CN
- China
- Prior art keywords
- photonic crystal
- band gap
- boundary
- matrix
- lattice photonic
- 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
Links
- 239000004038 photonic crystal Substances 0.000 title claims abstract description 68
- 238000000034 method Methods 0.000 title claims abstract description 56
- 238000013461 design Methods 0.000 title claims abstract description 17
- 239000011159 matrix material Substances 0.000 claims abstract description 55
- 238000004364 calculation method Methods 0.000 claims abstract description 32
- RKTYLMNFRDHKIL-UHFFFAOYSA-N copper;5,10,15,20-tetraphenylporphyrin-22,24-diide Chemical compound [Cu+2].C1=CC(C(=C2C=CC([N-]2)=C(C=2C=CC=CC=2)C=2C=CC(N=2)=C(C=2C=CC=CC=2)C2=CC=C3[N-]2)C=2C=CC=CC=2)=NC1=C3C1=CC=CC=C1 RKTYLMNFRDHKIL-UHFFFAOYSA-N 0.000 claims abstract description 24
- 238000005516 engineering process Methods 0.000 claims abstract description 5
- 239000013598 vector Substances 0.000 claims description 37
- 230000006870 function Effects 0.000 claims description 33
- 239000003637 basic solution Substances 0.000 claims description 8
- 230000003287 optical effect Effects 0.000 claims description 7
- 239000000758 substrate Substances 0.000 claims description 7
- 229910000831 Steel Inorganic materials 0.000 claims description 6
- 239000010959 steel Substances 0.000 claims description 6
- 230000005684 electric field Effects 0.000 claims description 5
- 238000006073 displacement reaction Methods 0.000 claims description 4
- 239000011248 coating agent Substances 0.000 claims description 3
- 238000000576 coating method Methods 0.000 claims description 3
- 238000010586 diagram Methods 0.000 claims description 3
- 239000002184 metal Substances 0.000 claims description 3
- 125000004430 oxygen atom Chemical group O* 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 239000000463 material Substances 0.000 claims description 2
- 230000008901 benefit Effects 0.000 abstract description 5
- 230000009467 reduction Effects 0.000 abstract description 2
- 239000000243 solution Substances 0.000 description 15
- 239000006185 dispersion Substances 0.000 description 5
- 230000010354 integration Effects 0.000 description 5
- 230000008569 process Effects 0.000 description 4
- 239000007787 solid Substances 0.000 description 4
- 239000002585 base Substances 0.000 description 3
- 238000012546 transfer Methods 0.000 description 3
- 239000013078 crystal Substances 0.000 description 2
- 230000000737 periodic effect Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 235000006629 Prosopis spicigera Nutrition 0.000 description 1
- 240000000037 Prosopis spicigera Species 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000000295 complement effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011438 discrete method Methods 0.000 description 1
- 238000011549 displacement method Methods 0.000 description 1
- 239000008204 material by function Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012916 structural analysis Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G02—OPTICS
- G02B—OPTICAL ELEMENTS, SYSTEMS OR APPARATUS
- G02B27/00—Optical systems or apparatus not provided for by any of the groups G02B1/00 - G02B26/00, G02B30/00
- G02B27/0012—Optical design, e.g. procedures, algorithms, optimisation routines
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Computing Systems (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Algebra (AREA)
- Optics & Photonics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Optical Integrated Circuits (AREA)
- Optical Modulation, Optical Deflection, Nonlinear Optics, Optical Demodulation, Optical Logic Elements (AREA)
Abstract
本发明公开了一种基于小波边界元模型的二维正方晶格光子晶体带隙设计方法，包括：S1:将区间B样条小波与边界元法相结合，用BSWI尺度函数取代传统边界元的多项式插值，结合单胞技术，获得了二维正方晶格光子晶体基体与散射体统一的离散化边界积分方程形式，进而获得代数方程组；S2:根据步骤S1得到的代数方程组，在频率域内，结合Bloch定理以及基体与散射体间的连续性条件，进一步建立二维正方晶格光子晶体带隙特性计算模型；S3:通过调整二维正方晶格光子晶体基体或散射体尺寸，获得所需要的带隙特性，最终完成二维正方晶格光子晶体带隙设计。本发明的方法既兼有边界元方法降维的特性，又具备小波多尺度逼近的优势，适合于二维正方晶格光子晶体带隙设计。
Description
技术领域
本发明属于光学功能材料结构设计领域，具体是指基于小波边界元模型的二维正方晶格光子晶体带隙设计方法。
背景技术
光子晶体是一种由基体和散射体组成的人造周期性结构。它最突出的特性就是光子带隙，即处于光子晶体带隙频率范围内的光波在晶体中被禁止传播。带隙又可分为完全带隙和不完全带隙(方向带隙)。若在一个频率范围内，光波在任何方向上都不能通过，则称为完全带隙，若在特定频率范围内只允许某个方向的光波传播，则称为不完全带隙。因此，许多学者致力于光子晶体带隙设计方面的研究旨在设计出一种具有良好带隙特性的光子晶体结构。
准确地计算出光子晶体带隙特性是带隙设计的基础。目前已经开发出许多光子晶体带隙特性的计算方法，比如平面波展开法、时域有限差分法、DtN-map法、有限元法、传递矩阵法、散射矩阵法和边界元法等。这些方法可以划分两大类，一类是Bloch波矢作为给定的参数，角频率ω是特征值方程的待求特征值，例如平面波展开法、时域有限差分法和有限元法。其中，平面波展开法可以简单地认为是将非均匀波动方程中的物理量展开为无穷级数形式，但是其收敛较慢，进一步限制了求解精度。时域有限差分法需要对时间和空间域进行足够的离散化，因此，非常耗时。此外，这种方法无法应用于形状复杂的晶体。有限元法是通过将单胞离散化来实现光子晶体能带结构的求解，并有良好的收敛性、兼容性和准确性。但是，它涉及大规模矩阵的求解，需要花费很多的时间成本。另外，对于上述情况，如果涉及非色散介质，则最终归结于线性特征值方程的求解问题，如果是涉及色散介质，则为较复杂非线性特征值方程的求解问题。另一种情况是假设角频率ω为已给定参数，特征值方程与Bloch波矢有关。对于此类处理方法，不管涉及色散介质还是非色散介质，最终都可以转化为矩阵规模小的线性特征值问题。比如说，DtN-map法，传递矩阵法，散射矩阵法。但是，DtN-map法仅局限于求解圆形散射体的情况。传递矩阵方法可以适用于各种形状的散射体，但过程中会存在求解不稳定的问题，需要用特殊的方法处理。有许多方法可以求解散射矩阵，但是涉及很复杂的求解技术，限制了其广泛应用。边界元法具有求解规模小，凭借基本解可以自动满足辐射条件的优势，因此，最近几年也被用于光子晶体能带结构的计算。
总之，虽然上述有些算法已经得到了大量应用，但共性缺点在于精度低、收敛性差等问题，这制约了二维正方晶格光子晶体带隙设计应用于工程实践。
发明内容
本发明的目的是为了克服现有技术存在的缺点和不足，而提供一种基于小波边界元模型的二维正方晶格光子晶体带隙设计方法。
为实现上述目的，本发明的技术方案是包括以下步骤：
S1:将区间B样条小波与边界元法相结合，用BSWI尺度函数取代传统边界元的多项式插值，结合单胞技术，获得了二维正方晶格光子晶体基体与散射体统一的离散化边界积分方程形式，进而获得代数方程组；
S2:根据步骤S1得到的代数方程组，在频率域内，结合Bloch定理以及基体与散射体间的连续性条件，进一步建立二维正方晶格光子晶体带隙特性计算模型；
S3:通过调整二维正方晶格光子晶体基体或散射体尺寸，获得所需要的带隙特性，最终完成二维正方晶格光子晶体带隙设计。
进一步设置是步骤S1包括以下步骤：
1)运用一维BSWI尺度函数作为插值函数，得离散化的边界积分方程
其中，P和Q分别代表源点和场点，Ne代表单元的个数，c(P)＝β/2π表示与源点P处边界形状有关的系数，β为P处切线张角，u*(P,Q)为本光学问题基本解，q*(P,Q)为基本解沿外法线方向的方向导数，li为单元的长度，
和
分别表示第i个单元的电场/磁场值和其法向导数值所组成的列向量，
指BSWI尺度函数所组成的行向量，Te为本光学问题所对应的转换矩阵，u*(P,Q)，q*(P,Q)，
及Te的表达式分别为：
其中，ε指材料的介电常数，k0＝ω/c0是自由空间波数，c0表示真空中的光速，ω是角频率，r＝|xP-xQ|表示源点与场点的距离，
表示第一类0阶汉克尔函数；
其中，
表示第一类1阶汉克尔函数，xi(Q)与xi(P)分别表示源点与场点的坐标点，ni(Q)表示场点处的方向余弦；
其中，m与j分别表示BSWI尺度函数的阶数和尺度，ξ∈[0,1]为局部坐标；
其中，ξi为第i个节点的局部坐标值，N表示小波单元节点的个数；
将每一个节点都设为源点，经过积分运算、矩阵组装可进一步得到代数方程组：
[H]NP×NP{U}NP×1+[G]NP×NP{Q}NP×1＝0。
其中，H和G为系统矩阵，U和Q分别表示所有节点位移和位移法向导数所组成的列向量，NP为节点的总数；
进一步设置是步骤S2包括以下步骤：
1)通过将计算基体与散射体的系统矩阵的子矩阵进行整合，得
其中，对于TM模式，η＝1，对于TE模式，η＝ε0/ε1，系数η是求解带隙特性的关键，kx,ky为第一布里渊区Bloch波矢，a为晶格常数，
和
表示场点在基体的边界Γ1上，将所有节点都视为源点时积分分别涉及u*(P,Q)和q*(P,Q)所得到的矩阵，其余类似，另外有：
分离出含待求Bloch波矢的项，可获得二维正方晶格光子晶体的BSWI边界元带隙计算模型
AX＝ξBX
其中，
或
2)由于第一布里渊区具有对称性，带隙计算通常沿简约布里渊区的边界进行，在正方晶格简约布里渊区Γ-X-M-Γ的每个边界上涉及到的矩阵A、B为
(1)Γ-X上，
(2)X-M上，
(3)M-Γ上，
根据上述矩阵A、B，可求得正方晶格简约布里渊区Γ-X-M-Γ的每个边界上未知的Bloch波矢，计算完毕之后，用简约波矢M、Γ、X作为横坐标，所求的Bloch波矢作为所在边界横坐标x上的值、归一化频率为ωa/(2πc)纵坐标y，根据归一化频率与所求Bloch波矢的对应关系，就可得到二维正方晶格光子晶体能带结构图，从而获得正方晶格光子晶体带隙特性。
本发明的有益效果：本发明将BSWI与边界元法相结合，用BSWI尺度函数取代传统边界元的多项式插值，结合单胞技术，建立了求解二维正方晶格光子晶体基体与散射体的代数方程组。之后，基于Bloch定理、基体与散射体连续性条件，构造了二维正方晶格光子晶体带隙计算模型，进而获得光子晶体带隙特性。具有下列区别于传统边界元求解方法的显著优势：
1)本发明结合传统边界元的降维特性与B样条函数的优良逼近特性进行结构分析，在光子晶体带隙计算过程中，用精确的BSWI尺度函数取代传统的多项式插值来构造形状函数，进而构成小波单元，最终能够以较少的单元与积分点数获得较高的计算精度。这使得小波边界元模型较传统边界元法具有更好的计算效率和收敛性；
2)本发明独创性地在边界元法中用多节点单元来进行二维正方晶格光子晶体带隙特性计算，通过给定的角频率ω，最终可以获得容易求解的广义线性特征值方程组，即二维正方晶格光子晶体带隙计算模型；
3)通过构建二维正方晶格光子晶体带隙设计的小波边界元模型，不断调整光子晶体基体或散射体尺寸，可高精度、快速收敛地获得所需要的带隙特性，最终完成正方晶格光子晶体带隙设计。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，根据这些附图获得其他的附图仍属于本发明的范畴。
图1是本发明的正方晶格第一布里渊区；
图2是本发明的正方晶格基体(a)和散射体(b)的域和边界；
图3是本发明的小波边界元模型计算的带隙特性(空心点)和传统边界元模型计算的带隙特性(实心点)；
图4是本发明的使用不同阶次尺度函数的小波边界元模型计算的带隙特性；
图5是本发明的使用不同积分点数的传统边界元模型计算的带隙特性；
图6本发明的实施流程图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步地详细描述。
如图6所示的本发明实施方案包括以下步骤：
S1:将区间B样条小波与边界元法相结合，用BSWI尺度函数取代传统边界元的多项式插值，结合单胞技术，获得了二维正方晶格光子晶体基体与散射体统一的离散化边界积分方程形式，进而获得代数方程组；
S2:根据步骤S1得到的代数方程组，在频率域内，结合Bloch定理以及基体与散射体间的连续性条件，进一步建立二维正方晶格光子晶体带隙特性计算模型；
S3:通过调整二维正方晶格光子晶体基体或散射体尺寸，获得所需要的带隙特性，最终完成二维正方晶格光子晶体带隙设计。
进一步地，本发明实施例中，本发明将BSWI与边界元法相结合，用BSWI尺度函数取代传统边界元的多项式插值。考虑到光子晶体的周期性，在一个单胞内建立了基体与散射体的边界积分方程组，此后结合Bloch理论、基体与散射体连续性条件，进而构造了二维正方晶格光子晶体带隙特性计算模型。
二维光子晶体问题的控制方程为：
其中，对于TM模式，u＝E，对于TE模式，u＝H，
指Laplace算子，Ω为求解域，k0＝ω/c0为自由空间波数，ω是角频率，ε是介电常数，c0表示真空中的光速。
基于基本解与格林公式，可以得到方程(1)的边界积分方程：
c(P)u(P)＝∫Γu*(P,Q)q(Q)-u(Q)q*(P,Q)dΓ(Q) (2)
其中，P和Q分别指边界上的源点和场点，c是与源点P处边界形状有关的系数，u*(P,Q)和q*(P,Q)分别是基本解及其法向导数：
其中，r＝|xP-xQ|表示源点与场点的距离，
和
分别是第一类0阶和1阶汉克尔函数。
由于光子晶体是周期型结构，只需要研究其一个单胞，但基体外边界所有边界变量必须满足Bloch理论，用统一的表达式表示为：
γ(x+l)＝eik·lγ(x) (5)
其中，l＝m1ax+m2ay，ax和ay表示晶格基向量，m1和m2是大小介于0和晶格常数a的系数，γ是所有边界变量的统称，k＝(kx,ky)是Bloch波矢。
为了区分与基体和散射体有关的量，本专利中所有带有下标“0”的量都与基体有关，带有下标“1”的量与散射体有关。
由于边界变量在基体和散射体之间的交界面处是连续的，因此必须满足以下关系式：
式(6)和(7)分别表示TM和TE模式的基体与散射体交界面连接条件，负号表示基体边界与散射体边界外法线方向相反。
采用一维BSWI尺度函数插值时，单元的电场/磁场值及其法向导数可以表示为：
其中，
表示m阶j尺度的尺度函数组成的行向量，Te为对应问题的转换矩阵，见公式(9)，ue，qe为物理空间自由度列向量，见公式(10)：
其中，大写字母N表示阶数为m，尺度为j的BSWI尺度函数
的维数，同时也表示一个BSWI单元的节点总数。
公式(2)被离散为Ne个单元可以表示为：
其中，li表示第i个单元的总长度，
和
分别表示第i个单元的电场/磁场值和其法向导数值所组成的列向量。
公式(11)中的项都可以直接通过积分得到。另外，转换矩阵Te可以根据所给的尺度函数表达式计算，将这些都考虑在内之后，得：
其中，h(P,i,1)与g(P,i,1)表示第一个尺度函数作为被考察函数时与基本解所得到的积分值，以此类推。将求和公式展开并将涉及到同一节点的系数进行叠加，可得：
其中，NP代表总的节点个数，U和Q分别表示所有节点处的电场/磁场值和其法向导数值所组成的列向量。如果将所有节点都视为源点，则可以形成代数方程组为：
[H]NP×NP{U}NP×1+[G]NP×NP{Q}NP×1＝0 (14)
其中：
通常情况下会采用刚体位移法确定H矩阵主对角线值，这样也可以间接求出系数c的值，即：
图2给出了一个正方晶格晶胞中基体和散射体的边界、求解域及求解路径。其中，Γ1～Γ4为基体的外边界，边界上相应的节点电场/磁场值和其法向导数分别用
和
表示。Γ0和Γ0'分别表示基体的内边界和散射体的外边界，边界上节点电场/磁场值和其法向导数分别为
和
另外，计算域分别由D0和D1表示。
为了构造二维正方晶格光子晶体带隙计算模型，假设边界Γ1～Γ4分别离散化为N1个BSWI单元，Γ0和Γ0'离散化为N2个BSWI单元。此后，公式(14)可以写成以下形式：
其中，上标α和γ分别表示与源点在基体边界Γα和散射体边界Γγ有关的量，β表示与场点在边界Γβ有关的量。根据基体和散射体交界面连续需要满足的关系式(6)、(7)，公式(17)和(18)可以进一步组装为：
其中，对于TM模式，η＝1，对于TE模式，η＝ε0/ε1，其它在公式(19)中出现的矩阵和向量可以表示为：
基于Bloch理论(5)，基体外边界上的边界变量具有以下关系：
对于边界变量法向导数，由于外法向量方向是相反的，因此，满足：
将公式(23)和(24)带入(19)，并经子矩阵之间的组合，得：
其中：
分离出所有含Bloch波矢的项，得：
由于第一布里渊区具有对称性，带隙特性的计算通常沿正方晶格简约布里渊区的边界进行，如图1中Γ-X-M-Γ所组成的三角形。假设
由于在Γ-X上ky＝0，因此，ξy＝1。同样的，在边界上X-M上，ξx＝-1，在边界M-Γ上，ξx＝ξy。对于以上三种情况，其广义的线性特征值方程组统一表示为：
AX＝ξBX (28)
式(28)为计算二维正方晶格光子晶体带隙计算模型，其中，ξ＝ξx或ξ＝ξy。假如使用NP0和NP1表示基体和散射体的离散节点数，在简约布里渊区每个边界上的A和B为：
(1)Γ-X上，
(2)X-M上，
(3)M-Γ上，
以ξ＝ξx为例，根据欧拉公式，得：
因此，对于每个给定角频率ω，都要保证计算出的特征值|ξ|＝1。考虑到误差的存在，则有：
||ξ|-1|≤δ (30)
其中，δ是充分小的正数。选择一个合理的δ值至关重要，当δ太小时会出现丢根的情况，而δ太大则会导致伪根的存在。
根据公式(30)可求得正方晶格简约布里渊区Γ-X-M-Γ的每个边界上未知的Bloch波矢，可用简约波矢M、Γ、X作为横坐标，所求的Bloch波矢作为所在边界横坐标x坐标值、归一化频率ωa/(2πc)为纵坐标y，根据归一化频率与所求Bloch波矢的对应关系，就可得到二维正方晶格光子晶体能带结构图，从而获得正方晶格光子晶体带隙特性。
实施例1：本实施例主要验证二维正方晶格光子晶体带隙计算的小波边界元数值求解模型的计算效率与可靠性。在真空中(ε0＝1)嵌入圆形的电介质散射体(ε1＝8.9)。其中，填充比f＝0.4489，δ＝5×10-2，最小归一化频率τ0＝5×10-3，最大归一化频率τmax＝0.5，计算步长Δτ＝5×10-3。此外，计算过程中使用2点积分公式，为了表达方便，在没有明确说明的情况下，所说的单元数为基体的离散单元数。
分别采用8个BSWI23单元和64个传统的线性单元离散基体边界计算该体系，基体内边界的单元数等于外边界的单元数。图3给出了两种离散方法在该体系TE模式下的求解结果，其中，c0＝3×108表示光在真空中的速度，实心点和空心点分别表示64个线性单元和8个BSWI23单元的求解结果。从图中可以发现，两种单元的计算结果很吻合，但在该频率范围内没有光子晶体带隙。此外，当采用BSWI多节点单元时，由于奇异积分的计算次数和相邻单元的叠加次数很少，这使得小波边界元模型的计算成本(39.8s)低于传统边界元的计算成本(55.0s)，相比之下，小波边界元模型的计算时间压缩了27.6％。另外，传统的边界元法已成功地应用于二维正方晶格光子晶体带隙特性计算，从而验证了小波边界元计算模型的可靠性和高效性。
实施例2：本实施例主要验证二维正方晶格光子晶体带隙计算的小波边界元数值求解模型的收敛性与灵活性。在真空中(ε0＝1)嵌入正方形的电介质散射体(ε1＝8.9)。其中，填充比f＝0.1257，δ＝5×10-2，最小归一化频率τ0＝5×10-3，最大归一化频率τmax＝0.8，计算步长Δτ＝5×10-3。此外，计算过程中使用2点积分公式。
基体和散射体的每条边由1个BSWI23单元离散，即基体单元数为8，散射体单元数为4。根据公式(30)，最终可以得到310个特征值，然后根据欧拉公式求得Bloch波矢，最终按照归一化频率和波矢的对应关系得到该体系在TM模式下的带隙特性，如图4空心点所示，阴影部分为该模式下的光子带隙。可以看出，在低频域中结果较好，但在高频域中的色散曲线不完整，这是因为求解精度不够导致在这些频率下计算出的特征值不满足公式(30)。为了补全这些色散曲线，将尺度函数的阶数从2提高到4，单元个数、积分点数和其他参数保持不变。最终，根据公式(30)可以得到322个特征值结果，高频域中的色散曲线也很完整，如图4实心点所示。因此，可以根据小波边界元模型多种小波基函数可选择的优势，自由地提升BSWI尺度函数的阶数(或尺度)来实现所需的精度要求。
另外，为了说明BSWI高阶单元的优越性，同时选择用传统三次单元求解该体系。其中，基体和散射体的每个边界由3个三次单元离散，即基体的总单元数为24，散射体总单元数为12，其他参数保持不变。按照同样的方法，获得了270个求解结果如图5所示(空心点)。可以看出，高频域的色散曲线不完整，这是因为2点积分公式无法满足求解的收敛性需求。将积分点从2个增加到6个，此后，323个求解结果如图5所示(实心点)，色散曲线也都表现完整。这证明了小波边界元良好的收敛性，与传统单元相比，仅采用2点积分公式和少量的单元就可以满足BSWI43单元的收敛性要求。
以上数值算例表明，与传统边界元模型相比，小波边界元模型在计算二维正方晶格光子晶体带隙方面具有计算效率高、灵活性好、收敛速度快的特点。最后，可通过构建的二维正方晶格光子晶体带隙设计的小波边界元模型，不断调整正方晶格光子晶体基体与散射体尺寸，可高性能(高效、灵活、收敛)地获得所需要的带隙特性，最终完成二维正方晶格光子晶体带隙设计，获得符合特定要求带隙特性的正方晶格光子晶体。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，所述的程序可以存储于一计算机可读取存储介质中，所述的存储介质，如ROM/RAM、磁盘、光盘等。
以上所揭露的仅为本发明较佳实施例而已，当然不能以此来限定本发明之权利范围，因此依本发明权利要求所作的等同变化，仍属本发明所涵盖的范围。
应当注意，本发明的实施例可以通过硬件、软件或者软件和硬件的结合来实现。硬件部分可以利用专用逻辑来实现；软件部分可以存储在存储器中，由适当的指令执行系统，例如微处理器或者专用设计硬件来执行。本领域的技术人员可以理解上述的设备和方法可以使用计算机可执行指令和/或包含在处理器控制代码中来实现，例如在可编程的存储器或者诸如光学或电子信号载体的数据载体上提供了这样的代码。
此外，尽管在附图中以特定顺序描述了本发明方法的操作，但是，这并非要求或者暗示必须按照该特定顺序来执行这些操作，或是必须执行全部所示的操作才能实现期望的结果。相反，流程图中描绘的步骤可以改变执行顺序。附加地或备选地，可以省略某些步骤，将多个步骤组合为一个步骤执行，和/或将一个步骤分解为多个步骤执行。还应当注意，根据本发明的两个或更多装置的特征和功能可以在一个装置中具体化。反之，上文描述的一个装置的特征和功能可以进一步划分为由多个装置来具体化。
虽然已经参考若干具体实施例描述了本发明，但是应当理解，本发明不限于所公开的具体实施例。本发明旨在涵盖所附权利要求的精神和范围内所包括的各种修改和等效布置。
Claims (3)
1.一种基于小波边界元模型的二维正方晶格光子晶体带隙设计方法，其特征在于包括以下步骤：
S1:将区间B样条小波与边界元法相结合，用BSWI尺度函数取代传统边界元的多项式插值，结合单胞技术，获得了二维正方晶格光子晶体基体与散射体统一的离散化边界积分方程形式，进而获得代数方程组；
S2:根据步骤S1得到的代数方程组，在频率域内，结合Bloch定理以及基体与散射体间的连续性条件，进一步建立二维正方晶格光子晶体带隙特性计算模型；
S3:通过调整二维正方晶格光子晶体基体或散射体尺寸，获得所需要的带隙特性，最终完成二维正方晶格光子晶体带隙设计。
2.根据权利要求1所述的一种基于小波边界元模型的二维正方晶格光子晶体带隙设计方法，其特征在于：步骤S1包括以下步骤：
1)运用一维BSWI尺度函数作为插值函数，得离散化的边界积分方程
其中，P和Q分别代表源点和场点，Ne代表单元的个数，c(P)＝β/2π表示与源点P处边界形状有关的系数，β为P处切线张角，u*(P,Q)为本光学问题基本解，q*(P,Q)为基本解沿外法线方向的方向导数，li为单元的长度，
和
分别表示第i个单元的电场/磁场值和其法向导数值所组成的列向量，
指BSWI尺度函数所组成的行向量，Te为本光学问题所对应的转换矩阵，u*(P,Q)，q*(P,Q)，
及Te的表达式分别为：
其中，ε指材料的介电常数，k0＝ω/c0是自由空间波数，c0表示真空中的光速，ω是角频率，r＝|xP-xQ|表示源点与场点的距离，
表示第一类0阶汉克尔函数；
其中，
表示第一类1阶汉克尔函数，xi(Q)与xi(P)分别表示源点与场点的坐标点，ni(Q)表示场点处的方向余弦；
其中，m与j分别表示BSWI尺度函数的阶数和尺度，ξ∈[0,1]为局部坐标；
其中，ξi为第i个节点的局部坐标值，N表示小波单元节点的个数；
将每一个节点都设为源点，经过积分运算、矩阵组装可进一步得到代数方程组：
[H]NP×NP{U}NP×1+[G]NP×NP{Q}NP×1＝0。
其中，H和G为系统矩阵，U和Q分别表示所有节点位移和位移法向导数所组成的列向量，NP为节点的总数。
3.根据权利要求2所述中的所述的一种基于小波边界元模型的二维正方晶格光子晶体带隙设计方法，其特征在于：步骤S2包括以下步骤：
1)通过将计算基体与散射体的系统矩阵的子矩阵进行整合，得
其中，对于TM模式，η＝1，对于TE模式，η＝ε0/ε1，系数η是求解带隙特性的关键，kx,ky为第一布里渊区Bloch波矢，a为晶格常数，
和
表示场点在基体的边界Γ1上，将所有节点都视为源点时积分分别涉及u*(P,Q)和q*(P,Q)所得到的矩阵，其余类似，另外有：
分离出含待求Bloch波矢的项，最终可获得二维正方晶格光子晶体的BSWI边界元带隙计算模型：
AX＝ξBX
其中，
或
2)由于第一布里渊区具有对称性，带隙计算通常沿简约布里渊区的边界进行，在正方晶格简约布里渊区Γ-X-M-Γ的每个边界上涉及到的矩阵A、B为
(1)Γ-X上，
(2)X-M上，
(3)M-Γ上，
根据上述矩阵A、B，可求得正方晶格简约布里渊区Γ-X-M-Γ的每个边界上未知的Bloch波矢，计算完毕之后，用简约波矢M、Γ、X作为横坐标，所求的Bloch波矢作为所在边界横坐标x上的值、归一化频率为ωa/(2πc)纵坐标y，根据归一化频率与所求Bloch波矢的对应关系，就可得到二维正方晶格光子晶体能带结构图，从而获得正方晶格光子晶体带隙特性。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110335850.7A CN113031263B (zh) | 2021-03-29 | 2021-03-29 | 基于小波边界元模型的二维正方晶格光子晶体带隙设计方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110335850.7A CN113031263B (zh) | 2021-03-29 | 2021-03-29 | 基于小波边界元模型的二维正方晶格光子晶体带隙设计方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113031263A true CN113031263A (zh) | 2021-06-25 |
| CN113031263B CN113031263B (zh) | 2022-10-14 |
Family
ID=76452960
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110335850.7A Active CN113031263B (zh) | 2021-03-29 | 2021-03-29 | 基于小波边界元模型的二维正方晶格光子晶体带隙设计方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113031263B (zh) |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1837780A (zh) * | 2006-04-10 | 2006-09-27 | 西安交通大学 | 区间b样条小波单元用于转子横向裂纹定量诊断的方法 |
| CN101114308A (zh) * | 2007-08-24 | 2008-01-30 | 西安交通大学 | 一种三类变量区间b样条小波梁单元构造方法 |
| CN103176272A (zh) * | 2011-12-21 | 2013-06-26 | 北京邮电大学 | 二维光子晶体最大绝对带隙结构优化方法 |
| CN106709202A (zh) * | 2017-01-09 | 2017-05-24 | 温州大学 | 一种基于小波有限元模型的一维声子晶体梁结构带隙设计方法 |
| CN106777771A (zh) * | 2017-01-09 | 2017-05-31 | 温州大学 | 基于小波有限元模型的二维声子晶体板结构带隙设计方法 |
| CN109783946A (zh) * | 2019-01-21 | 2019-05-21 | 河北工业大学 | 一种声子晶体带隙仿真的节点积分算法 |
-
2021
- 2021-03-29 CN CN202110335850.7A patent/CN113031263B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1837780A (zh) * | 2006-04-10 | 2006-09-27 | 西安交通大学 | 区间b样条小波单元用于转子横向裂纹定量诊断的方法 |
| CN101114308A (zh) * | 2007-08-24 | 2008-01-30 | 西安交通大学 | 一种三类变量区间b样条小波梁单元构造方法 |
| CN103176272A (zh) * | 2011-12-21 | 2013-06-26 | 北京邮电大学 | 二维光子晶体最大绝对带隙结构优化方法 |
| CN106709202A (zh) * | 2017-01-09 | 2017-05-24 | 温州大学 | 一种基于小波有限元模型的一维声子晶体梁结构带隙设计方法 |
| CN106777771A (zh) * | 2017-01-09 | 2017-05-31 | 温州大学 | 基于小波有限元模型的二维声子晶体板结构带隙设计方法 |
| CN109783946A (zh) * | 2019-01-21 | 2019-05-21 | 河北工业大学 | 一种声子晶体带隙仿真的节点积分算法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113031263B (zh) | 2022-10-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Kinzel et al. | Kinematic synthesis for finitely separated positions using geometric constraint programming | |
| Liu et al. | Shape optimization of sound barrier using an isogeometric fast multipole boundary element method in two dimensions | |
| Shojaee et al. | Free vibration analysis of thin plates by using a NURBS-based isogeometric approach | |
| Kim et al. | Isogeometric analysis for trimmed CAD surfaces | |
| Venkatakrishnan et al. | Implicit solvers for unstructured meshes | |
| Giraldo et al. | A nodal triangle-based spectral element method for the shallow water equations on the sphere | |
| Yu et al. | Three-dimensional transient heat conduction problems in FGMs via IG-DRBEM | |
| Chen et al. | A unified inverse design and optimization workflow for the Miura-oRing metastructure | |
| Ding et al. | Exact and efficient isogeometric reanalysis of accurate shape and boundary modifications | |
| Prince et al. | Application of proper generalized decomposition to multigroup neutron diffusion eigenvalue calculations | |
| Ying et al. | Fast geodesics computation with the phase flow method | |
| CN109101463A (zh) | 一种广域层状大地格林函数的多精度求解方法 | |
| CN108763777A (zh) | 基于泊松方程显式解的vlsi全局布局模型建立方法 | |
| Haghighi et al. | The fragile points method (FPM) to solve two-dimensional hyperbolic telegraph equation using point stiffness matrices | |
| Yang et al. | An open source, geometry kernel based high-order element mesh generation tool | |
| CN113031263B (zh) | 基于小波边界元模型的二维正方晶格光子晶体带隙设计方法 | |
| Zheng et al. | Three dimensional acoustic shape sensitivity analysis by means of adjoint variable method and fast multipole boundary element approach | |
| CN113050274B (zh) | 基于小波边界元模型的三角晶格声子晶体带隙设计方法 | |
| Liu et al. | Review of subdivision schemes and their applications | |
| Roy et al. | A finite volume method for viscous incompressible flows using a consistent flux reconstruction scheme | |
| Wolkov | Application of the multigrid approach for solving 3D Navier-Stokes equations on hexahedral grids using the discontinuous Galerkin method | |
| Weiyuan et al. | A feature points-based method for data transfer in fluid-structure interactions | |
| Zhang et al. | Isogeometric dual reciprocity BEM for solving time-domain acoustic wave problems | |
| Marshall et al. | A three-dimensional, p-version BEM: High-order refinement leveraged through regularization | |
| Wang et al. | Acoustic problems analysis of 3D solid with small holes by fast multipole boundary face method |
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 | ||
| EE01 | Entry into force of recordation of patent licensing contract |
Application publication date: 20210625 Assignee: ZHEJIANG GREATWALL MIXERS CO.,LTD. Assignor: Wenzhou University Contract record no.: X2023330000104 Denomination of invention: Design Method for Band Gap of Two-dimensional Square Lattice Photonic Crystals Based on Wavelet Boundary Element Model Granted publication date: 20221014 License type: Common License Record date: 20230311 |
|
| EE01 | Entry into force of recordation of patent licensing contract |