見る/開く - 琉球大学

Title
Author(s)
Citation
Issue Date
URL
Rights
沿岸開領域における非線形波動解析のための新しい無限
要素(第2報)−高次無限要素と誤差評価−
筒井, 茂明
琉球大学工学部紀要(59): 39-51
2000-03
http://ir.lib.u-ryukyu.ac.jp/handle/123456789/2220
琉球大学工学部紀要第59号,2000年
39
沿岸開領域における非線形波動解析のための新しい
無限要素(第2報)-高次無限要素と誤差評価一
筒井茂明*
NewlnfiniteElementfOrNonnnearWaveAnalysesinUnboundedCoastalDomains
Part2.-HigherOrderlnfiniteElementandErrorEstimateShigeakiTsuTsuI*
Abstract
Tieatmentoftheradiationboundaryconditionatinfinityisofgreatsignificancefbrwaveanalysesinunbounded
coastaldomains・TYlehigherorderinfiniteelementfbrthefiniteelementmethodisnewlydevelopedinorderto
improveaccuracyratherthanthatbythelinearinfiniteelementpresentedinthisserieSofreseamhes・Theshape
fimctionsfbrfUncUonvariatiomarethird-orderpolynomialsintheradialcoordinateandthevariationintheangular
coordinateislinear,Theelementhasthreenovelfeatures・First,themappingfbrfimctionvariationisableto
perfbnntheradiationboundalyconditionuptor-”,rbeingtheradialdistanceandm=l/ZfbrwaveprOblems・The
secondfeamleiseasinessinConnectingtheinfmiteelementtotheintelior,standardelements、T11ethirdoneisan
integrablefeatureovertheelementwiththeaidoftheinfinitemapping・Asatestofthenew,third-orderinfinitc
element,theproblemofwavediffisactionbyacylinderissolved,alongwithacompansonwiththeanalytical
solutionandwithnumericalresultsobtainedusmgotherfinitcelements・Numcricalerrorassessmcntshowsefficiency
ofthenewinfiniteelemenL
KeyWords:Infiniteelement,Finiteelementmethod,Unboundedwaveproblems,Nonlinearwaves,Mapping.
1.
緒言
無限遠での放射境界条件の処理は開領域での波動場の解
析にとって非常に重要である.前報(筒井,1999a)におい
ては,無限要素を用いた有限要素法による統一的な取り扱
いのため,開領域での非線形な波動場解析をも視野に入れ
て,従来の無限要素に代わる新しい無限要素が提案された.
ざらに,その適用性について数値計算例による検証が行わ
れ,次のような重要な特性を持つことが示された.
(a)無限要素と解析対象領域内の要素との結合は簡単で,
式はHeImholtzの方程式で,ハンケル関数はその基本解で
ある.ハイブリッド型の有限要素法(I形EM)(Chen・Mei,
1975)では,遠方場での波をハンケル関数を用いた級数解
で表し,これを境界解とし内部の有限要素解と接続する.し
たがって,級数解の項数を十分多く採ることにより高精度
な波動場解析が可能であり,その最終的な推算精度は一義
的に要素サイズに依存する.
前報で提案された無限要素は前述の(a),(b)のような優れ
た特性を持っているが,その精度はI盃EMの計算精度と比
かつ内部要素の形状に左右されない.
較してわずかに劣る.そこで本研究では,取り扱いの簡便
さなどの本来の特性を保ちつつ,波高推算の高精度化を図
(b)無限要素の要素積分においては,被積分関数に強い
特異性が存在せず,積分が陽に表示される.したがって,
るため,高次の無限要素を開発・提案する.さらに,数値
計算例により推算精度を調べ,その有用性を示す.
数値積分は容易で精度良<計算される.
無限要素が適用される外部領域での波動場の支配方程
受理:19,9年12月6日
*琉球大学工学部現u魔建設工学科
Dept・ofCivilEmgineeringandA1℃hitectur巳,Facul[yofEngrg
2.無限要素を用いた弱形式
解析対象海域を,図一1のように,仮想境界Tlにより領域。,
およびQ2に分割する.内部領域。lでの水深は変化し,そこ
には海岸線や構造物などの境界「Bあるいは水深不連続部r、
40
筒井:沿岸開領域における非線形波動解析のための新しい無限要素(第2報)-高次無限要素と誤差評価一
ite
ent
鶚lL仲w…蟄岬ル.:
-1.V"州"WD"
薑伽Mw“
(7)
さらに,式(7)の各項を離散化して得られる要素行列をそれ
図-1解析対象海域と定義
が存在する.外部領域Q2での水深は一定と仮定され,無限
要素はこの外部領域に適用される.
本章と次章における全ての物理量は,代表水深h6および
重力加速度gにより定められる基準長:h6,時間:↑/75両,速
度:、/詞による無次元量である.静水面に座標原点を置き,
水平方向に(jr,y)軸,鉛直上方にz軸を採る.また,(4m:領
域Q,およびQ2での水面変動量,刀G:入射波および反射波な
どの幾何光学的な波,が:散乱波,c:波速,cg:群速度,の:周波
数とする.
支配方程式は次式で与えられる.
が得られる.
z(K)(蜥禺(KW)
Ql
-昂(K圏)(由)-墨(K、)(`鰹)薑lli(9)(8)
支配方程式(2)は外部領域Q2における散乱波のみに対する
ものであるから,式(8)においては境界Tl上での波動量が不
連続となっている.そこで境界条件,式(3.2)を満たすため,
次のような処理を行う.無限要素の要素行列(K-)および境
界「,上の節点での幾何光学的な波(〃G}より{K~){〃c}なる
量を求め,式(8)の両辺の対応する節点に付加すると,次式
が得られる.
V(cc圏V勾十⑳2(c圏/c)`1W=g領域QlW(1)
v(・・愚W脆の2(。gノc)"`=o領域Q2内(2)
境界条件は次の通りである.
〃1.Vい"2.V(〃。+〃`)=0境界T1上(31)
畠=ワ。+刀s
ぞれ{K),(K~MKB),(K、}'(9)とすると,次の連立方程式
境界Tl上(3.2)
恩(k)(蜥禺(x~)(〃w)
一息(駈圃)(田)-昂(x・)(畠)
=胃(9)十部K~)(可.)
(9)
したがって,境界Tl上では全体の波動量ら=りG+77sが未知
、B、V畠=B魚境界TB上(4)
量となり,式(3.2)の境界条件が満される.
["ハv鳥]:i-咄鰯「・上(3)
る無限要素による波動場の定式化の際に,近似的ではある
無限遠でのSommerfUdの放射条件(6)
ただし,V=(3/、X,ヨノa”〃=±1,±2,…であり,(、1,,2),、B,
、Dはそれぞれ境界『1,TB,TDでの外向き法線を表す゛
式(1)は,0次のFourier成分波に対する非線形波動方程式で
あるく筒井・大木,1998;Tsutsuiら,1918).同式の左辺は緩
勾配方程式,右辺のQは高次成分波を含む非線形項を表す.
式(2)は一定水深のときにはHelmholtz方程式となる.式(4)
は反射境界条件(Tsutsui・Lewis,1992)であり,係数Bは
境界での波の反射率および境界への波の入射角の関数であ
る.水深不連続部における境界条件,式(5)は,水深の深い
Zamami,1993).ただし,50:水深不連続部の位置,D:無次
元係数である.なお,式(1)-(6)においては,Qを除く波動
量は全て"次成分波に対するものであるが,曲以外の添字〃
が自動的に満たされる.
3.高次の無限要素
(1)1次元無限要素
図-2に示すx座標系において,節点1,2で表される半無限
区間[x,,。。]から成る1次元無限要素を考える.撹乱現象の中
心点を0とし,区間[え。,ェ,]を無限要素に対する補要素と定義
する.無限要素内の1点と中心点との距離はr=x-xoであ
り,物理変数、が遠方場で「-腕の形で減衰することが判って
00功
側より浅い側へ波が伝播すると仮定して適用する(Tsutsui・
無限遠での放射境界条件,式(6)は,外部領域Q2におけ
1
3
2
エ
Xl
為魂= ̄
(a)1次元無限要素
は省略されている.
式(1),(2),(3.1),(4),(5)に対する弱形式は,形状関数γを用
いると,次式で与えられる.
lL("・wvい。仙仙)`。,
32
--o-し6
-101
(b)正規ご系
図-21次元無限要素と正規§系
琉球大学工学部紀要第59号,2000年
いるものとする.ただし,’'nは正の実数である.この無限要
素が区間[-1,1]の正規官系に写像できると仮定すると,節点
1,2はそれぞれ正規5系におけるご=-1およびご=]に対応す
る.x座標系と正規§系の間の写像,すなわち,無限写像は次
式で与えられる(菊池・岡部,1986i筒井,1999a,b).
‘-…-叩(¥「鶚
(10)
上式は補要素上の節点座標により定められていることが重
要である.
さらに,高次要素のために,5=Oに対応する位置に節点3
を採ると,その位置はx3-xo=2淵に,-J[o)となる.従来の無
限要素(Bettesら,1984;Zienkiewiczら,1985)は腕=1の場合
に対応し,そこでは所要の減衰モードが得られるように
種々の工夫がなされている.
ここで,中心点と節点lの間の距離「,=x1-xoによる無次
元距離をp=r/「,と定義すると,式(10)より次式が得られる.
れる.同様に,節点3では§=0,p‐腕=l/2,Ⅳ,=0,M=lで
あるから,β,=1ノ2,β2=】/8となる.したがって,形状関数
ⅣIw3は次式で与えられる.
M一十(\-4(\11
隅(\-(\))
は3節点を用いて決定されているにもかかわらず,その補
間多項式の次数は3次であり,O(p-5噸)の項は無視されてい
る.図-3はこれら3次の形状関数を示す.より高次の形状
関数は,中間節点を補い,式(14)と同様に高次モードの1次
結合を仮定して,簡単に定めることができる.
1
すなわち,物理変数の減衰特性r-mをp-mのモードで代替え
0.5
物理変数の遠方場での減衰特性として,その漸近展開に
0
おいてp-m,p-2m,p-3m,…のように1,2,3,…次と全ての項が
現れる現象では,高次の無限要素を求めるときに,有限要素
法で用いられる多項式補間を試行関数の定義に用いること
ができる.しかし,波動問題においては,以下のように特別
な配慮が必要である.
一定水深における波動場を表すHelmholtz方程式のグリー
ン関数は0次のハンケル関数Hoであり,その漸近展開は次
式で与えられる.
川…(趣(邊一:)}馬,に)
Pに)薑'十金一念士。。(z-3)
)
I
(15)
ただし,M=l-NI-Mである.式(15)で表される形状関数
LL子i'=p-"('1)
させることになる.
41
~--
-0.5
-1
-05
0
0.551
図-31次元無限要素に対する3次の形状関数
(2)2次元無限要素
図-4に示すように,有限要素および無限要素がそれぞれ
適用される領域。,とQ2との境界TI上における内部要素の節
点を1,2とする.無限要素Q;はこの有限辺と放射状の2直
(12)
線とで構成される.放射状の2直線の交点を中心点0と呼
ぶ.さらに,三角形012を無限要素Q;に対する補要素Qfと
ただし,複合(±)は第1,2種のハンケル関数に対応し,iは
虚数単位(‘=T)である.したがって,漸近展開の主要項か
ら,波動問題では係数m=1ノ2であり,散乱波は遠方場にお
y
【】
いてp-噸,p-3'",p-5m,…のように奇数次の項のみで表される
減衰特性を持っている.しかし,通常の形状関数にはこのよ
うな特性をもつ補間式は存在しないので,所要の形状関数
(
を新たに求める必要がある.
ここでは,物理変数に対する試行関数。hを節点1,3での形
状関数Ⅳ1,Mを用いて次式で表す.
‘h=。,Ⅳ,+俺M
(13)
図-43次の無限要素とその補要素
2(-1,1)4(0,1)6(1,1)
ただし,無限遠ではルーoとなることが考慮されている.こ
の試行関数によりp-mならびにp-3mの項が再現可能なよう
に,両者の1次結合でⅣ,,Ⅳ3を定める.すなわち,
に:;蝋#:)(噸
と固き,係数αj,β(i=1,2)を決定する.
節点1ではど=-1,p-m='であり,形状関数の満たすべき
特性としてⅣ,=1,Ⅳ3=Oであるから,係数α,=α2='が得ら
1(-1,-1)3(0,-1)5(1,-1)
図-5正規(5m)系
筒井:沿岸開領域における非線形波動解析のための新しい無限要素(第2報)-高次無限要素と誤差評価一
42
定義し,中心点(xb,yo)からの距離を「とする.無限要素。;は
図-5に示す正規(§,刀)系に写像できると仮定する.
波動問題における散乱波の振幅のように,減衰が散乱中
心からの距離,すなわち,動径方向にのみ生じる場合に
は,正規(§,刀)系においては5方向のみの変化となる.した
がって,ど方向には〆鰄およびp-3鰯で表される減衰項を満足
させることを考える.一方,〃方向には,内部要素との結合
の簡便さを保つため,物理変数の線形変化を仮定する.こ
のように,(師)方向にそれぞれ次数の異なる特殊な高次無
限要素を導く.
1次元の場合と同様に,正規(§,刀)系の2点(0,-1)および
(0,1)に対応する新たな2節点3,4を無限要素の放射線上に設
定する.無限要素の正規(§,〃)系への写像,すなわち,無限
写像は,補要素Qiの節点0,1,2のみで決定され,前報と同
様に次式(菊池・岡部,1986;筒井,1999a)で与えられる.
ii
MM
Cl00
xvジ
、Ⅲ、W
仁0
2図引2図引
I7
--一一
功光
工v〆
)
(16.1)
脇-号(学`(\)1旱。噸(-i鰯(パ\「鶚)
(19)
ただし,バリ)は中心点から無限要素の辺12上の任意点まで
の距離(筒井,1999a,b)である.
要素形状関数が満たすべき重要な条件は,内外の2領域
Q1,Q2を接続する仮想境界「,上,すなわち,§=-1において
内外の波の位相差がゼロとなることである.また,節点112
が半径、の円弧上にあり,内部領域の要素サイズが十分小
さいときには,7(〃)=const.=70と近似することができる.
したがって,これら2条件を考慮すると,
A=exp(iAro)(20)
となり,形状関数ⅣIは次式で与えられる.
恥-;(\-4(\)1\.,(-町(\「ボ)
(21.1)
同様に,補正係数Aにより形状関数の位相を合わせると,他
w-2w-2
Lm上価
些Zu2
l-一一
M〃
’
の形状関数Ⅳi(i=2,3,4)は次式となる.
(16.2)
試行関数,hは4節点での形状関数M(i=1,2,3,4)を用い
て次式で定められる.
。カー,,Ⅳ,+,2Ⅳ2+ハⅣ3+仇M
(17)
ただし,無限遠ではハーハ=oとなることが考慮されている.
形状関数Ⅳ'は,式(15)を参考にして,次式で与えられる.
Ⅱrll
一|’一
44’’
桿忰町一:¥
俘伴割引
凶2凶2凶Zu2
I111
一一
1-31-38-38-3
一一一一一一一一
MMM帆
(18.1)
M-;(\-鍬(\「}苧鐵,(-灘州鶚)
(21.2)
咋響(\-(\「)旱…(川|\「淵)
(21.3)
伽`薑響(\-(\)1苧…(,"佇「鶚)
(21.4)
(182)
(18.3)
3次の無限要素は,無限写像を表す式(16)および要素形状
関数(21)で構成される.無限写像は前報(筒井,1999a)の
式と同じであるから,ここでも変数変換のヤコピヤンは正
規系のりに無関係となる.その結果,要素積分を解析的に実
行することができ,要素行列はF,巳snelの正弦および余弦積
分により陽に表示される.
なお,付録には,3次の無限要素を用いて外部領域を離
(18.4)
散化したときの要素行列が示されている.
なお,2次元の場合においても,減衰特性を表す式(11)が成
立する.したがって,無限写像を表す式(16)の下で,中心点
から出る放射線上では,[email protected]
4.数値計算例および誤差解析
以上は物理変数が単調に減衰する場合の議論であるが,
波動問題においては補正が必要である.時間項をeXp(iのり
(r:時間)と仮定すると,散乱波はexp(-i(Ar-の『)}に比例
する振動特性をもっている.したがって,形状関数(18)は振
動特性を表すexp(-iADなる項を含まねばならない.いま,A
を未知係数とし,式(16)により定められる中心点からの距
素サイズおよび無限要素の次数との関係について検証す
る.用いる支配方程式は,式(1)において非線形項Q=oと
p-mおよびp-3mのモードが含まれている.
離rを用いると,形状関数Ⅳ1は次式で表すことができる.
ここでは,一定水深域での円柱による波の散乱を例に採
り,無限要素を用いた有限要素解析における計算精度と要
して得られる緩勾配方程式である.
数値モデルの諸元は,円柱の半径α=5.0m,仮想境界円「I
の半径「。=10.0m,水深ハー5.0mである.図毛は線形三角形
要素より成る有限要素網を示す.内部領域における総節点
数:404,要素数:687,要素の平均的な辺長:【=07mである.
琉球大学工学部紀要第59号,2000年
43
表-1数値実験諸元
RunT(SCC)L(、)M〃L
iiZy〉
5.143
31.42
□■●●
20.95
0050
より出る2放射線により構成される.前報の線形無限要素
と同様に,物理変数は仮想境界の辺上では線形変化すると
仮定され,辺上に中間節点が存在しないので,無限要素の
生成は内部要素の形状には左右されない.
このモデルに対して,表-1に示す4種類の波が入射する
場合の円柱の周りの波高分布について検討する.ただし,
T:波の周期,L;進行波の波長,k:波数である.
なお,円柱の周りの波動解は,入射波高を1とすると,次
式(田中,1956)で表される.
。□
‘r=ノ。('b)+Z豊1Wb(&…"8
-鵜M,)-2裏'鶚肌…,(22)
ただし,J感:ベッセル関数,ノッ、:第2種ハンケル関数,(脇の:
極座標,('):主変数に関する微分を表す.また,円柱表面
は完全反射壁と仮定する.
水面変動量は複素変数で表されるから,その絶対値を
採ると,任意点における入射波高に対する波高比が得ら
れる.以下では,前報および本報で提案した無限要素に基
づく有限要素法を,補間多項式の次数から,それぞれ1次
および3次のIFEMと呼ぶ.
図-7は円柱の周りの波高分布を示す.上半円が式(22)で
1/22.4
1ノ29.9
1/44.9
数値解析結果である.波は紙面の左側より入射している.
図中の実線は入射波高に対する波高比KZI,破線はx<1
となる等波高線を示し,数値は波高比Kの値を表す.等波
高線の間隔は全て0.2である.
Runlでの円柱前方のK<1となる領域における等波高
線には,1次のrEMの結果と比較して,3次のrEMに
よる改善が認められ,HFEMによる結果とほぼ同程度の
計算精度となっている.しかし,理論において円柱前方に
存在するK=0.4の等波高線がいずれの計算結果において
も見られない.この原因として,前方領域では部分重複波
が発生し,水面振動の間隔が進行波に比べて短くなって
いることが考えられる.したがって,円柱の前方領域にお
ける要素サイズは,この重複波を再現するにはやや過大
であると判断される.また,円柱背後におけるK=0.6の等
波高線に若干の差異が認められ,円柱後方においても推
算精度が劣ることが判る.
Run2における1次のIFEMよる結果では,円柱の背後
に生じるK=0.8の等波高線が再現されていない.しかし,
3次の配EMによる結果には,形状に差異が認められるも
のの,この等波高線が再現されている.円柱後方では明瞭
な重複波は発生しないので,IFEMにおけるこの差異は,
それぞれの計算に用いられた無限要素の次数に依存して
いるものと考えられる.3次の配EMによる結果には局所
的には理論とのわずかな差異が認められるが,全体とし
伽
IIllllI
I
IFEM(lst.)
3.852
1m5.0
IFEM,ハイブリッド型のH1屯M(Chen・Mei,1975)による
外部領域での無限要素は,仮想境界円r,上の辺と中心点
I
15.71
与えられる理論解,下半円は左から順に1次および3次の
図-6円柱の周りの有限要素網(半領域)
IlE
10.48
3.232
321I
1234
蕊ff淫
2.598
'IMllli
1FBM(3m.)
(a)Runl:い=Aα=3,k「o=6
図-7円柱の周りの波高分布
HFEM
44
筒井:沿岸開領域における非線形波動解析のための新しい無限要素(第2報)-高次無限要素と誤差評価一
-ごデー--------
- ̄1.2五
へ1.0-
jq8/--14
=クノ■
-~
1.21.o-
0.8-
-
呵艤6
ノ
、
IFEMUst.)
シ
、JLIミニ二一二〕
。、
ア}./
l<、、、
「』
、画、.、
『
へ、へ.へ/
町
1FEM(3㎡.)
HFEM
(b)Run2:M=Aα=2,Aro=4
nW
IFEMUst.)
1FEM(3I。.)
HFEM
(c)Run3:M=kα=1.5,lbb=3
11屯M(lst.)
IFEM(3m.)
H1毛M
(。)Run4:M=ハロ=1,lcro=2
図-7円柱の周りの波高分布(続き)
て両者は良好な一致を示している.この場合の平均的な
要素サイズは進行波の波長の1/22程度であり,要素サイ
ズの設定に対する1つの指標を示唆している.
次に,Run3においては,Ru、1,2と比べて要素サイズ
が小さく,IFEMおよびH1厄Mによる計算結果と理論との
一致度は良好である.ただし,円柱背後での等高線に極
わずかな差異がある.この誤差は3次のIFEMでは1次の
配EMより改善されている.
Run4においては相対的に要素網が密であり,いずれの
結果も理論結果を良く再現しており,1次のIFEMよる推
算でも十分であることが判る.
ここで,波高分布の詳細をみるため,円柱沿いおよび外
琉球大学工学部紀要第59号,2000年
方の仮想境界円上での水面変動量の実部および虚部の比
較を図-8に例示する.縦軸は入射波高に対する相対値K,
横軸は円柱の中心からの角座標である.実線は理論解(22)
を示す.○および+印は,それぞれ3次のIFEMおよびハ
イブリッド型のH1屯Mを適用して推算した結果である.
45
Runlにおいては,円柱沿いおよび仮想境界円上のいず
れにおいても,推算値と理論値には位相のずれが生じて
いる.この傾向は円柱沿いにおいて顕著である.
Run2に対する計算結果は,e〉160゜において理論解と比
べて最大3-4%程度の過小評価となっている.この影響
2
2
j;翼Wi亜i二
、。。。・U。。。。IOU
Onthecylinder:
kr=3
■
曰
■
ロ
I
ー
0121011
0
K.
K
1
1
■
q
■
■
2’
■
■
,…
1
■
●
■
 ̄
01
麹
1
r1
L」
1m ag
kh=ka=2.kFb=4
-Thcoreticm
oFiniteandinfiniteelemmn0& (31.).
---------+Finiteelementandextenorscnes-
「フミTTH百5575Fで55両7「面三万1
cal
・1m
:::誠:(.、夢,
 ̄
■
■
■
B
●
 ̄
 ̄
■
0
■
■
I.、.0Ⅱ....I..
■■ロ●
■
ロロ□01....I...△、_
50100150
50100150
0(。e獣CCS)
0(dcgmees)
(a)Runl:lUI=kα=3,A「o=6
■
■p■
I●・・●ID・・0100
1
K
K
0
■④〆巴〆 ̄
.kh=ka=15,kⅡb=3
-1 --ThemPqi唾!
.oFiniteandinfiniteclcmems(3Ⅱ。.)
2101
2
(b)Run2:」M=Acz=2,k「o=4
.●・・ID・・ロ0DDODIUU
kh=ka=1,kID=2
-T11eonpqi画I
oFiniteandinfmiteclcmems(31..)
+FimteFlmncntandelq麺ors亜es
●
巳
+Finitcerm便nlzmdextenorSC正s。
●
■■
■■句■■■
■
『
1
1
■
l■■
■
■
÷’
Ima
:1IjLl
kr=3
q
ロ■ロID■QQ
Uoooo
I。.
50100150
0(deglCes)
(c)Run3:ルノt=kα=1.5,1570=3
1,面、雨、=N::::塵麓:
Attheouterboundary:
AtlheouterboundaIy:
I
-1
国...、繧緯霞。
浸圖..§
0
0
kr=2
 ̄■
Q
0
■■■
DnQI■0ODlQOQQIロ
50100150
0(degl窪s)
(。)Run4:ハル=Aα=1,Aro=2
図-8円柱沿いおよび仮想境界円上での水面変動量の実部および虚部
■
46
筒井:沿岸開領域における非線形波動解析のための新しい無限要素(第2報)-高次無限要素と誤差評価一
が図-7(b)の円柱背後の等波高線のわずかな差異として現
どの程度の影響を及ぼすかについて調べた結果は,図-9
れている.しかし,本推算結果は十分な精度を有し,HFEM
および図-10に示す通りである.縦軸は入射波高に対する
Run3においても,円柱背後における波高比をやや過小
波高比の誤差,横軸は角座標である.図中の△,◇,およ
び○印はそれぞれ1次,3次のIFEMおよびハイブリッド
による結果との差異は1%以下である.
評価する傾向が認められる.その結果が図-7(c)での波高
比X=0.8の等波高線の差異となっている.Run4において
は実部および虚部ともに十分な推算精度を持っている.
さらに,これら水面変動量の実部および虚部の誤差が,
円柱沿いおよび仮想境界円上での波高分布の誤差として
Runlにおいては誤差の最大値は約0.04に達している.
誤差は円柱沿いではやや規則的に,仮想境界円上におい
ては不規則に変動している.ここでの特徴として,3次の
11毛MとHFEMの誤差の値および変動特性は,ともに,全
U・。。。U、・ロ。I・0
0
2
0
4
字一
1
・口・・0口
U、
モメヨ
exにriorseries
0口■■
0
I....、....1.
50100150
0(degICes)
0(degrees)
(a)Runl:ハハーハα=3,Aro=6
印側加
(a)Runl:ハル=kα=3,Arh=6
60
.。・・U・・・ロ0ロ●DuIuu
-o-Finiteandexteriorse正巳
Finitcandinfinitceにment§
●
40
0
』。E四
1
0
2
』
釦が
-60xl
●、●IロロDUU
Finiteandinfiniteelcments
1sI.,+3rd.
A卜、
l
Ⅱ.。..’....l
50100150
x
0
6
3
0
■■●■
0
山。■四
0
/’
-GFiniteandexにriorsenes
Finitcandinfinit⑧elemenUS
一一1sし;-←3J..
●
0000m㈹3
--0
642
印扣、
凶。■凶
-60xl
4[
型のHFEMを適用して推算した場合の誤差を示す.
-Q,`戸`属も日出一。ゴ杉8房e旨G影§ミ9-シ・P`肉
庁
-o-Finitcandexl璋riorsenes
●
§o
-20
Finitcandinfinitcclemen15
一一lst.;=←3【..
 ̄
0。ロ.I....、----1-
0
50100150
百歩1sL;-3'。
20
■
‐40
-60xlO3
0
50100150
0(degl巴Cs)
0(degI己es)
(b)Run2:」ヒル=kα=2,klb=4
印側、
-o-Finiteandextcnorsenes
FinitcandinfinitcclemenlB
-←1sし;-←31..
40
全一宮田一へへ_
1
如如が
-60xl
60
u・ ̄U□。・・U・c
0
凶。■四
4[
(b)Run2:ICA=Aα=2,klb=4
今鍼E2E8ヨ8邑邑へ-1邑亀目
「
l■■
 ̄
-o-Finiにandexにmorsenes
Finiteandinfiniteelements
-ヶlsL8→-3[。.
、...I...。Ⅱ.。--ロ
050100150
20
■
§o
-20
-40
-6OxlO-3
8(degr℃es)
050100150
0(degIEes)
(c)Run3:jch=kα=1.5,kro=3
(c)Run3:い=Aα=1.5,krb=3
』。届凶
X
l
印
⑪㈹、0知釦が
SE四
0
50100150
0(deg妃Cs)
-60xl
0
50100150
0(degrees)
(。)Run4:M=Aα=1,km=2
(。)Run4;ハル=Aα=1,k「o=2
図-9円柱沿いの波高比の誤差
図-10仮想境界円上での波高比の誤差
琉球大学工学部紀要第59号,2000年
偏角に渡りほぼ同じであることが挙げられる.
Run2における誤差は高々0.02の程度である.円柱沿い
47
について述べる.図-6に示すように領域は同心円で分割
されているので,各同心円上での誤差のRMSと解析対象
の誤差は3推算法ともにほぼ同じで,その変動はやや規
則的である.一方,仮想境界円上においては,e=120・前後
および0>160゜の領域での3次のIFEMによる誤差の改善
が顕著である.その結果,図-7(b)に示されているように,
1次のIFEMによる結果には見られない円柱背後の波高
領域の分割数〃=L〃(表-1参照)との関係を両対数紙上に
描くと,図-11が得られる.横軸は進行波の波長と要素サ
イズとの比であり,領域全体の平均的な分割数の指標を意
味する.図中の記号で,例えばR3は仮想境界円より第3番
目の同心円を意味する.その他の記号は図-9,1oと同じで
現れている.
ある.なお,第5同心円(R5)上ではR4上とほぼ同じ結果
が得られたので,R5上の結果は省略されている.
比K=0.8の等波高線が,3次のIFEMによる推算結果には
Run3の結果はRun2の推算精度が向上したもので,3次
のIFEMによる誤差が1次の誤差より小さく,HFEMによる
誤差とほぼ同じであることが判る.Run4においては,3推
算法における誤差に顕著な差異は認められない.
最後に,総合的な誤差評価の1つとして,3手法で推算
された入射波高に対する波高比の誤差の自乗平均(RMS)
回
Run3への分割数の増加の効果が少ないことが判る.
外方より第1,2,…,5同心円,円柱沿いと中心に近づくに
=GFiniにeIemcntand
exleriorserics
FinitcandinfhIitcclcmenls
一一Finit cclememand
cxtenorscries
FinitcandinfiniIce IcmcnIs
ヨトlSL
ヨトlSL;-←3rd.
明已屋餡
勺
3次のIFEMによる誤差はHFEMのものより小さく,境界条
件の処理が優れていることを示す.また,1次および3次の
IFEMによる誤差のRMSの変化は相似的であり,Run2から
43
43
43
Attheoutcrboundary
仮想境界円上の誤差についてみると,Ru、2,3に対する
=かIst.;-←3m.
Z
1
0
(U
9876543
、
可
1.5
2
1
1.5
仏
Z
1
23456789
0001
,=Lノ1100
234567810
1o
□
9876543
10
1.5
1
Z
-O-Fmitcclcmcntand
extcrioTscncs
1
0.001
1
心
(U
7
100
43
3
可
(U
4
,=L/1
2
四三厘
5
ヨト1s【.;-←3『。.
9876543
6
、)
7
23456781
(c)第2同心円上
-o-Finitcc1cmenIalld
Cxtcriorscrics
FiniteandinfinitecIemcnts
1
8
可一
I
0
1
9
回
めこ盛
⑫乏塵
1
0.01
43
43
コト1sし;-°-3「。.
10
(b)第1同心円上
-o-FiniにeIemcntand
cxtcriorsenes
Finileandinfini【cclcmeD【s
2
0.001
n=Lノ1100
(a)仮想境界円上
回
コト15t.;-←3m.
001
2
2
】0
-←FiniteeIementand
cxteriorsenes
Finileandinfini【cc1ementR
I
I
0.001
回
FinitcandinfiniteeIemenls
1
23456789
,=LノllOO
(。)第3同心円上
0.001
10
÷ISI.;-←3rd.
234567H9
n=L/1100
0.001
10
に)第4同心円上
図-11同心円上における波高比の誤差のRMSと分割数との関係
23456781
,=L/1100
(0円柱沿い
48
筒井:沿岸開領域における非線形波動解析のための新しい無限要素(第2報)一高次無限要素と誤差評価一
つれて,3次のIFEMによる誤差のRMSは,分割数の増大
とともに単調に減少する傾向を示す.また,Run3に対す
る1次のⅡ垣Mによる推算誤差は,仮想境界円上と同様に,
分割数の増加にもかかわらずRun2と比べてわずかに減少
する程度で,この傾向は円柱近くで特に顕著である.
内部領域で線形要素が用いられている場合には,要素内
の近似の最大誤差はo(42)(△:最大要素サイズ)を満たす
こではM,Aroがともに変化している.したがって,図-11に
基づく収束性に対する直接的な議論はできない.しかし,
IFEMおよびHFEMによる誤差のRMSの勾配は-1.5-2で
あり,上記の誤差特性の一端を示している.
IPEMを用いた場合のRun2からRun3における誤差の変
化に見られるように,与えられた地形および波浪条件の下
で良好な波高推算精度を得るためには,仮想境界円,すなわ
ち,“の設定,無限要素の次数および最適な要素サイズの
選択が重要であることが判る.
無限要素の次数の増加の効果により,当然ながら,3次の
IFEMが1次の11屯Mより優れている.Run2に対する3次
のIFEMによる結果をHFEMの結果と比較すると,3次の
IFEMは仮想境界円付近ではHFEMより優れ,円柱付近では
やや劣る.H1垣Mにおいては散乱振幅を十分満たすように
仮想境界円上での境界解を定めているが,3次の11屯Mに
おいてはO(「-5ノユ)以上の項は無視されている.このような
仮想境界円上での境界条件の処理方法の差異が,両者の
誤差のRMSの差異として現れていると考えられる.
図-12は,解析領域内の全節点での誤差のRMSと分割数
との関係を示す.基本特`性は図-11と同じであるが,3次の
IFEMはRun2においてHFEMと同等の推算精度を達成し,
Run3においてもHFEMの結果に近い.領域全体の誤差の
程度から判断すると,3次の11屯Mによる解析は,図-11
にも示されているように,局所的にはHFEMより若干劣る
と考えられる.しかし,新しい無限要素はHFEMにおける
ある.さらにiRun2,3における3次のIFEMによる推算
波高比の誤差のRMSは,入射波高の高々1%程度である.
したがって,新しい無限要素を用いた有限要素解析に対
する1指標として,領域を進行波の波長のl/20-1/25の
要素サイズで分割し,3次のIFEMを採用すれば,実用的
には十分な精度で波高推算が可能である.
5.結
牽車
と考えられ,推算誤差のRMSはO("-2)(〃:領域の分割数)
で与えられる.分割数"による解の収束性を議論するには,
蛾およびkrOを一定として〃を変化させるべきであるが,こ
取り扱いの煩雑さなどを無くすために提案きれたもので
本研究では,波動場解析のために前報で提案された新し
い無限要素による解析精度の向上を図ろため,3次の無
限要素を誘導し,その推算誤差の特性を数値計算例によ
り検証した.これらの結果は次のように要約される.
(1)散乱波の振幅は散乱源から動径方向に減衰するので,
正規(鋤)系において,5方向にはβ-1Ⅲおよびp-3111のモードが
再現可能であり,刀方向には物理変数の線形変化を仮定した
特殊な無限要素が提案された.
(2)しかし,無限要素と解析領域内の要素との結合が簡
単で,かつ内部要素の形状に左右されないこと,および要
素積分が陽に表示されるという可積分性などの基本的な
特性は保たれている.
(3)3次の無限要素に基づく有限要素法は,境界解を用
いるHFEMと比較すると推算精度は極わずかに劣るが,
仮想境界円の設定および要素サイズとの整合性を図るこ
とにより,十分実用に供することができる.
参考文献
菊池文雄・岡部政之(1986〕:有限要素システム入門,日科技連,東
京,l91pp
田中清(】956):円形島による波の回折,土木学会,第3回海岸工
学講演会講演集,Pp33-35・
筒井茂明・大木洋典(1998):スロープおよびステップ型リーフ上で
の波の非線形挙動,海岸工学論文集,Vol45,JSCE,pP4l-45、
筒井茂明(l999a):沿岸開領域における非線形波動解析のための新
しい無限要素,琉球大学工学部紀要,Vol58,ppl7-27・
筒井茂明(1999b):沿岸開領域における波動場解析のための新しい
無限要素,海岸工学論文集,VOL46,JSCE,pp81-85、
望乞函
Bettes,P、CEmsonandT、CChjam(1984〕:AnewmappedinImiteelement
fbrexに【iorwaveproblcms、NumcricalMethodsinCoup】cdSysにms,
JohnWileySonsLtd,pp489-504、
Chen,HSandCCMci(1975):Hybrid毛lemcntmethodfbrwaterwaves,
0.01
Pmc・ModelIimgTbchniquesConf(Model1ingl975),VoLLpp,63-8L
Tsutsui,S・andDP・Lewis(1992):WaveheightplCdic[ioninunboundcd
coasta】domainswithbathymetricdiscontmuity,CoastalEnginJapan,
JSCE,V01.34,pp・l45-I58
Tsutsui,SandKZKlmami(1993):jumpconditionofcnergynuxatlhe
lineofbathymctricdiscon(inui[yandwavebr巳akingonthereefnaL
CoastalEng・injapan,JSCE,Vol36,pp・'55-175.
Tsutsui,S,KSuzuyamaandHOhki(1998):ModeIequationsofnonlinear
dispe応ivewavesinshallowwaterandanapp1icationofitssimpIesl
versiontowavcevolutionomlhestep-type泥efCoastalEngJour..JSCE,
0001
0
図-12全域での波高比の誤差のRMSと分割数との関係
VOL40,N0.1,pp、41-60.
Zicnkiewicz,0.C,K・Band0,P、Bettess,C・EmsonaJ1dTCChian(1985):
Mappedinfiniteelemen[sfbrextcriorwavcproblems・I、(.』our、fbr
Nemer・Me[hodsinEng.,V01.21,pp」229-1251.
琉球大学工学部紀要第59号,2000年
49
付録.3次の無限要素による要素行列
外部領域Q2での水深を一定と仮定し,式(7)の左辺第2項を無限要素により離散化すると,要素行列{バー)の成分,
要素行列{バー)の成分は次式で
表される.
K;=・關lL;wwIv'“-"wcLM“
(A1)
ここで,式(16)および(21)で構成される3次の無限要素を用いると,要素行列の成分K;は以下の諸式で与えられる.
1111
鋤I111
61111
十十
蛎凡卯8s
1鯨鰔脇蝸
+Ⅱ江兜叫
似十十
++2Z
山山叶伜
蛾恥恥輸納
》“恥叱恥恥11
ら一一
牌伽血澁悴柳叩
“いいトト乢皿
・村句&
C||
rGq33刀⑭
一“&
くjj
“22℃℃いく
j11.’2’’2
I加
’|++
句、“++
83’’
+||
凸1J
州刎輌いい肌“
+++十
++十+
6釧狸調弧66
138
十++
十++
34お配
川碗川、、輸納
山山I1
JJ1
西羽44
十+++
一|’|
I1s8
恥いい恥
34咽岻
22222S3
838
822
一’一
、加川
+++
3SS
111トト-1
22二一一
一一一一一一
m鯉⑭q
服一型雛一m識|麺陞巫澱一m罪一狐配一狸
一一一一一一
⑭q四
yXy
ddd
22222.d
e2c202c2c2
利刺珂呵呵弧一弘肌一町
側いい帆帆川一鉱汎一可
ⅢⅢⅢⅢⅢⅨⅢ
S7412伽-4S84J4燗+16s9416m
'L鰐岬擶ト…川M…、…卜……)I
1111
仰馴岬興長命トル「咄恥ⅢM……トル鶚…wl
岬M腿w端トル川…箒w……'……w)
仰…論|M…川川…■M…'…w艫wl
(A2.1)
(A2.2)
(A2.3)
(A2.4)
(A2.5)
(A2.6)
(A2.7)
(A2.8)
(A2,)
(A2.10)
(A2.11)
50
筒井:沿岸開領域における非線形波動解析のための新しい無限要素(第2報)-高次無限要素と誤差評価一
')B瀞峠売lM…に……'…州ハI
(A2.12〕
11k{等等`・奉読|馴化……,…鱸!…-iM胸l
(A2.13)
lL鶚岬読lM………壱川鈩'Mol
(A2.14)
l]H鰐岬針,……鑿ww…肋I
(A2.15)
IL鶚岬読|`州・伽…ハ…き(…川蝉l
(A2.16)
lilH鰐岬一読卜…-…州FM…(肌…鵬11
鵬岬擶トル禧………トル…蝿I
(A2.17)
(A2.18)
lln鰐岻-諾|`……川佶Mow')
(A2.19)
lL鶚`峠諾ト…ハ.州M……l
(A2.20)
Lvf`・譽薑Ⅲ;'w;…山帆罎紳阿:Q
(A2.21)
lLjv:`・興薑几;lvi…liL:1W・罎豐岬l雲。
(A2.22)
仰帆藝薑岬…z山…ユル仙贄薑一洲豹,「塾c,
(A2.23)
ただし,△=きい!(y1-y2)+難!(L-yI1)+麹(汁,,))は補要素。;の面積であり,ムーュー妬,,△y=y2-y1および
6,=("[+DAy,
。,=(、+l)△x,
62=(3"j+1)△y,
。Z=(3川十1)△x
α】=2(ルー)'0)-6,, α2=2(y1-yO)+6',α3=2(y2-y(1)-62, a4=2(yl-yO)+b2
Cl=2(x2-xo)-.,, c2=2(×1-和)+d,,c3=2(x2-x(】)-.2, c4=2(xl-xo)+d2
と置くと,式(A2)中の係数Sij,Ci,およびDiO,ノー1,2,…)は以下の諸式で与えられる.
I
(A3)
琉球大学工学部紀要第59号,2000年
51
S1,=3。;+6f,S12=3Cf+df,S,3章3。;+6f,S14=3c;+diSl,=3QIq2+b;,S16=3c,c2+df
(A4.1)
S2,=3α,α3+6,62,s22=3c1c3+‘,d2
S23=3CZ2Q4+6162,s24=3C2C4+dld2
S25=3ala4+3a2a3+26162,s26=3cIc4+3c2c3+Zdld2
S27=3qlq4+12a2q3+56162,S28=3clc4+12c2c3+Sdldz
Sm=12α,α4+3a2fz3+56,62,s2,0=12c1c4+3c2c3+5.,.2
(A42)
I
S3]=3。;+b;,S32=3c;十.;,S33=3α:+6;,S34=3.;十.;,S35=3a3a4+b;&`=3c3c4+d;
(A4.3)
S7l=3α1-61,s72=3Cl-.1,373=3α2+61,S74=3c2+dl
(A4.4)
S81=3口,+3α3-6,-62,s82=3Cl+3c3-.1-.2,s83=3α2+3α4+61+62,sM=3c2+3c4+dI+d2
(A45)
S9,=3α3-62,s92=3℃3-.2,s93=3α4+62,S94=3c4+d2
(A46)
Cl=J2m-l-8I4m-l+16I6m-l
C2=ノ2,-,-514,-】+416m-1
C3='2,-,-214,-,+恥"-,
)
D,=(2腕-1)J2m-8(4m-1)川+16(6,-1)ノ6,,,
,2=(2,-1)'2,-5(4,-1)'4m+4(6m-1)J“
D3=(2,-1)'2,,,-2(4腕-1)'4m+(6,-1)'6m
(A4.7)
ハ薑11;;+州一α(川-1(:-s!(,。))い`ユ)
)
"三2に対しては,次の漸化式が成立する.
(A4.8)
また,積分ノノ+Aは次式で定義される.
’…lli;;…州仏鋤
ただし,PC=2km;ノー2,,4,,6腕;A=-1,0,1である.
さらに,式(A5)を算定する際に必要な積分J"("zO)は以
下の諸式で与えられる.
'・薑I;;。-峨薑処I:;州-i…〔A`ロ
'蠅臺I;;士…臺蒜請一☆午]い")
なお,積分I,中のSiおよびCjは次式で定義されるFresnelの
正弦および余弦積分である.
s(2)薑lf響`雌薑f-lr平`“
(A7.1)
。②-r竿。〃
(A7.2)