2. Bascule non-intrusive 1D-3D#
Une bascule d’un modèle poutre à un modèle mixte poutre-3D permet de gagner considérablement en temps de calcul pour une précision équivalente à un modèle 3D entier [bib2].
Image 1.1-a fonctions de pondération dans le cadre Arlequin
Code_Aster permet d’initialiser un calcul dynamique avec des champs de déplacements, vitesses et accélérations non nulles (ETAT_INIT dans les opérateurs de calcul dynamique). Ces champs d’initialisation peuvent être préalablement calculés, puis transformés en un champ dans Code_Aster .
2.1. Principe théorique de la bascule#
La bascule consiste à passer d’un modèle de Poutre ( Up est la solution poutre du système \({M}_{p}\ddot{{U}_{p}}+{C}_{p}\dot{{U}_{p}}+{K}_{p}{U}_{p}={f}_{p}\) ) à un modèle 3D (solution U3D à initialiser) par l’intermédiaire d’une initialisation pertinente de la solution 3D. Le principe de base est le suivant :
à partir de la solution poutre Up , créer une solution 3D PUp pour l’hypothèse de section rigide
écrire le champs déplacement 3D comme étant la somme du champ extrudé de type Poutre Pup et d’un champ correcteur en déplacement noté U3Dc , prenant en compte la déformation de la section :
\({U}_{\mathrm{3D}}={\mathit{PU}}_{p}+{U}_{\mathrm{3Dc}}\)
La dynamique du modèle 3D s’écrira alors comme suit :
\({M}_{\mathrm{3D}}(\ddot{{\mathit{PU}}_{p}}+\ddot{{U}_{\mathrm{3Dc}}})+{C}_{\mathrm{3D}}(\dot{{\mathit{PU}}_{p}}+\dot{{U}_{\mathrm{3Dc}}})+{K}_{\mathrm{3D}}({\mathit{PU}}_{p}+{U}_{\mathrm{3Dc}})={f}_{\mathrm{3D}}\)
corriger les déplacements, en statique , sous l’effet d’un vecteur force f 3Dc calculé comme suit :
\({f}_{\mathrm{3Dc}}={K}_{\mathrm{3D}}{U}_{\mathrm{3Dc}}={f}_{\mathrm{3Dc}}(t={t}_{b})-{M}_{\mathrm{3D}}\ddot{{\mathit{PU}}_{p}}-{C}_{\mathrm{3D}}\dot{{\mathit{PU}}_{p}}-{K}_{\mathrm{3D}}{\mathit{PU}}_{p}\)
La correction statique est effectuée sur 3 pas de temps successifs : l’instant de bascule tbainsi que celui d’avant (tb-1 = tb- Δt) et celui d’après (tb+1 = tb+ Δt) où Δt est le pas de temps de calcul. On en déduit une correction en vitesse qui s’écrit comme suit :
\(\dot{{U}_{\mathrm{3D}}}=\frac{{\left[{\mathit{PU}}_{p}+{U}_{\mathrm{3Dc}}\right]}_{({t}_{b+1})}-{\left[{\mathit{PU}}_{p}+{U}_{\mathrm{3Dc}}\right]}_{({t}_{b-1})}}{2\Delta t}\) (Différences Finies Centrées)
Les accélérations, quant à elles, sont soit calculées de la même façon selon un schéma de type Différences Finies Centrées,
\(\ddot{{U}_{\mathrm{3D}}}=\frac{{\left[{\mathit{PU}}_{p}+{U}_{\mathrm{3Dc}}\right]}_{({t}_{b-1})}+{\left[{\mathit{PU}}_{p}+{U}_{\mathrm{3Dc}}\right]}_{({t}_{b+1})}-2{\left[{\mathit{PU}}_{p}+{U}_{\mathrm{3Dc}}\right]}_{({t}_{b})}}{{\Delta t}^{2}}\)
soit déterminées de façon automatique par Code_Aster dans le cas de l’utilisation d’un schéma d’intégration implicite en temps.
La correction comporte une correction de la flèche ainsi qu’une correction pour prendre en compte la déformation de la section.
2.2. Exemple de mise en œuvre#
La méthodologie est explicitée ci-dessous pour un problème de poutre bi-appuyée soumise à une excitation dynamique sinusoïdale durant 3s.
La simulation démarre par un modèle de poutre. Le calcul dynamique sur le mono-modèle de poutre a lieu de l’instant initial t0 jusqu’à l’instant de bascule au modèle 3D (ici, t0=0s, tb=2.4s).
Ti=0;
Tf=3.0;
Tb=2.4;
dt=0.0015;
Tb_1=Tb-dt;
Tb_2=Tb+dt;
Les opérations suivantes seront utilisées pour la bascule.
# Création du maillage Python vide
mail_py = MAIL_PY()
# Conversion du maillage Aster en maillage Python
mail_py.FromAster(“Mail”)
# Coordonnées des nœuds
COORD_3D = mail_py.cn
iPouAllNode = mail_py.gno.get(“AllNode”)
iPouRP1 = mail_py.gno.get(“RP1”)
iPouRP2 = mail_py.gno.get(“RP2”)
C3DZ = zeros((len(COORD_3D),1))
for i in range(len(C3DZ)):
C3DZ[i] = COORD_3D[i][2]
C1DZ1 = []
min=1112
while (min!=1111):
i = 0
s = int((abs(C3DZ[i])+ajust)*Pow)/float(Pow)
for i in range(len(C3DZ)):
if (int((C3DZ[i]+ajust)*Pow)/float(Pow))==min:
C3DZ[i] = 1111
elif( (int((C3DZ[i]+ajust)*Pow)/float(Pow)) < s ):
s = int((C3DZ[i]+ajust)*Pow)/float(Pow)
if (min!=s):
min=s
C1DZ1.append(min)
C1D = zeros((len(C1DZ1)-1,3))
for i in range(len(C1D)):
C1D[i][0] = 0
C1D[i][1] = 0
C1D[i][2] = C1DZ1[i]
# Fonction pour chercher la section qui correspond a une valeur de z
def OrdreSect(z):
jk=0
while(C1D[jk][2]!=z):
jk=jk+1
return jk
A l’instant de bascule tb, on récupère les déplacements Up, les vitesses Vp et accélérations Ap de la fibre neutre de la poutre par l’opérateur CREA_CHAMP, opération “EXTR”. Ensuite, les composantes (DX,DY,DZ,DRX,DRY,DRZ) sont extraites dans des tables Python.
Temps = NP.array([Tb,Tb_1,Tb_2])
nbpas = temps.shape[0]
U0=[]
U1=[]
U2=[]
for i in range(nbpas):
pas = temps[i]
DEP=CREA_CHAMP(OPERATION=”EXTR”,
TYPE_CHAM=”NOEU_DEPL_R”,
RESULTAT=DLT,
NOM_CHAM=”DEPL”,
INST=pas,);
Upx = DEP.EXTR_COMP(“DX”,[“AllNode”]).valeurs;
Upy = DEP.EXTR_COMP(“DY”,[“AllNode”]).valeurs;
Upz = DEP.EXTR_COMP(“DZ”,[“AllNode”]).valeurs;
Uprx = DEP.EXTR_COMP(“DRX”,[“AllNode”]).valeurs;
Upry = DEP.EXTR_COMP(“DRY”,[“AllNode”]).valeurs;
Uprz = DEP.EXTR_COMP(“DRZ”,[“AllNode”]).valeurs;
ACC=CREA_CHAMP(OPERATION=”EXTR”,
TYPE_CHAM=”NOEU_DEPL_R”,
RESULTAT=DLT,
NOM_CHAM=”ACCE”,
INST=pas,);
ddotUpx = ACC.EXTR_COMP(“DX”,[“AllNode”]).valeurs;
ddotUpy = ACC.EXTR_COMP(“DY”,[“AllNode”]).valeurs;
ddotUpz = ACC.EXTR_COMP(“DZ”,[“AllNode”]).valeurs;
ddotUprx = ACC.EXTR_COMP(“DRX”,[“AllNode”]).valeurs;
ddotUpry = ACC.EXTR_COMP(“DRY”,[“AllNode”]).valeurs;
ddotUprz = ACC.EXTR_COMP(“DRZ”,[“AllNode”]).valeurs;
Par une hypothèse de type corps rigide, on calcule des tables Python contenant les déplacements PUp, les vitesses PVp et accélérations PAp du modèle mixte par extrusion des champs poutre dans l’épaisseur.
Ensuite, des tables Python PUp et ddotPUp sont calculées et contiennent les champs de déplacement et d’accélération des nœuds du modèle mixte en fonction de la topologie (1D ou 3D).
dictLNo = []
dictPUpx = []
dictPUpy = []
dictPUpz = []
dictPUprx = []
dictPUpry = []
dictPUprz = []
dictAccx = []
dictAccy = []
dictAccz = []
dictAccrx = []
dictAccry = []
dictAccrz = []
for i in range(len(PUp)):
j=i+1
if len(PUp[i])==3:
dictLNo.append(“N%i”%j);
dictPUpx.append(PUp[i][0]);
dictPUpy.append(PUp[i][1]);
dictPUpz.append(PUp[i][2]);
dictAccx.append(ddotPUp[i][0]);
dictAccy.append(ddotPUp[i][1]);
dictAccz.append(ddotPUp[i][2]);
else:
dictLNo.append(“N%i”%j);
dictPUpx.append(PUp[i][0]);
dictPUpy.append(PUp[i][1]);
dictPUpz.append(PUp[i][2]);
dictPUprx.append(PUp[i][3]);
dictPUpry.append(PUp[i][4]);
dictPUprz.append(PUp[i][5]);
dictAccx.append(ddotPUp[i][0]);
dictAccy.append(ddotPUp[i][1]);
dictAccz.append(ddotPUp[i][2]);
dictAccrx.append(ddotPUp[i][3]);
dictAccry.append(ddotPUp[i][4]);
dictAccrz.append(ddotPUp[i][5]);
Ces tables Python sont par la suite transformées en listes puis en champs de Code_Aster.
TD=CREA_TABLE(LISTE=(_F(LISTE_K=dictLNo,PARA=”NOEUD”),
_F(LISTE_R=dictPUpx,PARA=”DX”),
_F(LISTE_R=dictPUpy,PARA=”DY”),
_F(LISTE_R=dictPUpz,PARA=”DZ”),
_F(LISTE_R=dictPUprx,PARA=”DRX”),
_F(LISTE_R=dictPUpry,PARA=”DRY”),
_F(LISTE_R=dictPUprz,PARA=”DRZ”),))
TA=CREA_TABLE(LISTE=(_F(LISTE_K=dictLNo,PARA=”NOEUD”),
_F(LISTE_R=dictAccx,PARA=”DX”),
_F(LISTE_R=dictAccy,PARA=”DY”),
_F(LISTE_R=dictAccz,PARA=”DZ”),
_F(LISTE_R=dictAccrx,PARA=”DRX”),
_F(LISTE_R=dictAccry,PARA=”DRY”),
_F(LISTE_R=dictAccrz,PARA=”DRZ”),))
DeplPUp=CREA_CHAMP(TYPE_CHAM=”NOEU_DEPL_R”,
OPERATION=”EXTR”,
MAILLAGE=Mail,
TABLE=TD,);
AccIni=CREA_CHAMP(TYPE_CHAM=”NOEU_DEPL_R”,
OPERATION=”EXTR”,
MAILLAGE=Mail,
TABLE=TA,);
Les produits matrices-champs K3D*PUp et M3D*PAp sont effectués par PROD_MATR_CHAM.
FORCE1=PROD_MATR_CHAM(MATR_ASSE=RIGIDITE,
CHAM_NO=DeplPUp,
TITRE=”PROD_MATR_CHAM1”,);
TBLf1=POST_RELEVE_T(ACTION=_F(OPERATION=”EXTRACTION”,
INTITULE=”f1”,
CHAM_GD=FORCE1,
GROUP_NO=”AllNode”,
NOM_CMP=(“DX”,”DY”,”DZ”,”DRX”,”DRY”,”DRZ”,),),);
VARf1=TBLf1.EXTR_TABLE()
FORCE2=PROD_MATR_CHAM(MATR_ASSE=MASSE,
CHAM_NO=AccIni,
TITRE=”PROD_MATR_CHAM2”,);
TBLf2=POST_RELEVE_T(ACTION=_F(OPERATION=”EXTRACTION”,
INTITULE=”f2”,
CHAM_GD=FORCE2,
GROUP_NO=”AllNode”,
NOM_CMP=(“DX”,”DY”,”DZ”,”DRX”,”DRY”,”DRZ”,),),);
VARf2=TBLf2.EXTR_TABLE()
En fonction des forces en présence dans le problème traité, le chargement f3D(tb)-K3D*PUp-M3D*PAp est calculé sous forme d’un dictionnaire par nœuds puis puis stocké comme suit :
CHARGE3D=AFFE_CHAR_MECA(MODELE=MODELE,
FORCE_NODALE=(dictCharge),);
Ensuite, un champ de correction de la déformation de la section U3Dc est calculé par MECA_STATIQUE comme la réponse, à l’instant de bascule, du modèle 3D de référence au chargement f3D(tb)-K3D*PUp-M3D*PAp .
Mstat=MECA_STATIQUE(MODELE=MODELE,
CHAM_MATER=CHMAT,
CARA_ELEM=Carel3D,
EXCIT=(_F(CHARGE=CondLim,),
_F(CHARGE=CHARGE3D,),),);
La réponse qui en résulte est relevée puis stockée dans une table :
tableUc=POST_RELEVE_T(ACTION=_F(OPERATION=”EXTRACTION”,
INTITULE=”U3Dc”,
RESULTAT=Mstat,
NOM_CHAM=”DEPL”,
PRECISION=1e-15,
GROUP_NO=”AllNode”,
NOM_CMP=(“DX”,”DY”,”DZ”, “DRX”, “DRY”, “DRZ”,),),);
TC=tableUc.EXTR_TABLE()
Enfin, le champ de déplacement du patch 3D du modèle mixte à l’instant tb sera la somme des deux champs U3Dc et PUp .
# dictionnaire des déplacements au moment de la bascule
dictDeplx = []
dictDeply = []
dictDeplz = []
dictDeplrx = []
dictDeplry = []
dictDeplrz = []
for i in range(len(U0)):
j=i+1
if len(U0[i])==3:
dictDeplx.append(U0[i][0]);
dictDeply.append(U0[i][1]);
dictDeplz.append(U0[i][2]);
else:
dictDeplx.append(U0[i][0]);
dictDeply.append(U0[i][1]);
dictDeplz.append(U0[i][2]);
dictDeplrx.append(U0[i][3]);
dictDeplry.append(U0[i][4]);
dictDeplrz.append(U0[i][5]);
Dini=CREA_TABLE(LISTE=(_F(LISTE_K=dictLNo,PARA=”NOEUD”),
_F(LISTE_R=dictDeplx,PARA=”DX”),
_F(LISTE_R=dictDeply,PARA=”DY”),
_F(LISTE_R=dictDeplz,PARA=”DZ”),
_F(LISTE_R=dictDeplrx,PARA=”DRX”),
_F(LISTE_R=dictDeplry,PARA=”DRY”),
_F(LISTE_R=dictDeplrz,PARA=”DRZ”),))
Dep0=CREA_CHAMP(TYPE_CHAM=”NOEU_DEPL_R”,
OPERATION=”EXTR”,
MAILLAGE=Mail,
TABLE=Dini,);
Les champs de vitesse et accélération sont également calculées selon un schéma de Différences Centrées, en utilisant les champs Poutre et corrigés des déplacements, vitesses et accélérations à l’instant de bascule tb et aux instants tb – Δt et tb + Δt .
# Calcul des vitesses corrigées
dictVitex = []
dictVitey = []
dictVitez = []
dictViterx = []
dictVitery = []
dictViterz = []
for i in range(len(U1)):
j=i+1
if len(U1[i])==3:
dictVitex.append((U2[i][0]-U1[i][0])/(2*dt));
dictVitey.append((U2[i][1]-U1[i][1])/(2*dt));
dictVitez.append((U2[i][2]-U1[i][2])/(2*dt));
else:
dictVitex.append((U2[i][0]-U1[i][0])/(2*dt));
dictVitey.append((U2[i][1]-U1[i][1])/(2*dt));
dictVitez.append((U2[i][2]-U1[i][2])/(2*dt));
dictViterx.append((U2[i][3]-U1[i][3])/(2*dt));
dictVitery.append((U2[i][4]-U1[i][4])/(2*dt));
dictViterz.append((U2[i][5]-U1[i][5])/(2*dt));
Vini=CREA_TABLE(LISTE=(_F(LISTE_K=dictLNo,PARA=”NOEUD”),
_F(LISTE_R=dictVitex,PARA=”DX”),
_F(LISTE_R=dictVitey,PARA=”DY”),
_F(LISTE_R=dictVitez,PARA=”DZ”),
_F(LISTE_R=dictViterx,PARA=”DRX”),
_F(LISTE_R=dictVitery,PARA=”DRY”),
_F(LISTE_R=dictViterz,PARA=”DRZ”),))
Vit0=CREA_CHAMP(TYPE_CHAM=”NOEU_DEPL_R”,
OPERATION=”EXTR”,
MAILLAGE=Mail,
TABLE=Vini,);
# Calcul des accélérations corrigées
dictAccex = []
dictAccey = []
dictAccez = []
dictAccerx = []
dictAccery = []
dictAccerz = []
for i in range(len(U1)):
j=i+1
if len(U1[i])==3:
dictAccex.append((U2[i][0]-2*U0[i][0]+U1[i][0])/(dt**2));
dictAccey.append((U2[i][1]-2*U0[i][1]+U1[i][1])/(dt**2));
dictAccez.append((U2[i][2]-2*U0[i][2]+U1[i][2])/(dt**2));
else:
dictAccex.append((U2[i][0]-2*U0[i][0]+U1[i][0])/(dt**2));
dictAccey.append((U2[i][1]-2*U0[i][1]+U1[i][1])/(dt**2));
dictAccez.append((U2[i][2]-2*U0[i][2]+U1[i][2])/(dt**2));
dictAccerx.append((U2[i][3]-2*U0[i][3]+U1[i][3])/(dt**2));
dictAccery.append((U2[i][4]-2*U0[i][4]+U1[i][4])/(dt**2));
dictAccerz.append((U2[i][5]-2*U0[i][5]+U1[i][5])/(dt**2));
Aini=CREA_TABLE(LISTE=(_F(LISTE_K=dictLNo,PARA=”NOEUD”),
_F(LISTE_R=dictAccex,PARA=”DX”),
_F(LISTE_R=dictAccey,PARA=”DY”),
_F(LISTE_R=dictAccez,PARA=”DZ”),
_F(LISTE_R=dictAccerx,PARA=”DRX”),
_F(LISTE_R=dictAccery,PARA=”DRY”),
_F(LISTE_R=dictAccerz,PARA=”DRZ”),))
Acc0=CREA_CHAMP(TYPE_CHAM=”NOEU_DEPL_R”,
OPERATION=”EXTR”,
MAILLAGE=Mail,
TABLE=Aini);
Les champs d’initialisation ainsi construits sont utilisés pour initialiser le calcul dynamique (ETAT_INIT dans les opérateurs de calcul dynamique).
Bascule=DYNA_LINE_TRAN(MODELE=MODELE,
CHAM_MATER=CHMAT,
MATR_MASS=MASSE,
MATR_RIGI=RIGIDITE,
SCHEMA_TEMPS=_F(…),
ETAT_INIT=_F(DEPL=Dep0,
VITE=Vit0,
ACCE=Acc0,),
EXCIT=(_F(VECT_ASSE=F_Xass,
FONC_MULT=FSIN),),
INCREMENT=_F(LIST_INST=LIST,),);
La procédure non intrusive de bascule d’un modèle de poutre à un modèle 3D en dynamique transitoire, développée et validée dans Code_Aster (cas-test SDNV139), permet de gagner en temps de calcul car le modèle global de type structure mince peut remplacer le modèle local 3D en l’absence de non linéarités localisées en temps. Par ailleurs, l’utilisation du raccord 3D-POU préexistant dans Code_Aster permet de réduire la taille du modèle pour des cas de non linéarités localisées en espace.