JP3657928B2

JP3657928B2 - Method for analyzing granular material and continuum

Info

Publication number: JP3657928B2
Application number: JP2002191135A
Authority: JP
Inventors: 裕史恒川; 忠彦塩見
Original assignee: Takenaka Corp
Current assignee: Takenaka Corp
Priority date: 2002-06-28
Filing date: 2002-06-28
Publication date: 2005-06-08
Anticipated expiration: 2022-06-28
Also published as: JP2004037114A

Description

【０００１】
【発明の属する技術分野】
本発明は、粒状体及び連続体の解析方法に係り、特に、地盤と地盤中の杭基礎等の構造物とを含む地盤構造物の地震等による破壊を解析するための粒状体及び連続体の解析方法に関する。
【０００２】
【従来の技術】
地盤と地中の構造物、例えば杭基礎などとの相互作用を解析する場合、従来は有限要素法（以下、ＦＥＭという）が使われてきた。しかしながら、ＦＥＭは連続体の解析手法であり、物体が連続していて伸び縮みによって変形する場合には適した解析手法であるが、物体間に剥離やすべりが起こるような大きな破壊が生じる場合には適切な解析とは言えなかった。
【０００３】
一方、個別要素法（以下、ＤＥＭという）は、各々独立な粒々の動きを逐次追いかける解析手法で、前述のような大きな破壊を生じる解析に適している。このＤＥＭを用いて解析を行う場合に、杭などの連続体として扱った方が良い対象に関しても地盤と同様なモデル化をする手法が提案されているが、解析精度が悪く実用的ではない。
【０００４】
そこで、地盤などの非連続体として扱った方が良い部分はＤＥＭで、それ以外の連続体として扱った方が良い部分は連続体の解析手法でモデル化して解析するハイブリッド法が提案されているが、収束計算の繰り返し回数が多大なため計算時間がかかり過ぎ、実用的ではなかった。また、計算時間を短縮するために誤差を許容すると、精度が落ちて非現実的な解を導く原因となっていた。
【０００５】
ＤＥＭとその他の連続体解析手法とのハイブリッド解析法としては、例えば特開平９−２２１７５４号公報に記載されているように、連続体解析手法としてＦＥＭを採用し、ＦＥＭ要素を矩形のＤＥＭ要素として他のＤＥＭ要素との間の接触力を計算し、その力をＦＥＭ解析の境界条件として与えて解析を行い、収束するまでこの手順を繰り返すものがある。この手法では、ＤＥＭ要素からＦＥＭ要素に作用する接触力に基づいてＦＥＭ要素の変形を計算する。
【０００６】
【発明が解決しようとする課題】
しかしながら、特開平９−２２１７５４号公報に記載された方法では、得られたＦＥＭ要素の変形によってＤＥＭ要素との接触力が変化してしまうため、ＦＥＭ解析で求めた変形状態では安定せず、再度ＦＥＭ解析を実施しなければならない、という問題があった。
【０００７】
また、こうして求めたＦＥＭ要素の変形によってもＤＥＭ要素との間の接触力は再度変化してしまうため、かなりの回数の繰り返し計算が必要になり、解を得るのに大変時間がかかる、という問題もあった。
【０００８】
さらに、ＤＥＭでは要素と要素の新たな接触が生じたり、すべりによって接触が外れたりするためにＦＥＭなどに比べて解析が不安定であり、このため更に繰り返しに時間がかかったり、最終的な収束状態が得られないなどの問題が生じる場合があった。
【０００９】
本発明は、上記問題を解決すべく成されたものであり、計算時間を大幅に短縮することができると共に精度良く地盤構造物の挙動を解析することができる粒状体及び連続体の解析方法を提供することを目的とする。
【００１０】
【課題を解決するための手段】
上記課題を解決するため、請求項１記載の発明は、複数の粒状体である個別要素と当該複数の個別要素中に設けられた連続体である構造物との挙動を解析する粒状体及び連続体の解析方法であって、前記複数の個別要素同士の第１の接触力を個別要素法により各々算出する第１ステップと、複数の構造物要素及び当該複数の構造物要素同士を接続する節点要素から構成された前記構造物と前記個別要素との第２の接触力を個別要素法により各々算出する第２ステップと、前記第２の接触力の各々を前記節点要素に集約した第３の接触力を算出する第３ステップと、前記構造物要素及び前記節点要素に前記第２の接触力を及ぼす個別要素をバネとして表した時の前記バネのバネ定数を前記第２の接触力に基づいて算出する第４ステップと、前記第３の接触力と前記バネ定数とに基づいて、前記節点要素と前記バネとが釣り合うと仮定した時の前記節点要素の変位を算出する第５ステップと、前記節点要素と前記バネとが釣り合うと仮定した時の前記節点要素の変位が予め定めた収束条件を満たさない場合には、前記第５ステップで求めた変位のもとで前記収束条件を満たすまで前記第２ステップ〜第５ステップの処理を繰り返し、前記変位が予め定めた収束条件を満たす場合には、前記第５ステップで求めた変位のもとで前記第１ステップ〜第５ステップの処理を繰り返すことを特徴とする。
【００１１】
この発明によれば、複数の粒状体である個別要素と当該複数の個別要素中に設けられた連続体である構造物との挙動を解析する。例えば、複数の粒状体である個別要素で表すことができる地盤内に、地盤内に設けられた構造物、例えば力を加えると撓むような杭等の構造物の挙動を解析する。
【００１２】
第１ステップでは、地盤を構成する複数の個別要素同士の第１の接触力を個別要素法により各々算出する。
【００１３】
また、構造物は連続体であり、複数の構造物要素及び当該複数の構造物要素同士を接続する節点要素によって表すことができる。このため、第２ステップでは、複数の構造物要素及び当該複数の構造物要素同士を接続する節点要素から構成された構造物と個別要素との第２の接触力を個別要素法により各々算出する。これにより、個別要素と、構造物を構成する構造物要素及び節点要素の各々の挙動を解析することができる。
【００１４】
次の第３ステップでは、第２の接触力の各々を節点要素に集約した第３の接触力を算出する。
【００１５】
そして、第４ステップでは、構造物要素及び節点要素に第２の接触力を及ぼす個別要素をバネとして表した時のバネのバネ定数を第２の接触力に基づいて算出する。
【００１６】
このバネ定数は、例えば請求項２に記載したように、節点要素の変位量と第３の接触力の変化量とから、節点要素の変位量に対する第３の接触力の変化率である傾きを求め、これをバネ定数とすることができる。
【００１７】
また、請求項３に記載したように、バネ定数は、構造物要素及び節点要素に接触する個別要素の個別バネ定数を節点要素に集約することにより求めてもよい。
【００１８】
そして、第５ステップでは、第３の接触力とバネ定数とに基づいて、節点要素とバネとが釣り合うと仮定した時の節点要素の変位を算出する。
【００１９】
そして、節点要素とバネとが釣り合うと仮定した時の節点要素の変位が予め定めた収束条件を満たさない場合には、第５ステップで求めた変位のもとで収束条件を満たすまで収束ループである第２ステップ〜第５ステップの処理を繰り返す。収束条件は、バネと節点要素とが釣り合うと判断できる程度に、節点要素の移動量が小さくなった場合である。
【００２０】
一方、変位が予め定めた収束条件を満たす場合には、第５ステップで求めた変位のもとで時間ループである第１ステップ〜第５ステップの処理を繰り返す。
【００２１】
このように、構造物に接触する個別要素をバネで表して解析を行うため、構造物が変形することによって新たに生じる接触力がバネの伸び縮みによって考慮される。従って、収束ループにおける繰り返し計算の回数を従来と比較して大幅に少なくすることができる。
【００２２】
また、バネ定数を各計算ステップにおいてそのときの複数の個別要素の状態に応じてその都度算出するため、精度よく解析を行うことができる。
【００２３】
なお、解析モデルとしては、地盤及び地盤内に設けられた杭等に限らず、複数の個別要素中に連続体が設けられたモデルであればよい。
【００２４】
【発明の実施の形態】
以下、図面を参照して本発明の実施の形態について説明する。図１には、本発明に係る粒状体及び連続体の解析方法を実施するための地盤構造物解析装置１０の概略ブロック図を示した。
【００２５】
地盤構造物解析装置１０は、図１に示すように、第１の接触力算出部１２、第２の接触力算出部１４、接触力集約部１６、変位算出部１８、及び収束判定部２０を含んで構成されている。
【００２６】
第１の接触力算出部１２は、図２に示すように、地盤を構成する個別要素（ＤＥＭ要素）としての複数の地盤要素３０の挙動をＤＥＭにより解析し、地盤要素３０同士の第１の接触力を算出する。
【００２７】
第２の接触力算出部１４は、第１の接触力算出部１２のＤＥＭによる解析結果に基づき、杭３２と接触する地盤要素３０が杭３２に及ぼす第２の接触力を算出する。
【００２８】
なお、本実施の形態では、詳細は後述するが、杭３２は、図２に示すように杭３２を複数に分割したＤＥＭ要素としての円筒形要素（構造物要素）３４及び、この円筒形要素３４を結ぶＤＥＭ要素としての節点要素３６を含む杭要素で構成されているものとして扱う。
【００２９】
接触力集約部１６は、第２の接触力算出部１４で求めた地盤要素３０が杭３２に及ぼす第２の接触力を杭３２の節点要素３６に集約し、第３の接触力を算出する。
【００３０】
変位算出部１８は、接触力集約部１６で算出した第３の接触力を用いて、杭３２の節点要素３６の変位を求める。このとき、地盤要素３０は、図３に示すように、地盤バネ３８で表されるものとして扱う。
【００３１】
収束判定部２０は、変位算出部１８で求めた節点要素３６の変位が予め定めた収束条件を満たす場合には、その変位を第１の接触力算出部１２へ出力する。これにより、次の時間ループに移行し、変位が生じた節点要素３６のもとで第１の接触力算出部１２、第２の接触力算出部１４、接触力集約部１６、変位算出部１８、収束判定部２０による上記と同様の処理が行われる。
【００３２】
一方、節点要素３６の変位が予め定めた収束条件を満たさない場合には、その変位を第２の接触力算出部１４へ出力する。これにより、変位が生じた節点要素３６のもとで第２の接触力算出部１４、接触力集約部１６、変位算出部１８、収束判定部２０による上記と同様の処理が行われる。この処理は、節点要素３６の変位が予め定めた収束条件を満たすまで繰り返される。
【００３３】
次に、本実施の形態の作用として、地盤構造物解析装置１０の各部の動作について、図４に示すフローチャートを参照して説明する。なお、ここでは、図２に示すように、解析領域３７の地盤中に設けられた杭３２の杭頭に強制的に変位を加えた場合における杭３２の各部の挙動を解析する場合を例に説明する。
【００３４】
まず、ステップ１００では、第１の接触力算出部１２により、地盤を構成する複数の地盤要素３０及び地盤内に設けられた杭３２を構成する円筒形要素３４及び節点要素３６の各々の挙動をＤＥＭにより解析し、地盤要素３０同士の第１の接触力を算出する。
【００３５】
ＤＥＭでは、図２に示すように、地盤は離散化された地盤要素３０の集合体で表される。そして、各地盤要素３０は剛体であり、各地盤要素３０間の力の伝達は各地盤要素３０の接点においてのみ行われると仮定して、各要素毎に、次式で示す要素間に働く法線方向と接線方向の力を表す式を用いて時刻歴解析により各要素の挙動が解析される。また、この解析においては、杭要素である円筒形要素３４及び節点要素３６も解析領域３７の境界部として扱われ、固定されているものとして解析が行われる。
【００３６】
【数１】

【００３７】
なお、地盤に関する材料定数、例えば地盤要素間及び地盤要素−杭要素間のバネ定数や粘性係数、摩擦係数、単位体積重量などは予め設定される。また、時間ループの時間間隔等も予め設定される。
【００３８】
そして、第１の接触力算出部１２は、ＤＥＭによる解析結果から、地盤要素３０同士の第１の接触力を算出する。
【００３９】
次に、ステップ１０２において、第２の接触力算出部１４により、杭要素、すなわち円筒形要素３４及び節点要素３６に接する地盤要素３０が、円筒形要素３４及び節点要素３６に及ぼす第２の接触力を算出する。この接触力は、地盤要素−杭要素間のバネ定数や地盤要素及び杭要素の相対変位等に基づいて算出される。
【００４０】
次に、ステップ１０４において、接触力集約部１６により、図５（Ａ）に示すように、第２の接触力算出部１４で算出された円筒形要素３４及び節点要素３６に接する地盤要素３０が円筒形要素３４及び節点要素３６に及ぼす第２の接触力Ｆ₂を、各節点要素３６に集約して第３の接触力Ｆ₃を各々算出する。
【００４１】
第３の接触力Ｆ₃は、例えば、各節点要素３６について予め定めた範囲内の第２の接触力を、第２の接触力が及ぼされる位置から節点要素３６までの距離に応じた比例配分で各節点要素３６に各々配分し、各節点要素３６に配分された各々の第２の接触力を各節点要素３６毎に合計することにより算出することができる。ここで、予め定めた範囲は、例えば算出対象となる節点要素３６に接続された両隣の円筒形要素３４までの範囲とすることができる。
【００４２】
例えば、各々の第２の接触力をＦ_2iとした場合（ｉは予め定めた範囲内の各々の第２の接触力の番号を表す添え字）、各節点要素３６に集約される第３の接触力Ｆ₃は、次式で求めることができる。
【００４３】
【数２】

【００４４】
ここで、ａ_iは算出対象の節点要素３６に接続された円筒形要素３４の、算出対象の節点要素３６と反対側の端部から、その円筒形要素３４に第２の接触力Ｆ_2iが及ぼされる位置までの距離であり、ｌは節点要素３６に接続された円筒形要素３４の長さである。例えば図５（Ｂ）に示すように、算出対象が節点要素３６₂の場合には、ａ_iは円筒形要素３４₂の節点要素３６₂と反対側の端部から第２の接触力Ｆ_2iが及ぼされる位置までの距離又は円筒形要素３４₃の節点要素３６₂と反対側の端部から円筒形要素３４₃に第２の接触力Ｆ_2iが及ぼされる位置までの距離となり、ｌは円筒形要素３４₂又は円筒形要素３４₃の長さとなる。
【００４５】
次に、ステップ１０６において、変位算出部１８により、地盤要素３０を地盤バネ３８として表したときの地盤バネ３８と節点要素３６とが釣り合う釣合点（安定点）を求めることにより、節点要素３６の変位を算出する。
【００４６】
この節点要素３６の変位の算出の詳細について説明する。まず、杭３２の杭頭に強制的に変位が加えられていない状態、すなわち時間ループの第１ステップにおける杭３２の各部の変位の算出について説明する。
【００４７】
杭頭に強制的に変位が加えられていない状態であっても、地盤要素３０のつまり具合の不均衡から、節点要素３６にかかる力の合力はゼロにはならない。従って、初期状態で杭３２が安定する点を見つけだす必要がある。
【００４８】
ここでは、図６に示すように、中間に２つの節点要素３６₁、３６₂を持つ杭３２について考え、この節点要素３６₁、節点要素３６₂に第３の接触力を及ぼす地盤要素を地盤バネ３８₁、３８₂で表す。この節点要素３６₁と地盤バネ３８₁、節点要素３６₂と地盤バネ３８₂が釣り合う安定点は、節点要素３６₁、３６₂を幾らか動かすことによって見つけ出すものとし、節点要素３６₁をδ₁ ⁽¹⁾だけ動かし、節点要素３６₂をδ₂ ⁽¹⁾だけ動かして地盤バネ３８と釣り合ったと仮定すると、次式が成り立つ。
【００４９】
Ｐ₁ ^'(1)−Ｐ₁ ⁽¹⁾＝ｋ₁δ₁ ⁽¹⁾ …（３）
Ｐ₂ ^'(1)−Ｐ₂ ⁽¹⁾＝ｋ₂δ₂ ⁽¹⁾ …（４）
ただし、Ｐ₁ ⁽¹⁾は、移動前に節点要素３６₁に作用している力、Ｐ₁ ^'(1)は、移動後に節点要素３６₁に作用している力である。同様に、Ｐ₂ ⁽¹⁾は、移動前に節点要素３６₂に作用している力、Ｐ₂ ^'(1)は、移動後に節点要素３６₂に作用している力である。ｋ₁、ｋ₂は地盤バネ３８₁、３８₂の地盤バネ定数である。この地盤バネ定数の算出については後述する。この地盤バネ定数は本発明のバネ定数に相当する。
【００５０】
節点要素３６₁、３６₂を移動させた後に、杭３２の剛性と地盤からの力が釣り合うのであるから、杭３２の変位と地盤からの力の関係は、次式で表される。
【００５１】
【数３】

【００５２】
ただし、Ｅは杭３２の弾性係数、Ｉは杭３２の断面２次モーメントである。
【００５３】
上記（３）〜（６）式をまとめると以下のようになる。
【００５４】
【数４】

【００５５】
上記（８）式をＰ₁ ^'(1)は、Ｐ₂ ^'(1)について解くことにより、安定点における地盤からの力を算出することができる。そして、これを上記（５）、（６）式に代入することにより、節点要素３６₁、３６₂をどれくらい移動させればよいかを求めることができる。
【００５６】
しかしながら、こうして決定した点に節点要素３６₁、３６₂を移動させても、ＤＥＭによる要素間の相互作用は非線形であるため、移動させた位置で杭３２が安定することは少ない。これは、実際の地盤バネが非線形であることに対応する。このため、引き続き安定点を見つけるための収束計算を行う必要がある。
【００５７】
ここで、前述した最初の計算と同様に、節点要素３６₁をΔδ₁ ⁽²⁾だけ動かし、節点要素３６₂をΔδ₂ ⁽²⁾だけ動かして地盤バネ３８と釣り合ったと仮定すると、次式が成り立つ。
【００５８】
Ｐ₁ ^'(2)−Ｐ₁ ⁽²⁾＝ｋ₁Δδ₁ ⁽²⁾ …（９）
Ｐ₂ ^'(2)−Ｐ₂ ⁽²⁾＝ｋ₂Δδ₂ ⁽²⁾ …（１０）
ところで、Δδ_i ⁽²⁾とその時点での節点要素３６_iとのたわみδ_i ⁽²⁾とには、以下のような関係がある。
【００５９】
δ_i ⁽²⁾＝δ_i ⁽¹⁾＋Δδ_i ⁽²⁾ …（１１）
また、地盤バネ定数を後述する第１の算出方法で算出した場合は、１回目の繰り返し計算で算出した節点要素３６_iに作用する力と２回目の繰り返し計算で算出した節点要素３６_iに作用する力との差と、１回目に移動させた節点要素３６_iとの関係は以下のようになる。
【００６０】
Ｐ_i ⁽²⁾−Ｐ_i ⁽¹⁾＝ｋ_iδ_i ⁽¹⁾ …（１２）
上記（１１）、（１２）式を上記（９）、（１０）式に代入すると、以下のような関係が得られる。
【００６１】
【数５】

【００６２】
【数６】

【００６３】
これらをまとめると以下のようになる。
【００６４】
【数７】

【００６５】
上記（１５）式をＰ₁ ^'(2)、Ｐ₂ ^'(2)について解くことにより、安定点での地盤からの力を算出することができる。そして、これを次式に代入することにより、節点要素３６₁、３６₂をどこへ移動させればよいかを求めることができる。
【００６６】
【数８】

【００６７】
なお、地盤バネ定数を後述する第２の算出方法で算出した場合には、上記（１２）式の関係が成り立たないため、上記（１５）式は以下のような形になる。
【００６８】
【数９】

【００６９】
以下同様に収束計算を繰り返し、以下のような条件を満たす場合に安定点が見つかったと判断する。
【００７０】
【数１０】

【００７１】
次に、杭３２の杭頭に強制的に変位を加えた状態、すなわち時間ループの第ｎステップ（ｎは自然数）における杭３２の各部の変位の算出について説明する。
【００７２】
ここでは、図７に示すように、中間に１つの節点要素３６を持つ杭３２について考え、この節点要素３６に第３の接触力を及ぼす地盤要素を地盤バネ３８で表す。
【００７３】
第１ステップと同様に考えると、節点要素３６と地盤バネ３８が釣り合ったと仮定すると、次式が成り立つ。
【００７４】
Ｐ^'(1)−Ｐ⁽¹⁾＝ｋΔδ_t ⁽¹⁾ …（２０）
ここで、Δδ_t ⁽¹⁾は節点要素３６の移動量であり、このΔδ_t ⁽¹⁾と杭３２のたわみとの関係は、図７から以下のようになる。
【００７５】

これを上記（２０）式に代入すると、次式のような関係が得られる。
【００７６】

杭３２を移動させると、杭３２の剛性と地盤からの力が釣り合うため、こうした節点要素３６が複数ある場合、その時の地盤からの力と杭３２のたわみとの関係は力のベクトルＰ^'(1)とたわみのベクトルδ_t ⁽¹⁾を用いて次式で表される。
【００７７】
【数１１】

【００７８】
ここで、［Ｈ₀］は節点要素３６に働く力とたわみとの関係を表す行列であり、例えば図６に示すように２つの節点要素３６₁，３６₂を持つ杭３２の場合には、次式で表される。
【００７９】
【数１２】

【００８０】
【数１３】

【００８１】
上記（２５）式をＰ^'(1)に関して解くことにより、安定点での地盤からの力を算出することができる。そして、求めたＰ^'(1)を上記（２３）式に代入することにより、節点要素３６をどこへ移動させればよいかを求めることができる。
【００８２】
第ｎステップの収束ループにおける２回目の繰り返し計算以降は、第１ステップの収束ループにおける２回目の繰り返し計算以降と同様のアルゴリズムで安定点を求めることができる。
【００８３】
なお、杭３２に強制的に変位を与えない第ｎステップの処理に関しては、ｙ_t＝ｙ_t-Δ_tとして上記のアルゴリズムと同様の処理により安定点を求めることができる。
【００８４】
次に、地盤バネ３８の地盤バネ定数の算出について説明する。地盤バネ定数の算出方法としては、収束ループによる繰り返し計算において、図８に示すように杭が変形することにより節点要素３６が移動した距離と、節点要素３６が移動したことによる地盤要素３０からの力の変化に基づいて求める第１の算出方法と、図９に示すように、杭３２に接触している地盤要素３０の個々を表すバネ４０の個別バネ定数を合成して地盤バネ３８の地盤バネ定数を算出する第２の方法とがある。
【００８５】
まず、第１の算出方法について説明する。
【００８６】
時間ループの第１ステップの場合は、収束ループの１回目の計算では杭３２を動かさないため、節点要素３６を移動した距離と、節点要素３６が移動したことによる地盤からの力の変化に基づいて求めることはできない。収束ループによる２回目の繰り返し計算では、節点要素３６₁に対して変位δ₁ ⁽¹⁾を与えたことによって、節点要素３６₁に加わる力は、Ｐ₁ ⁽¹⁾からＰ₁ ⁽²⁾に変化する。従って、フックの法則から、地盤バネ定数ｋ₁を求めることは、図１０に示すように、杭３２の変形の変化、すなわち節点要素３６₁の変位と節点要素３６₁に加わる力の変化とから力の変化率を表す傾きを求めることに等しい。これにより、地盤バネ定数ｋ₁は次式で求めることができる。
【００８７】

この関係は、第１ステップのアルゴリズムで仮定した上記（１２）式に等しい。
【００８８】
第ｎステップの収束ループによる１回目の計算では、ｔ−２Δｔからｔ−Δｔまでの節点要素３６の移動量とその間の地盤から節点要素３６に加わる力の変化量を用いて、次式により地盤バネ定数ｋを求めることができる。
【００８９】
【数１４】

【００９０】
こうして求めた地盤バネ定数は、実際の地盤の状態をそのまま反映しており、地盤バネを比較的良質の非線形バネと考えた場合には、ある時点からある時点までの割線合成と考えることができる。
【００９１】
このように、それぞれの計算ステップにおいて、そのときの節点要素３６の変位と節点要素３６に加わる力の変化から地盤バネ定数を求めるため、より精度よく地盤構造物の状態を解析することができる。
【００９２】
次に、第２の算出方法について説明する。
【００９３】
ＤＥＭでは、要素と要素とは法線方向と接線方向のバネを介して接触力を伝達している。第２の算出方法では、このバネのバネ定数を直接地盤バネ定数に変換する。円筒形要素３４及び節点要素３６には複数の地盤要素３０が接触しており、これらの要素のバネが並列に接続しているとして、これらのバネのバネ定数（個別バネ定数）を合算することにより地盤バネ定数を求める。これは、それぞれの計算ステップにおいて行われる。
【００９４】
なお、各要素の運動及び節点要素の移動は各方向成分毎に計算されるため、以下では個々の要素のバネを各方向成分に分解する方法について説明する。節点要素３６と接触している地盤要素３０との位置関係が図１１ような関係にあったとする。この要素のバネ定数をｋ、分解したｘ方向のバネ定数をＫ_xとする。θ＜＜１とすると、以下の関係が得られる。
【００９５】

また、図１１から、以下の関係となることが判る。
【００９６】
【数１５】

【００９７】
このように、第１の算出方法及び第２の算出方法では、各計算ステップにおいて、そのときの杭３２と地盤要素３０との接触力から地盤バネ定数を求めるため、精度よく地盤構造物の解析を行うことができる。また、解析の結果として、図１２に示すような杭３２の変位と地盤バネ定数との関係を得ることができる。
【００９８】
このようにして、地盤バネ定数を求めて節点要素３６の変位を求めた後は、ステップ１０８において、収束判定部２０により、変位算出部１８で算出した節点要素３６の変位が予め定めた収束条件を満たすか否か、すなわち、上記（１９）式を満たすか否かを判断する。
【００９９】
そして、変位が収束しておらず、収束条件を満たさない場合には、ステップ１０８の判断が否定され、ステップ１０２へ移行し、ステップ１０２〜ステップ１０８の処理を繰り返す。
【０１００】
一方、変位が収束し、収束条件を満たす場合には、ステップ１０８の判断が肯定され、ステップ１１０へ以降する。
【０１０１】
ステップ１１０では、強制的に変位を加えていた杭３２の杭頭の移動が停止したか否かを判断する。そして、杭頭が停止していない場合には、ステップ１１０の判断が否定されてステップ１００へ移行し、時間ループにより次の時間ステップの処理を上記と同様にてして行う。一方、杭頭が停止した場合には、ステップ１１０の判断が肯定され、本ルーチンを終了する。
【０１０２】
このように、本実施の形態では、杭３２に接触する地盤要素３０を地盤バネで表して解析を行うため、杭３２が変形することによって新たに生じる接触力が地盤バネの伸び縮みによって考慮される。従って、収束ループにおける繰り返し計算の回数が従来と比較して大幅に少なくすることができる。
【０１０３】
例えば、ＤＥＭでは1/10,000秒から1/100,000秒という短い時間間隔を単位として時々刻々と要素の動きを追っていく解析手法であるため、１０秒程度の解析を行うだけで１００万ステップに及ぶ解析を行わなければならず、解析には非常に時間がかかっていた。
【０１０４】
また、ＤＥＭと他の解析手法とのハイブリッド解析では、この時間ステップ毎に前述した釣り合い計算をしなければならず、仮に１００回の繰り返し計算が必要だとすると、１００万回の時間ステップの解析をするためには全体では１億回の計算が必要になる。
【０１０５】
これに対し、本発明を利用すれば、この繰り返し回数が数回から多くても１０回以下となるため、計算時間を大幅に短縮することができ、迅速に計算結果を得ることができる。
【０１０６】
また、図１３には、本発明により杭３２の挙動を解析した結果を、図１４には、収束条件を満足しない誤差を含んだ状態で時間ループの処理を進めた場合における解析結果を示した。図１３、図１４共に、横軸は杭３２の杭頭の移動距離を示し、縦軸は杭３２に加わる力を示す。
【０１０７】
図１４から明らかなように、誤差を残したまま時間ループを進めた場合には、杭３２に加わる力が大きく変動しており、全体の解析が安定していないのが判る。一方、図１３から明らかなように、本発明による解析では、杭３２に加わる力の変動が少なく、全体の解析が安定しているのが判る。
【０１０８】
【発明の効果】
以上説明したように、本発明によれば、計算時間を大幅に短縮することができると共に精度良く地盤構造物の挙動を解析することができる、という効果を有する。
【図面の簡単な説明】
【図１】地盤構造物解析装置の概略ブロック図である。
【図２】地盤構造物の解析モデルを表す図である。
【図３】地盤要素及び杭要素のモデルを表す図である。
【図４】地盤構造物解析装置による処理の流れを表すフローチャートである。
【図５】第２の接触力の集約について説明するためのイメージ図である。
【図６】節点要素の変位について説明するための図である。
【図７】節点要素の変位について説明するための図である。
【図８】第１の算出方法における地盤バネ定数の算出について説明するための図である。
【図９】第２の算出方法における地盤バネ定数の算出について説明するための図である。
【図１０】第１の算出方法における地盤バネ定数の算出について説明するための図である。
【図１１】第２の算出方法における地盤バネ定数の算出について説明するための図である。
【図１２】本発明の解析による杭頭の変位と地盤バネ定数との関係を表すグラフである。
【図１３】本発明の解析による杭頭の変位と杭に加わる力との関係を表すグラフである。
【図１４】従来の解析による杭頭の変位と杭に加わる力との関係を表すグラフである。
【符号の説明】
１０地盤構造物解析装置
１２接触力算出部
１４接触力算出部
１６接触力集約部
１８変位算出部
２０収束判定部
３０地盤要素
３１地盤要素
３２杭
３４円筒形要素
３６節点要素
３７解析領域
３８地盤バネ
４０バネ[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a method for analyzing a granular material and a continuum, and in particular, a granular material and a continuum for analyzing a ground structure including a ground and a structure such as a pile foundation in the ground due to an earthquake or the like. It relates to the analysis method.
[0002]
[Prior art]
Conventionally, the finite element method (hereinafter referred to as FEM) has been used to analyze the interaction between the ground and an underground structure such as a pile foundation. However, FEM is a continuum analysis method, and is a suitable analysis method when an object is continuous and deforms due to expansion and contraction. However, when a large fracture occurs that causes separation or slippage between objects. Was not an appropriate analysis.
[0003]
On the other hand, the individual element method (hereinafter referred to as DEM) is an analysis method that sequentially tracks the movement of each independent grain, and is suitable for the analysis that causes the above-described large destruction. In the case of performing analysis using this DEM, a technique for modeling the same object as the ground is proposed for an object that should be handled as a continuum such as a pile, but the analysis accuracy is poor and impractical.
[0004]
Therefore, a hybrid method has been proposed in which the portion that should be treated as a discontinuous body such as the ground is DEM, and the portion that should be treated as a continuum other than that is modeled and analyzed by a continuum analysis method. However, since the number of repetitions of the convergence calculation is large, it takes too much calculation time and is not practical. In addition, if an error is allowed in order to shorten the calculation time, the accuracy is lowered, leading to an unrealistic solution.
[0005]
As a hybrid analysis method of DEM and other continuum analysis methods, for example, as described in JP-A-9-221754, FEM is adopted as a continuum analysis method, and FEM elements are used as rectangular DEM elements. In some cases, the contact force between other DEM elements is calculated, the force is given as a boundary condition for FEM analysis, analysis is performed, and this procedure is repeated until convergence. In this method, the deformation of the FEM element is calculated based on the contact force acting on the FEM element from the DEM element.
[0006]
[Problems to be solved by the invention]
However, in the method described in Japanese Patent Application Laid-Open No. 9-221754, the contact force with the DEM element changes due to the deformation of the obtained FEM element. There was a problem that FEM analysis had to be performed.
[0007]
In addition, since the contact force with the DEM element changes again due to the deformation of the FEM element thus obtained, a considerable number of repeated calculations are required, and it takes a very long time to obtain a solution. There was also.
[0008]
In addition, in DEM, new contact between elements occurs, or contact is lost due to slipping, so analysis is unstable compared to FEM, etc., so it takes more time to repeat and final convergence There were cases where problems such as inability to obtain the state occurred.
[0009]
The present invention has been made to solve the above problems, and provides a method for analyzing a granular material and a continuum that can greatly reduce the calculation time and can accurately analyze the behavior of a ground structure. The purpose is to provide.
[0010]
[Means for Solving the Problems]
In order to solve the above-mentioned problem, the invention according to claim 1 is directed to a granular body for analyzing the behavior of an individual element that is a plurality of granular bodies and a structure that is a continuous body provided in the plurality of individual elements. A body analysis method, wherein a first step of calculating a first contact force between a plurality of individual elements by an individual element method, a plurality of structure elements, and a node connecting the plurality of structure elements A second step of calculating a second contact force between the structure composed of elements and the individual element by an individual element method, and a third step of integrating each of the second contact forces into the node element. A third step of calculating a contact force, and a spring constant of the spring when the individual element that exerts the second contact force on the structure element and the node element is represented as a spring based on the second contact force A fourth step of calculating Contact force of 3 and on the basis of the spring constant, and the fifth step of calculating the displacement of the node element when assuming said node element and said spring are balanced, theOf the nodal element when it is assumed that the nodal element and the spring are balancedIf the displacement does not satisfy a predetermined convergence condition, the processes in the second to fifth steps are repeated until the convergence condition is satisfied under the displacement obtained in the fifth step, and the displacement is determined in advance. If the convergence condition is satisfied, the processing from the first step to the fifth step is repeated under the displacement obtained in the fifth step.
[0011]
According to this invention, the behavior of the individual elements which are a plurality of granular bodies and the structure which is a continuous body provided in the plurality of individual elements is analyzed. For example, the behavior of a structure provided in the ground, for example, a structure such as a pile that bends when a force is applied, is analyzed in the ground that can be represented by individual elements that are a plurality of granular bodies.
[0012]
In the first step, a first contact force between a plurality of individual elements constituting the ground is calculated by the individual element method.
[0013]
The structure is a continuum, and can be represented by a plurality of structure elements and node elements connecting the plurality of structure elements. For this reason, in the second step, the second contact force between the structure composed of the plurality of structure elements and the node elements connecting the plurality of structure elements and the individual elements is calculated by the individual element method. . Thereby, it is possible to analyze the behaviors of the individual elements and the structure elements and node elements constituting the structure.
[0014]
In the next third step, a third contact force is calculated by integrating each of the second contact forces into the node elements.
[0015]
In the fourth step, the spring constant of the spring when the individual element that exerts the second contact force on the structure element and the node element is represented as a spring is calculated based on the second contact force.
[0016]
For example, as described in claim 2, the spring constant is obtained by calculating a slope, which is a rate of change of the third contact force with respect to the amount of displacement of the node element, from the amount of displacement of the node element and the amount of change of the third contact force. This can be determined as a spring constant.
[0017]
According to a third aspect of the present invention, the spring constant may be obtained by aggregating the individual spring constants of the individual elements contacting the structure element and the nodal element into the nodal elements.
[0018]
In the fifth step, based on the third contact force and the spring constant, the displacement of the nodal element when it is assumed that the nodal element and the spring are balanced is calculated.
[0019]
AndAssuming that the node element and the spring are balanced,When the displacement does not satisfy the predetermined convergence condition, the processes of the second step to the fifth step that are the convergence loop are repeated until the convergence condition is satisfied under the displacement obtained in the fifth step. The convergence condition is when the amount of movement of the node element is small enough to determine that the spring and the node element are balanced.
[0020]
On the other hand, when the displacement satisfies a predetermined convergence condition, the processes of the first step to the fifth step, which are time loops, are repeated under the displacement obtained in the fifth step.
[0021]
In this way, since the individual elements that contact the structure are represented by the spring for analysis, the contact force newly generated by the deformation of the structure is taken into account by the expansion and contraction of the spring. Therefore, the number of repeated calculations in the convergence loop can be significantly reduced compared to the conventional case.
[0022]
Further, since the spring constant is calculated each time according to the state of the plurality of individual elements at each calculation step, the analysis can be performed with high accuracy.
[0023]
The analysis model is not limited to the ground and a pile provided in the ground, and may be a model in which a continuous body is provided in a plurality of individual elements.
[0024]
DETAILED DESCRIPTION OF THE INVENTION
Embodiments of the present invention will be described below with reference to the drawings. In FIG. 1, the schematic block diagram of the ground structure analysis apparatus 10 for enforcing the analysis method of the granular material which concerns on this invention, and a continuous body was shown.
[0025]
As shown in FIG. 1, the ground structure analysis apparatus 10 includes a first contact force calculation unit 12, a second contact force calculation unit 14, a contact force aggregation unit 16, a displacement calculation unit 18, and a convergence determination unit 20. It is configured to include.
[0026]
As shown in FIG. 2, the first contact force calculation unit 12 analyzes the behavior of a plurality of ground elements 30 as individual elements (DEM elements) constituting the ground using a DEM, and first ground elements 30 to each other. Calculate the contact force.
[0027]
The second contact force calculation unit 14 calculates a second contact force exerted on the pile 32 by the ground element 30 in contact with the pile 32 based on the analysis result by the DEM of the first contact force calculation unit 12.
[0028]
Although details will be described later in this embodiment, the pile 32 includes a cylindrical element (structure element) 34 as a DEM element obtained by dividing the pile 32 into a plurality of pieces as shown in FIG. 34 is treated as a pile element including a node element 36 as a DEM element connecting the elements 34.
[0029]
The contact force aggregating unit 16 aggregates the second contact force exerted on the pile 32 by the ground element 30 obtained by the second contact force calculating unit 14 on the node element 36 of the pile 32, and calculates a third contact force. .
[0030]
The displacement calculation unit 18 obtains the displacement of the node element 36 of the pile 32 using the third contact force calculated by the contact force aggregation unit 16. At this time, the ground element 30 is treated as represented by a ground spring 38 as shown in FIG.
[0031]
The convergence determination unit 20 outputs the displacement to the first contact force calculation unit 12 when the displacement of the node element 36 obtained by the displacement calculation unit 18 satisfies a predetermined convergence condition. As a result, the process proceeds to the next time loop, and the first contact force calculation unit 12, the second contact force calculation unit 14, the contact force aggregation unit 16, and the displacement calculation unit 18 under the nodal element 36 where the displacement has occurred. The same processing as described above by the convergence determination unit 20 is performed.
[0032]
On the other hand, when the displacement of the nodal element 36 does not satisfy a predetermined convergence condition, the displacement is output to the second contact force calculation unit 14. As a result, the same processing as described above is performed by the second contact force calculation unit 14, the contact force aggregation unit 16, the displacement calculation unit 18, and the convergence determination unit 20 under the nodal element 36 where the displacement has occurred. This process is repeated until the displacement of the nodal element 36 satisfies a predetermined convergence condition.
[0033]
Next, as an operation of the present embodiment, the operation of each part of the ground structure analyzing apparatus 10 will be described with reference to the flowchart shown in FIG. Here, as shown in FIG. 2, as an example, the behavior of each part of the pile 32 when the displacement is forcibly applied to the pile head of the pile 32 provided in the ground of the analysis region 37 is taken as an example. explain.
[0034]
First, in step 100, the first contact force calculation unit 12 determines the behavior of each of the plurality of ground elements 30 constituting the ground and the cylindrical elements 34 and the nodal elements 36 constituting the piles 32 provided in the ground. The first contact force between the ground elements 30 is calculated by analyzing with DEM.
[0035]
In the DEM, as shown in FIG. 2, the ground is represented by an aggregate of discrete ground elements 30. The local board element 30 is a rigid body, and it is assumed that the transmission of force between the local board elements 30 is performed only at the contact points of the local board elements 30. The behavior of each element is analyzed by time history analysis using expressions representing the forces in the line direction and the tangential direction. Further, in this analysis, the cylindrical element 34 and the nodal element 36 which are pile elements are also handled as boundaries of the analysis region 37 and are analyzed as being fixed.
[0036]
[Expression 1]

[0037]
In addition, the material constant regarding the ground, for example, the spring constant between the ground elements and between the ground elements and the pile elements, the viscosity coefficient, the friction coefficient, the unit volume weight, and the like are set in advance. Also, the time interval of the time loop is set in advance.
[0038]
And the 1st contact force calculation part 12 calculates the 1st contact force of the ground elements 30 from the analysis result by DEM.
[0039]
Next, in step 102, the second contact force calculation unit 14 causes the second contact that the pile element, that is, the ground element 30 in contact with the cylindrical element 34 and the node element 36, exerts on the cylindrical element 34 and the node element 36. Calculate the force. This contact force is calculated based on the spring constant between the ground element and the pile element, the relative displacement of the ground element and the pile element, and the like.
[0040]
Next, in step 104, as shown in FIG. 5 (A), the contact force aggregating unit 16 causes the ground element 30 in contact with the cylindrical element 34 and the nodal element 36 calculated by the second contact force calculating unit 14 to move. Second contact force F acting on the cylindrical element 34 and the nodal element 36₂Are integrated into each node element 36 and the third contact force F_ThreeAre calculated respectively.
[0041]
Third contact force F_ThreeFor example, the second contact force within a predetermined range for each node element 36 is distributed to each node element 36 in proportion to the distance from the position where the second contact force is applied to the node element 36, respectively. It can be calculated by distributing and summing each second contact force allocated to each node element 36 for each node element 36. Here, the predetermined range can be, for example, a range up to the adjacent cylindrical element 34 connected to the node element 36 to be calculated.
[0042]
For example, each second contact force is F_2i(I is a subscript indicating the number of each second contact force within a predetermined range), the third contact force F aggregated in each node element 36_ThreeCan be obtained by the following equation.
[0043]
[Expression 2]

[0044]
Where a_iIs the second contact force F applied to the cylindrical element 34 from the end of the cylindrical element 34 connected to the node element 36 to be calculated on the side opposite to the node element 36 to be calculated._2iIs the distance to the position where l is applied, and l is the length of the cylindrical element 34 connected to the nodal element 36. For example, as shown in FIG.₂In the case of_iIs a cylindrical element 34₂Nodal element 36₂2nd contact force F from the end opposite to the_2iThe distance to the position at which it is exerted or the cylindrical element 34_ThreeNodal element 36₂Cylindrical element 34 from the opposite end_Three2nd contact force F_2iIs the distance to the position where l is applied, l is the cylindrical element 34₂Or cylindrical element 34_ThreeIt becomes the length.
[0045]
Next, in step 106, the displacement calculating unit 18 obtains a balance point (stable point) where the ground spring 38 and the node element 36 are balanced when the ground element 30 is represented as the ground spring 38. Calculate the displacement.
[0046]
Details of the calculation of the displacement of the node element 36 will be described. First, the calculation of the displacement of each part of the pile 32 in the state in which no displacement is forcibly applied to the pile head of the pile 32, that is, the first step of the time loop will be described.
[0047]
Even if the displacement is not forcibly applied to the pile head, the resultant force of the force applied to the node element 36 does not become zero due to the imbalance of the ground element 30, that is, the condition. Therefore, it is necessary to find a point where the pile 32 is stable in the initial state.
[0048]
Here, as shown in FIG. 6, there are two node elements 36 in the middle.₁, 36₂Considering a pile 32 having a nodal element 36₁, Node element 36₂The ground element that exerts the third contact force on the

ground spring

38₁, 38₂Represented by This node element 36₁And ground spring 38₁, Node element 36₂And ground spring 38₂The stable point that balances is the

nodal element

36₁, 36₂Is found by moving some of the nodal elements 36₁Δ₁ ⁽¹⁾Move only the nodal element 36₂Δ₂ ⁽¹⁾Assuming that the ground spring 38 is balanced by moving only the following equation, the following equation is established.
[0049]
P₁ ^'(1)-P₁ ⁽¹⁾= K₁δ₁ ⁽¹⁾ ... (3)
P₂ ^'(1)-P₂ ⁽¹⁾= K₂δ₂ ⁽¹⁾ ... (4)
However, P₁ ⁽¹⁾Is the nodal element 36 before moving.₁Force acting on P₁ ^'(1)Is the nodal element 36 after movement.₁It is the force acting on. Similarly, P₂ ⁽¹⁾Is the nodal element 36 before moving.₂Force acting on P₂ ^'(1)Is the nodal element 36 after movement.₂It is the force acting on. k₁, K₂Is

ground spring

38₁, 38₂Of the ground spring constant. The calculation of the ground spring constant will be described later. This ground spring constant corresponds to the spring constant of the present invention.
[0050]

Node element

36₁, 36₂Since the rigidity of the pile 32 balances with the force from the ground after the movement, the relationship between the displacement of the pile 32 and the force from the ground is expressed by the following equation.
[0051]
[Equation 3]

[0052]
However, E is an elastic coefficient of the pile 32 and I is a cross-sectional second moment of the pile 32.
[0053]
The above formulas (3) to (6) are summarized as follows.
[0054]
[Expression 4]

[0055]
The above equation (8)₁ ^'(1)Is P₂ ^'(1)By solving for, the force from the ground at the stable point can be calculated. Then, by substituting this into the equations (5) and (6) above, the

nodal element

36₁, 36₂You can ask how much you want to move.
[0056]
However, the node element 36 is not added to the point thus determined.₁, 36₂Since the interaction between the elements by DEM is non-linear even if it is moved, the pile 32 is rarely stabilized at the moved position. This corresponds to the fact that the actual ground spring is non-linear. For this reason, it is necessary to continue convergence calculation to find a stable point.
[0057]
Here, as in the first calculation described above, the nodal element 36₁Δδ₁ ⁽²⁾Move only the nodal element 36₂Δδ₂ ⁽²⁾Assuming that the ground spring 38 is balanced by moving only the following equation, the following equation is established.
[0058]
P₁ ^'(2)-P₁ ⁽²⁾= K₁Δδ₁ ⁽²⁾      ... (9)
P₂ ^'(2)-P₂ ⁽²⁾= K₂Δδ₂ ⁽²⁾      (10)
By the way, Δδ_i ⁽²⁾And the node element 36 at that time_iDeflection with δ_i ⁽²⁾Have the following relationship.
[0059]
δ_i ⁽²⁾= Δ_i ⁽¹⁾+ Δδ_i ⁽²⁾    ... (11)
Further, when the ground spring constant is calculated by the first calculation method described later, the node element 36 calculated by the first iterative calculation is used._iAnd the nodal element 36 calculated by the second iterative calculation_iThe difference between the force acting on the node element 36 and the node element 36 moved for the first time_iThe relationship is as follows.
[0060]
P_i ⁽²⁾-P_i ⁽¹⁾= K_iδ_i ⁽¹⁾      (12)
Substituting the above equations (11) and (12) into the above equations (9) and (10), the following relationship is obtained.
[0061]
[Equation 5]

[0062]
[Formula 6]

[0063]
These are summarized as follows.
[0064]
[Expression 7]

[0065]
The above equation (15) is changed to P₁ ^'(2), P₂ ^'(2)By solving for, the force from the ground at the stable point can be calculated. Then, by substituting this into the following equation, the

nodal element

36₁, 36₂You can ask where to move.
[0066]
[Equation 8]

[0067]
When the ground spring constant is calculated by the second calculation method described later, the relationship of the above equation (12) does not hold, so the above equation (15) has the following form.
[0068]
[Equation 9]

[0069]
Similarly, the convergence calculation is repeated, and it is determined that a stable point has been found when the following conditions are satisfied.
[0070]
[Expression 10]

[0071]
Next, calculation of the displacement of each part of the pile 32 in a state in which displacement is forcibly applied to the pile head of the pile 32, that is, the nth step (n is a natural number) of the time loop will be described.
[0072]
Here, as shown in FIG. 7, a pile 32 having one node element 36 in the middle is considered, and a ground element that exerts a third contact force on the node element 36 is represented by a ground spring 38.
[0073]
Considering the same as in the first step, assuming that the node element 36 and the ground spring 38 are balanced, the following equation is established.
[0074]
P^'(1)-P⁽¹⁾= KΔδ_t ⁽¹⁾ ... (20)
Where Δδ_t ⁽¹⁾Is the amount of movement of the nodal element 36, and this Δδ_t ⁽¹⁾The relationship between the stake and the deflection of the pile 32 is as follows from FIG.
[0075]

If this is substituted into the above equation (20), the following relationship is obtained.
[0076]

When the pile 32 is moved, the rigidity of the pile 32 balances with the force from the ground. Therefore, when there are a plurality of such nodal elements 36, the relationship between the force from the ground and the deflection of the pile 32 is the force vector P^'(1)And deflection vector δ_t ⁽¹⁾Is represented by the following formula.
[0077]
## EQU11 ##

[0078]
Where [H₀] Is a matrix representing the relationship between the force acting on the node element 36 and the deflection. For example, as shown in FIG.₁, 36₂In the case of a pile 32 having
[0079]
[Expression 12]

[0080]
[Formula 13]

[0081]
The above equation (25) is changed to P^'(1)By solving for, the force from the ground at the stable point can be calculated. And the calculated P^'(1)Is substituted into the above equation (23), it is possible to determine where the nodal element 36 should be moved.
[0082]
After the second iteration calculation in the n-th step convergence loop, a stable point can be obtained by the same algorithm as the second iteration computation after the first step convergence loop.
[0083]
In addition, regarding the process of the nth step which does not give displacement to the pile 32 compulsorily, y_t= Y_t-Δ_tAs described above, a stable point can be obtained by the same processing as the above algorithm.
[0084]
Next, calculation of the ground spring constant of the ground spring 38 will be described. As the calculation method of the ground spring constant, in the iterative calculation by the convergence loop, as shown in FIG. 8, the distance from which the nodal element 36 is moved by the deformation of the pile, and the ground element 30 from the nodal element 36 is moved. The first calculation method obtained based on the change in force and the ground of the ground spring 38 by combining the individual spring constants of the springs 40 representing the individual ground elements 30 in contact with the pile 32 as shown in FIG. There is a second method for calculating the spring constant.
[0085]
First, the first calculation method will be described.
[0086]
In the first step of the time loop, since the pile 32 is not moved in the first calculation of the convergence loop, it is based on the distance moved by the node element 36 and the change in the force from the ground due to the movement of the node element 36. Cannot be requested. In the second iteration calculation by the convergence loop, the nodal element 36₁Displacement δ₁ ⁽¹⁾By giving the node element 36₁The force applied to P is P₁ ⁽¹⁾To P₁ ⁽²⁾To change. Therefore, from the Hooke's law, the ground spring constant k₁Is obtained by changing the deformation of the pile 32, that is, the nodal element 36, as shown in FIG.₁Displacement and nodal element 36₁It is equivalent to finding the slope representing the rate of change of force from the change in force applied to the. As a result, the ground spring constant k₁Can be obtained by the following equation.
[0087]

This relationship is equal to the above equation (12) assumed in the first step algorithm.
[0088]
In the first calculation by the convergence loop of the nth step, using the amount of movement of the node element 36 from t−2Δt to t−Δt and the amount of change in the force applied to the node element 36 from the ground in the meantime, The spring constant k can be obtained.
[0089]
[Expression 14]

[0090]
The ground spring constant obtained in this way reflects the actual ground condition as it is, and when the ground spring is considered to be a relatively good quality nonlinear spring, it can be considered as a secant synthesis from a certain point to a certain point. .
[0091]
As described above, in each calculation step, the ground spring constant is obtained from the displacement of the nodal element 36 and the change in the force applied to the nodal element 36 at that time, so that the state of the ground structure can be analyzed with higher accuracy.
[0092]
Next, the second calculation method will be described.
[0093]
In DEM, elements transmit contact forces via springs in the normal and tangential directions. In the second calculation method, the spring constant of this spring is directly converted to the ground spring constant. Assuming that a plurality of ground elements 30 are in contact with the cylindrical element 34 and the node element 36 and the springs of these elements are connected in parallel, the spring constants (individual spring constants) of these springs are added together. Obtain the ground spring constant. This is done at each calculation step.
[0094]
Since the movement of each element and the movement of the nodal element are calculated for each direction component, a method for decomposing the spring of each element into each direction component will be described below. It is assumed that the positional relationship with the ground element 30 that is in contact with the node element 36 is as shown in FIG. The spring constant of this element is k, and the resolved spring constant in the x direction is K._xAnd When θ << 1, the following relationship is obtained.
[0095]

Moreover, it turns out that it becomes the following relationships from FIG.
[0096]
[Expression 15]

[0097]
As described above, in the first calculation method and the second calculation method, in each calculation step, the ground spring constant is obtained from the contact force between the pile 32 and the ground element 30 at that time. It can be performed. Further, as a result of the analysis, a relationship between the displacement of the pile 32 and the ground spring constant as shown in FIG. 12 can be obtained.
[0098]
After determining the ground spring constant and determining the displacement of the nodal element 36 in this way, in step 108, the convergence determining unit 20 determines the displacement of the nodal element 36 calculated by the displacement calculating unit 18 in a predetermined convergence condition. Whether or not the above equation (19) is satisfied is determined.
[0099]
If the displacement has not converged and the convergence condition is not satisfied, the determination in step 108 is denied, the process proceeds to step 102, and the processes in steps 102 to 108 are repeated.
[0100]
On the other hand, when the displacement converges and the convergence condition is satisfied, the determination in step 108 is affirmed and the process proceeds to step 110.
[0101]
In step 110, it is determined whether or not the movement of the pile head of the pile 32 that has been forcibly displaced has stopped. And when the pile head has not stopped, judgment of Step 110 is denied and it shifts to Step 100, and processing of the next time step is performed in the same manner as above by a time loop. On the other hand, if the pile head stops, the determination in step 110 is affirmed, and this routine ends.
[0102]
Thus, in this Embodiment, since the ground element 30 which contacts the pile 32 is represented by a ground spring and analyzed, the contact force newly generated when the pile 32 is deformed is considered by the expansion and contraction of the ground spring. The Therefore, the number of repeated calculations in the convergence loop can be significantly reduced compared to the conventional case.
[0103]
For example, DEM is an analysis method that tracks the movement of elements from time to time in units of a short time interval of 1 / 10,000 seconds to 1 / 100,000 seconds. And the analysis was very time consuming.
[0104]
Further, in the hybrid analysis of DEM and other analysis methods, the above-described balance calculation must be performed for each time step. If 100 repetitions are necessary, analysis is performed for 1 million time steps. This requires 100 million calculations overall.
[0105]
On the other hand, if the present invention is used, since the number of repetitions is from several to 10 or less, the calculation time can be greatly shortened and the calculation result can be obtained quickly.
[0106]
FIG. 13 shows the result of analyzing the behavior of the pile 32 according to the present invention, and FIG. 14 shows the result of analysis when the time loop process is advanced in a state including an error that does not satisfy the convergence condition. . In both FIG. 13 and FIG. 14, the horizontal axis indicates the moving distance of the pile head of the pile 32, and the vertical axis indicates the force applied to the pile 32.
[0107]
As can be seen from FIG. 14, when the time loop is advanced with an error remaining, the force applied to the pile 32 varies greatly, and it can be seen that the overall analysis is not stable. On the other hand, as is apparent from FIG. 13, in the analysis according to the present invention, it is understood that the fluctuation of the force applied to the pile 32 is small and the entire analysis is stable.
[0108]
【The invention's effect】
As described above, according to the present invention, it is possible to greatly reduce the calculation time and to analyze the behavior of the ground structure with high accuracy.
[Brief description of the drawings]
FIG. 1 is a schematic block diagram of a ground structure analyzing apparatus.
FIG. 2 is a diagram showing an analysis model of a ground structure.
FIG. 3 is a diagram illustrating models of ground elements and pile elements.
FIG. 4 is a flowchart showing the flow of processing by the ground structure analysis apparatus.
FIG. 5 is an image diagram for explaining aggregation of second contact force.
FIG. 6 is a diagram for explaining displacement of a nodal element.
FIG. 7 is a diagram for explaining displacement of a nodal element.
FIG. 8 is a diagram for explaining calculation of a ground spring constant in the first calculation method.
FIG. 9 is a diagram for explaining calculation of a ground spring constant in the second calculation method.
FIG. 10 is a diagram for explaining the calculation of the ground spring constant in the first calculation method.
FIG. 11 is a diagram for explaining calculation of a ground spring constant in the second calculation method.
FIG. 12 is a graph showing the relationship between pile head displacement and ground spring constant according to the analysis of the present invention.
FIG. 13 is a graph showing the relationship between the displacement of the pile head and the force applied to the pile according to the analysis of the present invention.
FIG. 14 is a graph showing the relationship between the displacement of the pile head and the force applied to the pile according to the conventional analysis.
[Explanation of symbols]
10 Ground structure analysis device
12 Contact force calculator
14 Contact force calculator
16 Contact force concentrator
18 Displacement calculator
20 Convergence judgment unit
30 Ground elements
31 Ground element
32 piles
34 Cylindrical elements
36 node elements
37 Analysis area
38 Ground spring
40 Spring

Claims

複数の粒状体である個別要素と当該複数の個別要素中に設けられた連続体である構造物との挙動を解析する粒状体及び連続体の解析方法であって、
前記複数の個別要素同士の第１の接触力を個別要素法により各々算出する第１ステップと、
複数の構造物要素及び当該複数の構造物要素同士を接続する節点要素から構成された前記構造物と前記個別要素との第２の接触力を個別要素法により各々算出する第２ステップと、
前記第２の接触力の各々を前記節点要素に集約した第３の接触力を算出する第３ステップと、
前記構造物要素及び前記節点要素に前記第２の接触力を及ぼす個別要素をバネとして表した時の前記バネのバネ定数を前記第２の接触力に基づいて算出する第４ステップと、
前記第３の接触力と前記バネ定数とに基づいて、前記節点要素と前記バネとが釣り合うと仮定した時の前記節点要素の変位を算出する第５ステップと、
前記節点要素と前記バネとが釣り合うと仮定した時の前記節点要素の変位が予め定めた収束条件を満たさない場合には、前記第５ステップで求めた変位のもとで前記収束条件を満たすまで前記第２ステップ〜第５ステップの処理を繰り返し、前記変位が予め定めた収束条件を満たす場合には、前記第５ステップで求めた変位のもとで前記第１ステップ〜第５ステップの処理を繰り返すことを特徴とする粒状体及び連続体の解析方法。A method of analyzing a granular material and a continuum for analyzing the behavior of an individual element that is a plurality of granular materials and a structure that is a continuum provided in the plurality of individual elements,
A first step of calculating a first contact force between the plurality of individual elements by an individual element method;
A second step of calculating, by an individual element method, a second contact force between the structure composed of a plurality of structure elements and node elements connecting the plurality of structure elements and the individual elements;
A third step of calculating a third contact force in which each of the second contact forces is aggregated into the node element;
A fourth step of calculating a spring constant of the spring based on the second contact force when the individual element that exerts the second contact force on the structure element and the node element is represented as a spring;
A fifth step of calculating a displacement of the nodal element when it is assumed that the nodal element and the spring are balanced based on the third contact force and the spring constant;
When the displacement of the nodal element when it is assumed that the nodal element and the spring are balanced does not satisfy the predetermined convergence condition, the convergence condition is satisfied under the displacement obtained in the fifth step. When the processing of the second step to the fifth step is repeated and the displacement satisfies a predetermined convergence condition, the processing of the first step to the fifth step is performed under the displacement obtained in the fifth step. A method of analyzing a granular material and a continuum characterized by repeating.

前記バネ定数は、前記節点要素の変位量と前記第３の接触力の変化量とから求めることを特徴とする請求項１記載の粒状体及び連続体の解析方法。 2. The granular body and continuum analysis method according to claim 1, wherein the spring constant is obtained from a displacement amount of the nodal element and a change amount of the third contact force.

前記バネ定数は、前記構造物要素及び前記節点要素に接触する前記個別要素の個別バネ定数を前記節点要素に集約することにより求めることを特徴とする請求項１記載の粒状体及び連続体の解析方法。 The granular body and continuum analysis according to claim 1, wherein the spring constant is obtained by aggregating the individual spring constants of the individual elements contacting the structure element and the nodal element into the nodal elements. Method.