CN102073796B - 一种模拟溶质三维运移过程的格子行走方法 - Google Patents

一种模拟溶质三维运移过程的格子行走方法 Download PDF

Info

Publication number
CN102073796B
CN102073796B CN201110041567XA CN201110041567A CN102073796B CN 102073796 B CN102073796 B CN 102073796B CN 201110041567X A CN201110041567X A CN 201110041567XA CN 201110041567 A CN201110041567 A CN 201110041567A CN 102073796 B CN102073796 B CN 102073796B
Authority
CN
China
Prior art keywords
delta
solute
lattice
concentration
partiald
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.)
Expired - Fee Related
Application number
CN201110041567XA
Other languages
English (en)
Other versions
CN102073796A (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.)
Nanjing University
Original Assignee
Nanjing 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 Nanjing University filed Critical Nanjing University
Priority to CN201110041567XA priority Critical patent/CN102073796B/zh
Publication of CN102073796A publication Critical patent/CN102073796A/zh
Application granted granted Critical
Publication of CN102073796B publication Critical patent/CN102073796B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Abstract

本发明公开了一种模拟三维溶质运移过程的格子行走方法，该方法首先投放溶质源；确定边界条件，设定网格间距
Figure DEST_PATH_IMAGE002
、
Figure DEST_PATH_IMAGE004
、
Figure DEST_PATH_IMAGE006
和时间间隔
Figure DEST_PATH_IMAGE008
；根据流场流速和弥散系数确定转移矩阵
Figure DEST_PATH_IMAGE010
，在
Figure DEST_PATH_IMAGE012
时刻格点
Figure DEST_PATH_IMAGE014
处的浓度
Figure DEST_PATH_IMAGE016
，由前一时刻的浓度就可以求出后一时刻的浓度。如果弥散系数
Figure DEST_PATH_IMAGE018
、
Figure DEST_PATH_IMAGE020
或随空间位置连续变化，则流速各分量分别修正为，
Figure DEST_PATH_IMAGE026
，
Figure DEST_PATH_IMAGE028
，然后用修正的速度分量求出转移矩阵
Figure 386123DEST_PATH_IMAGE010
。通过对三维定通量连续注入点源在均匀介质中的输运模拟和瞬时投放点源在随空间位置连续变化的介质中的溶质运移模拟，发现该方法模拟得到的结果与解析解以及有限差分方法的结果吻合得很好。该方法最大的优点是可以在给定的网格间距下模拟任意流速的溶质运移，因此它特别适用于模拟对流占优问题。

Description

一种模拟溶质三维运移过程的格子行走方法
技术领域
本发明涉及水力学领域，具体涉及一种模拟溶质三维运移过程的格子行走方法。
背景技术
水环境问题是当前和人类生存关系至为密切的一个重要问题。随着经济的发展和人口的增多，用水量越来越大，而水资源的污染和破坏使得用水矛盾更加突出。在我国北方地区用水大部分采自地下水，因此对地下水污染的监测和控制有重要的意义。地下水中的污染物大多以溶质的形式存在，地下水中的溶质运移过程通常由对流弥散方程描述，其三维形式为：
Figure 201110041567X100002DEST_PATH_IMAGE001
  (1)，
这里
Figure 872754DEST_PATH_IMAGE002
是溶质浓度，、
Figure 539358DEST_PATH_IMAGE004
、
Figure DEST_PATH_IMAGE005
分别为x、y和z方向的弥散系数，
Figure 618173DEST_PATH_IMAGE006
、
Figure DEST_PATH_IMAGE007
和
Figure 306031DEST_PATH_IMAGE008
分别为x、y和z方向的流速。一般地说弥散系数张量还有
Figure DEST_PATH_IMAGE009
等非主轴项。在局域坐标系里，只要我们把水流方向设为坐标轴，非主轴项就变为0。对流弥散方程可以用常规的有限差分方法和有限单元方法求解。但有限差分方法处理对流占优问题时的效率比较低，如在x方向上，它要求
Figure 749781DEST_PATH_IMAGE010
，如果较大，则格点间距
Figure 152950DEST_PATH_IMAGE012
需要取得比较小。同时这些方法常常会遇到两个困难，即数值弥散和数值振荡。产生这种不良数值现象的原因是对流扩散方程同时包含双曲项和抛物项。[薛禹群，谢春红  地下水数值模拟 北京：科学出版社，2007]。为了克服这个困难，人们做了大量的努力，发明了很多方法如上风方法、特征线方法、随机行走方法等。上风方法可以减少数值振荡，但同时会增加数值弥散[N. Sun, and W. Yeh, A proposed upstream weight numerical method for simulating pollutant transport in ground water, Water Resour. Res. 19 1489 (1983).]，且在复杂流场中上风方法的应用会变得很麻烦。特征线方法是比较常用的有效方法，它用拉格朗日法处理对流项，用欧拉法处理扩散项[S. P. Neumann, Adaptive Eulerian Lagrangian finite-element method for advection dispersion,Int. J. Numer. Meth. Eng. 20 321 (2003)]，但缺点是不能保证溶质的质量守恒。随机行走方法是典型的拉格朗日方法[A. Tompson and L. w. Gelhar, Numerical simulation of solute transport in three-dimensional,randomly heterogeneous porous media, Water Resour. Res. 26 2541 (1990), E. M. LaBolle, et al, Random-Walk Simulation of Transport in Heterogeneous Porous Media:Local Mass-Conservation Problem and Implementation Methods, Water Resour. Res. 32 583  (1996)], 它消除了数值弥散，但它的适用条件是流场中流速和扩散系数必须变化缓慢，且需要足够多的粒子数来表示质量分布。如果边界条件比较复杂，也会给随机行走方法带来很大的麻烦。
发明内容
发明目的：针对现有技术中存在的不足，本发明的目的是提供一种模拟溶质三维运移过程的格子行走方法。在弥散系数张量不变的情况下，该方法模拟得到的结果与解析解非常吻合。在弥散系数随空间位置连续变化的情况下，该方法模拟得到的结果与有限差分方法计算的结果一致。
技术方案：为了实现上述发明目的，本发明采用的技术方案为：
一种模拟溶质三维运移过程的格子行走方法，包括以下步骤：
（1）在
Figure DEST_PATH_IMAGE013
的时刻，投放溶质源；
（2）确定边界条件，设定网格长度
Figure 973138DEST_PATH_IMAGE012
和时间间隔
Figure 325622DEST_PATH_IMAGE014
；
（3）在
Figure DEST_PATH_IMAGE015
时刻格点
Figure DEST_PATH_IMAGE017
处的浓度
Figure 455121DEST_PATH_IMAGE018
，
Figure DEST_PATH_IMAGE019
为转移矩阵；先考虑t时刻在格点
Figure 830739DEST_PATH_IMAGE020
的粒子经过
Figure 438307DEST_PATH_IMAGE014
时间的运动达到的平均位置，格点是距离此平均位置最近的格点，
Figure 848559DEST_PATH_IMAGE022
定义为，
Figure 696430DEST_PATH_IMAGE024
定义为
Figure DEST_PATH_IMAGE025
，
Figure 823259DEST_PATH_IMAGE026
定义为：
Figure DEST_PATH_IMAGE027
，定义自变量为整数的函数
Figure 985250DEST_PATH_IMAGE028
、
Figure DEST_PATH_IMAGE029
和
Figure 499277DEST_PATH_IMAGE030
，
，
Figure 721311DEST_PATH_IMAGE032
，
Figure DEST_PATH_IMAGE033
，
则
Figure 320789DEST_PATH_IMAGE034
；这样由前一时刻的浓度就可以求出后一时刻的浓度；
（4）若介质为非均质介质，
Figure 83208DEST_PATH_IMAGE003
、
Figure 264791DEST_PATH_IMAGE004
、
Figure 657726DEST_PATH_IMAGE005
连续变化，则在求解时采用的速度分别为：
Figure DEST_PATH_IMAGE035
，，
Figure DEST_PATH_IMAGE037
；
（5）以此迭代，则可以求出任意时刻的浓度，得出结果。
步骤（1）中，所述的溶质源是连续均匀投放的点源。
步骤（2）中，网格长度
Figure 310610DEST_PATH_IMAGE012
和时间间隔，优先使其满足条件：
Figure 910536DEST_PATH_IMAGE038
,
Figure 487535DEST_PATH_IMAGE003
、
Figure 591758DEST_PATH_IMAGE004
、
Figure 685615DEST_PATH_IMAGE005
分别是x、y、z 方向的弥散系数分量，这种间隔选取可模拟任意流速的溶质运移；
所述的弥散系数分量
Figure 482670DEST_PATH_IMAGE003
连续变化，
Figure DEST_PATH_IMAGE039
。
有益效果：本发明提出了一种全新的模拟溶质三维运移过程的方法，它将溶质随水的移流过程与溶质的扩散过程分开考虑，即移流决定平均位置，扩散决定从平均位置向附近格点的运动。通过0-2阶矩相等，用平均位置附近格点上的溶质的分立分布表示连续空间中的高斯分布。通过对连续均匀投放的点源在三维均质介质中的运移模拟，发现该方法模拟得到的结果与解析解吻合得很好。然后在非均质介质中瞬时投放点源，通过速度修正，得到的模拟结果与有限差分方法得到的结果一致。格子行走方法最显著的特点是格点距离的选择可以与流场速度无关，在处理对流占优问题时，效率比有限差分方法高很多，因为有限差分方法在计算时要求格点的间隔满足
Figure 278457DEST_PATH_IMAGE040
、
Figure DEST_PATH_IMAGE041
和
Figure 389632DEST_PATH_IMAGE042
。
附图说明
图1是粒子从
Figure 400314DEST_PATH_IMAGE020
点到
Figure DEST_PATH_IMAGE043
邻近格点处的运动轨迹。
图2是连续注入点源在三维均质介质中运移的溶质分布等浓度线图。
图3是连续注入点源在三维均质介质中运移的溶质分布剖面图。
图4是瞬时投放点源在三维非均质介质中运移时的溶质分布等浓度线图。
图5是瞬时投放点源在三维非均质介质中运移时的溶质分布剖面图。
具体实施方式
如图1所示，t时刻在格点的粒子经过时间的移流与扩散，平均位置在处。离平均位置最近的格点记为。由于弥散作用，粒子可以运动到
Figure 491635DEST_PATH_IMAGE021
及其附近的格点上。我们假定粒子只能弥散到以格点
Figure 12746DEST_PATH_IMAGE021
为中心，边长为
Figure 262462DEST_PATH_IMAGE044
、
Figure DEST_PATH_IMAGE045
、
Figure 437616DEST_PATH_IMAGE046
的立方体上的所有格点上，包括格点，共27个格点。
Figure 755782DEST_PATH_IMAGE022
定义为
Figure 809189DEST_PATH_IMAGE023
，
Figure 835919DEST_PATH_IMAGE024
定义为
Figure 316579DEST_PATH_IMAGE025
，
Figure 609020DEST_PATH_IMAGE026
定义为
Figure 403801DEST_PATH_IMAGE027
，运动到格点
Figure 98087DEST_PATH_IMAGE017
的几率记为。
先考虑弥散系数为常数时的情况，由方程(1)，可以得到格点
Figure 936599DEST_PATH_IMAGE020
上的粒子经过时间的运动后的分布几率：
 2)
根据质量守恒，分立分布
Figure 863601DEST_PATH_IMAGE047
的
Figure DEST_PATH_IMAGE049
阶矩应该与连续分布
Figure 138594DEST_PATH_IMAGE050
相等，得到：
                                          （3）
根据平均位置相等，则分立分布的一阶矩和连续分布相等，可以得到：
        
Figure 343310DEST_PATH_IMAGE052
                   （4）
由分立分布和连续分布的二阶矩相等，得到：
Figure DEST_PATH_IMAGE053
    （5），
考虑到（2）式中的
Figure 466511DEST_PATH_IMAGE050
可以写为x、y、z方向上的独立分布乘积，先考虑各个方向的分布，定义自变量为整数的函数(要求格点
Figure DEST_PATH_IMAGE055
满足
Figure 332016DEST_PATH_IMAGE056
)。由x方向0-2阶矩相等，则：
Figure DEST_PATH_IMAGE057
,
可以得到：
                 
Figure 7717DEST_PATH_IMAGE031
，                  （6）
类似地，定义并得到
Figure 10308DEST_PATH_IMAGE029
、
Figure 471376DEST_PATH_IMAGE030
，表达式与（6）类似，、和
Figure 604920DEST_PATH_IMAGE012
要作相应的替换。于是可以得到：
                ，                （7）
它满足(3)-(5)式，同时还使得分立分布与连续分布的部分三阶矩和四阶矩相等。浓度随时间的变化写为转移矩阵的形式：
Figure DEST_PATH_IMAGE059
，
（7）式中的
Figure 521241DEST_PATH_IMAGE047
即为
Figure 171534DEST_PATH_IMAGE019
，格点
Figure 515928DEST_PATH_IMAGE020
决定(7)中的格点
Figure 951588DEST_PATH_IMAGE021
的坐标。几率必须为非负数，否则负的转移几率没有物理意义。
    对于任意流速，只要满足：
          
Figure 852035DEST_PATH_IMAGE060
  ，                         （8）
Figure 734540DEST_PATH_IMAGE047
不会出现负值。这种网格破分不依赖于流速。所做的大量模拟也表明这种网格破分可以模拟任意流速下的溶质运移而不会出现明显的数值弥散和数值振荡。
如果污染物在弥散系数连续变化的非均质介质中传播，需要做一些修正。在推导转移矩阵
Figure 87024DEST_PATH_IMAGE019
时，用到了连续空间的几率分布（2），而这个分布函数是下面的微分方程：
Figure DEST_PATH_IMAGE061
的解[Kinzelbach W, Groundwater Modeling: an Introduction With Sample Programs in BASIC, Elsevier. (1986).]，此偏微分方程右边比对流弥散方程（1）多出一项
Figure 701676DEST_PATH_IMAGE062
，
把它们减去并与移流项合并，就会得到，
Figure 12758DEST_PATH_IMAGE036
和，然后以这些修正的速度代入（6）式和（7）式计算转移矩阵。
实施例1
以下结合具体实施例验证本发明的效果。
先以三维定通量连续投入点源模型验证弥散系数为常数的格子行走方法，并与解析解结果做对比。从
Figure 270881DEST_PATH_IMAGE013
的时刻起，在单位时间内向原点处连续投放单位质量的污染源。本例中
Figure DEST_PATH_IMAGE063
，
Figure 383062DEST_PATH_IMAGE064
，，，
Figure DEST_PATH_IMAGE067
，
Figure 186644DEST_PATH_IMAGE068
。格点间距选为
Figure DEST_PATH_IMAGE069
，时间间隔取
Figure 408678DEST_PATH_IMAGE070
。虽然y和z方向的间距不满足（8），但是
Figure 555625DEST_PATH_IMAGE047
在本例中的流速下不会出现负值。图2和图3是在t=500时的溶质分布图。其中图2给出的是
Figure DEST_PATH_IMAGE071
时x-y平面的等浓度线，图3描述的是y、z为常数时，浓度随x的变化。从这两个图中我们可以看到格子行走方法给出的结果与解析解符合得很好。
实施例2
再考虑弥散系数连续变化时的模拟。令
Figure 770575DEST_PATH_IMAGE039
，
Figure 952157DEST_PATH_IMAGE072
，
Figure DEST_PATH_IMAGE073
，
Figure 345093DEST_PATH_IMAGE074
，，
Figure 997977DEST_PATH_IMAGE065
。格点间距选为
Figure 34066DEST_PATH_IMAGE069
，时间间隔取
Figure 597902DEST_PATH_IMAGE070
。转移矩阵不会出现负值。在
Figure 719442DEST_PATH_IMAGE013
时向原点投入单位质量的瞬时溶质点源。图4和图5描述的是
Figure DEST_PATH_IMAGE075
时的溶质分布图。其中图4描述的是
Figure 276194DEST_PATH_IMAGE071
的x-y平面的等浓度线，图5描述的是y、z为常数时，浓度随x的变化。从图4和图5中我们可以看到格子行走方法给出的结果与有限差分方法给出的结果非常接近。因此，格子行走方法适用于非均质介质中的溶质输运模拟。

Claims (3)

1.一种模拟溶质三维运移过程的格子行走方法，其特征在于，包括以下步骤：
（1）在t=0的时刻，投放溶质源；
（2）确定边界条件，分别设定x，y，z方向的网格间距Δx、Δy、Δz，设定时间间隔Δt；
（3）在t+Δt时刻格点i处的浓度
Figure FDA00001960553100011
Kij为转移矩阵，Cj(t)为t时刻格点j处的溶质浓度；先考虑t时刻在格点j的溶质分子经过Δt时间的运动达到的平均位置；格点j的坐标记为(xj,yj,zj)，则溶质分子的平均位置为(xj+uxΔt,yj+uyΔt,zj+uzΔt)，其中ux、uy和uz分别为水流速度在x、y和z方向的分量；格点m是距离此平均位置最近的格点，rx定义为xj+uxΔt-xm，ry定义为yj+uyΔt-ym，rz定义为zj+uzΔt-zm；弥散张量D在x、y和z方向的分量记为Dxx、Dyy和Dzz；定义自变量为整数的函数px(n)、py(n)和pz(n)，
px ( - 1 ) = ( 2 D xx Δt + r x 2 - r x Δx ) / ( 2 Δ x 2 ) px ( 0 ) = ( - 2 D xx Δt + Δx 2 - r x 2 ) / Δ x 2 px ( 1 ) = ( 2 D xx Δt + r x 2 + r x Δx ) / ( 2 Δx 2 ) px ( n ) = 0 , n ≠ ± 1,0 ,
py ( - 1 ) = ( 2 D yy Δt + r y 2 - r y Δy ) / ( 2 Δ y 2 ) py ( 0 ) = ( - 2 D yy Δt + Δy 2 - r y 2 ) / Δ y 2 py ( 1 ) = ( 2 D yy Δt + r y 2 + r y Δy ) / ( 2 Δy 2 ) py ( n ) = 0 , n ≠ ± 1,0 ,
pz ( - 1 ) = ( 2 D zz Δt + r z 2 - r z Δy ) / ( 2 Δ z 2 ) pz ( 0 ) = ( - 2 D zz Δt + Δz 2 - r z 2 ) / Δ z 2 pz ( 1 ) = ( 2 D zz Δt + r z 2 + r z Δz ) / ( 2 Δz 2 ) pz ( n ) = 0 , n ≠ ± 1,0 ,
则 K ij = px ( x i - x m Δx ) py ( y i - y m Δy ) pz ( z i - z m Δz ) ; 这样由前一时刻的浓度就可以求出后一时刻的
浓度；
（4）如果溶质在弥散系数连续变化的非匀质介质中传播，则在求解时采用的速度分量分别为： u x = u x + ∂ D xx / ∂ x , u y = u y + ∂ D yy / ∂ y ,
u z = u z + ∂ D zz / ∂ Z ;
（5）以此迭代，则可以求出任意时刻的浓度，得出结果。
2.根据权利要求1所述的模拟溶质三维运移过程的格子行走方法，其特征在于：步骤（1）中，所述的溶质源是连续均匀投放的点源。
3.根据权利要求1所述的模拟溶质三维运移过程的格子行走方法，其特征在于：步骤（2）中，网格间距和时间间隔Δt，满足条件：
8 3 D xx Δt ≤ Δx ≤ 8 D xx Δt 8 3 D yy Δt ≤ Δy ≤ 8 D yy Δt 8 3 D zz Δt ≤ Δz ≤ 8 D zz Δt .
CN201110041567XA 2011-02-21 2011-02-21 一种模拟溶质三维运移过程的格子行走方法 Expired - Fee Related CN102073796B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201110041567XA CN102073796B (zh) 2011-02-21 2011-02-21 一种模拟溶质三维运移过程的格子行走方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201110041567XA CN102073796B (zh) 2011-02-21 2011-02-21 一种模拟溶质三维运移过程的格子行走方法

Publications (2)

Publication Number Publication Date
CN102073796A CN102073796A (zh) 2011-05-25
CN102073796B true CN102073796B (zh) 2012-11-21

Family

ID=44032335

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201110041567XA Expired - Fee Related CN102073796B (zh) 2011-02-21 2011-02-21 一种模拟溶质三维运移过程的格子行走方法

Country Status (1)

Country Link
CN (1) CN102073796B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102509009B (zh) * 2011-11-02 2014-12-17 赵健伟 一种基于受限空间内随机行走的仿真方法
CN104865165B (zh) * 2015-06-05 2018-01-12 武汉大学 全过程通量守恒的现场土壤溶质运移弥散系数测定的方法
CN105512417B (zh) * 2015-12-17 2018-09-18 中国环境科学研究院 基于粒子追踪的孔隙地下水污染物三维运移模拟方法
CN108614942B (zh) * 2018-05-10 2020-05-05 河海大学 一种表征多孔介质中溶质运移时空尺度关联的方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101866392A (zh) * 2010-05-21 2010-10-20 南京大学 一种模拟溶质一维运移过程的格子行走方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8147615B2 (en) * 2004-11-05 2012-04-03 Infineon Technologies Ag Method of fabricating semiconductor cleaners

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101866392A (zh) * 2010-05-21 2010-10-20 南京大学 一种模拟溶质一维运移过程的格子行走方法

Also Published As

Publication number Publication date
CN102073796A (zh) 2011-05-25

Similar Documents

Publication Publication Date Title
Pollock User guide for MODPATH Version 7—A particle-tracking model for MODFLOW
Sotiropoulos et al. Immersed boundary methods for simulating fluid–structure interaction
Mirzaei CFD modeling of micro and urban climates: Problems to be solved in the new decade
CN102945295B (zh) 一种格子玻尔兹曼方法的并行加速方法及系统
Morgenthal et al. An immersed interface method for the vortex-in-cell algorithm
CN102073796B (zh) 一种模拟溶质三维运移过程的格子行走方法
CN107341318B (zh) 一种基于全河的月径流时间位移二维矩阵的模拟方法
CN106651765A (zh) 一种采用深度神经网络的缩略图自动生成的方法
CN110119590B (zh) 一种基于多源观测数据的水质模型粒子滤波同化方法
CN105512417A (zh) 基于粒子追踪的孔隙地下水污染物三维运移模拟方法
Yue et al. A multi‐grid method of high accuracy surface modeling and its validation
CN104035096B (zh) 一种基于多普勒天气雷达的垂直风廓线非线性反演方法
CN107657075A (zh) 模拟地下水介质交界面处达西速度的区域分解有限元法
CN110990926B (zh) 一种基于面积修正率的城市地表建筑水动力学仿真方法
CN102692491A (zh) 基于分阶段禁忌寻优算法的土壤水分特征参数计算方法
CN105974471B (zh) 一种基于异步流的地震数据多gpu快速正演计算方法
CN103077274A (zh) 高精度曲面建模智能化方法及装置
Zhang et al. The CA model based on data assimilation
CN109408838A (zh) 一种对缝洞型油藏剩余油进行快速分析的方法及系统
CN107832482A (zh) 致密储层多尺度裂缝网络建模及模拟方法
CN101866392B (zh) 一种模拟溶质一维运移过程的格子行走方法
Li et al. An improved non-iterative surface layer flux scheme for atmospheric stable stratification conditions
Ghisu et al. A level-set algorithm for simulating wildfire spread
CN102094402B (zh) 一种模拟溶质二维运移过程的格子行走方法
CN104599318B (zh) 一种植物三维模型网格无缝融合的方法及系统

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20121121