CISAILLEMENT ET AUTRES CONTRAINTES DANS LA PLAQUE À TROIS COUCHES

Préliminaires

Initialisation

On a créé dans le "home" un dossier nommé "bibliothequeBN " et dans lequel s'enregistreront les fichiers à conserver.

Block[{racine},racine="~/bibliothequeBN";
Off[SetDelayed::write];
SetDirectory[racine];
$TextStyle={FontFamily->"Times",FontSize->14};
SetDirectory[StringJoin[racine,"/Fichiers plaques 2009"]]]

Off[General::"spell1"];

Commandes d'instanciation

(*Désinstanciation des raideurs "raid" et des coefficients "q"*)
inst0:=
Block[{},
If[ToString[raid[1,1]]≠"raid[1, 1]",Do[raid[i,j]= .;,{i,6},{j,6}]];
If[ToString[q[1]]≠"q[1]",Do[q[i]= .;,{i,3 (ddz+1) dixi}]]];
(*Instanciation des raideurs "raid"*)
inst1:=
Block[{},Do[raid[i,j]=k3d1[[i,j]],{i,6},{j,6}]];
inst2:=
Block[{},Do[raid[i,j]=k3d2[[i,j]],{i,6},{j,6}]];

Données générales

deg=2;ddz=5;

géométrie des  couches (couches de même épaisseur)

nc=3;nci=Ceiling[nc/2];ncp=nc-nci;
zc=Table[0,{i,nc+1}];
zc[[1]]=-0.5;
Do[zc[[n+1]]=zc[[n]]+If[OddQ[n],1/(2 nci),1/(2 ncp)],{n,nc}];

zc

Constantes mécaniques

Raideur formelle

k3d=Table[If[i<j,raid[i,j],raid[j,i]],{i,6},{j,6}];

Raideurs particulières pour un multicouche. 
Couches orthotropes (plan 1-2)- isotropes (plan 2-3)
(Batoz-Dhatt, 1-101)

Couches impaires : x est la direction des fibres

k3d1={{25,0.25,0.,0.,0.,0.},{0.25,1,0.,0.,0.,0.},{0.,0.,1,0.,0.,0.},{0.,0.,0.,0.5,0.,0.},{0.,0.,0.,0.,0.5,0},{0.,0.,0.,0.,0.,0.2}}

Couches paires : y est la direction des fibres

k3d2={{1,0.25,0.,0.,0.,0.},{0.25,25,0.,0.,0.,0.},{0.,0.,1,0.,0.,0.},{0.,0.,0.,0.5,0.,0.},{0.,0.,0.,0.,0.2,0.},{0.,0.,0.,0.,0.,0.5}}

inst1;Print[MatrixForm[k3d]];
inst2;Print[MatrixForm[k3d]];inst0;

Fonctions diverses

Suffixe pour les fichiers enregistrés

suffixe="_M7_xy"<>ToString[deg]<>"_z"<>ToString[ddz]<>"_c"<>ToString[nc]<>"_09.m";
Print[suffixe];

Raccourcis d'écriture

dim[x_]:=Dimensions[x];
dim1[x_]:=Dimensions[x][[1]];
dim2[x_]:=Dimensions[x][[2]];
codim[a_]:=Dimensions[NullSpace[a]][[1]];
tr[x_]=Transpose[x];
CE[x_]:=Chop[Expand[x]];

sub0={0.->0,1.->1,-1.->-1};

gradient d'une fonction scalaire de(x,y)

gradsca2[w_]:={D[w,x],D[w,y]}

Divergence d'un tenseur de contraintes représenté par une colonne de 6 éléments

Le tenseur des  contraintes est représenté par le vecteur : ten = {(Subscript[σ, 11]), (Subscript[σ, 22]), (Subscript[σ, 33]), (Subscript[σ, 12]), (Subscript[σ, 13]), (Subscript[σ, 23])}

div[ten_]:=
{D[ten[[1]],x]+D[ten[[4]],y]+D[ten[[5]],z],D[ten[[4]],x]+D[ten[[2]],y]+D[ten[[6]],z],D[ten[[5]],x]+D[ten[[6]],y]+D[ten[[3]],z]};

gradient symétrique d'un vecteur du plan

gradsy2[vec_]:={D[vec[[1]],x],D[vec[[2]],y],(D[vec[[1]],y]+D[vec[[2]],x])/2}

gradient d'un tenseur plan représenté par une colonne de 3 éléments

grad2[ten_]:={D[ten[[1]],x],D[ten[[1]],y],D[ten[[2]],x],D[ten[[2]],y],D[ten[[3]],x],D[ten[[3]],y]}

Divergence plane d'un tenseur  représenté par une colonne de 3 éléments

div2[ten_]:=
{D[ten[[1]],x]+D[ten[[3]],y],D[ten[[3]],x]+D[ten[[2]],y]};

Divergence plane d'un vecteur

divvec2[vec_]:=D[vec[[1]],x]+D[vec[[2]],y];

Opérateur de déformation

Le tenseur des  déformations est représenté par le vecteur : ε = {(Subscript[ε, 11]), (Subscript[ε, 22]), (Subscript[ε, 33]), 2(Subscript[ε, 12]), 2(Subscript[ε, 13]), 2(Subscript[ε, 23])}, d'où l'opérateur grads :

grads[u_]:=
{D[u[[1]],x],D[u[[2]],y],D[u[[3]],z],
D[u[[1]],y]+D[u[[2]],x],D[u[[1]],z]+D[u[[3]],x],D[u[[2]],z]+D[u[[3]],y]};

outils d'intégration en z

On intègre les monômes une fois pour toutes

intzz = Table[Integrate[z^(k-1), {z,zc[[n]], zc[[n+1]]}], {n, nc}, {k, ddz+2}];
intsub = Table[z^k->intzz[[n, k]], {n, nc}, {k, ddz+2}];
integ[a_] := Block[{}, Chop[Expand[Sum[(Expand[z a[[n]]])/.intsub[[n]], {n,nc}]]]];

Monômes en (x,y)

(*Liste des monômes en (x,y) de degré inférieur ou égal à n*)
monxy[n_]:=
Block[{pol,p,q},pol={1};Do[pol=Append[pol,x^(p-q) y^q],{p,n},{q,0,p}];pol];
mon=monxy[deg];
(*Liste des paires d'exposants correspondantes*)
exposant[n_]:=
Block[{pol,p,q},
pol={{0,0}};
Do[pol=Append[pol,{(p-q),q}],{p,n},{q,0,p}];pol];
(*nombre de ces monômes*)
dixi=dim1[mon];
dixi1=dim1[monxy[deg+1]];

Présentation des degrés de zeta

Au lieu de s[m], avec m=1 à 8 dixi,on écrit :
s[i,k,l] où i est le numéro de composante, k la puissance de x, l=j-k le degré de y.

tabs1=Flatten[Table[s1[i,k,j-k],{i,8},{j,0,deg},{k,j,0,-1}],2];
substizeta=Table[s[j]->tabs1[[j]],{j,8dixi}];
Print[substizeta];

Déplacements et déformations polynomiaux

Déplacements généralisés

membrane

m=0;vq=Table[Sum[m=m+1;m1=1+Mod[m-1,dixi1];q1[m]monxy[3][[m1]],{j,0,deg+1},{k,j,0,-1}],{i,2}]

flèche

m=0;wq=Table[Sum[m=m+1;m1=1+Mod[m-1,dixi1];q1[m+2dixi1]monxy[3][[m1]],{j,0,deg+1},{k,j,0,-1}],{i,1}]

inclinaison de la normale

m=0;
beta=Table[Sum[m=m+1;m1=1+Mod[m-1,dixi1];q1[m+3dixi1]monxy[3][[m1]],{j,0,deg+1},{k,j,0,-1}],{i,2}]

qgene={vq,wq,beta};

Déformations généralisées engendrées par les précédents

membrane

epsim=Expand[gradsy2[vq]]

courbure

chi=Expand[gradsy2[beta]]

cisaillement

gamma=Expand[gradsca2[w]+beta]

Déformations généralisées "abstraites"

membrane

m=0;epsim1=Table[Sum[m=m+1;m1=1+Mod[m-1,dixi];s[m]monxy[2][[m1]],{j,0,deg},{k,j,0,-1}],{i,3}]

courbure

m=3dixi;chi1=Table[Sum[m=m+1;m1=1+Mod[m-1,dixi];s[m]monxy[2][[m1]],{j,0,deg},{k,j,0,-1}],{i,3}]

cisaillement

m=6dixi;gamma1=Table[Sum[m=m+1;m1=1+Mod[m-1,dixi];s[m]monxy[2][[m1]],{j,0,deg},{k,j,0,-1}],{i,2}]

substitution des réelles aux abstraites

On remplace les déformations généralisées "abstraites" précédentes par leur expressions en fontion des déplacements généralisés q1

procédure d'identification polynomiale

idenpol[deg1_,pol1_,pol2_]:=Table[Coefficient[x y pol1,x y monxy[deg+1][[k]]]==Coefficient[x y pol2,x y monxy[deg1][[k]]],{k,dim1[monxy[deg1]]}]

substidef=
Module[{epsid, chiid,gammaid,sol1,sol2,sol3},
epsid=idenpol[2,epsim1,epsim];
sol1=Flatten[Expand[Solve[epsid,Table[s[i],{i,3dixi}]]]];
chiid=idenpol[2,chi1,chi];
sol2=Flatten[Expand[Solve[chiid,Table[s[i],{i,3dixi+1,6dixi}]]]];
gammid=idenpol[2,gamma1,gamma];
sol3=Flatten[Expand[Solve[gammid,Table[s[i],{i,6dixi+1,8dixi}]]]];
Sort[Flatten[{sol1,sol2,sol3}]]];

substidef

substidef>>"substidef"<>suffixe;

substidef1=Get["substidef"<>suffixe];
substidef==substidef1

Vérification de la solution

inst0; {uposol1,secmulti1} = Chop[Get["ucomp"<>suffixe]];

Déformations et contraintes produites par ce déplacement complémentaire

defco=grads[uposol1];
sigco=k3d.defco;
divsigco=Expand[div[sigco]];

Déformations principales

Les s[i], pour i de 1 à 8 dixi, sont les coefficients des déformations généralisées zeta mises sous forme polynomiales ;
membrane :	i=1,2,3 ;
flexion-torsion : i=4,5,6 ;
cisaillement : i=7,8 ;
defpr : déformations tridimensionnelles  (SuperDagger[B]) ζ de la théorie

zeta=Table[Sum[s[(i-1) dixi+j] mon[[j]],{j,dixi}],{i,8}];
defprabs={zeta[[1]]+z zeta[[4]],zeta[[2]]+z zeta[[5]],0,zeta[[3]]+z zeta[[6]],zeta[[7]],zeta[[8]]};
(*contraintes tridimensionnelles principales*)
sigpr=k3d.defprabs;
(*divergence des contraintes tridimensionnelles principales*)
divsigpr=Expand[div[sigpr]];
(* affichage éventuel, pour information, de la première composante de cette divergence formelle*)

Liaisons

Print[
Block[{integ1,integ2,integ3,integ4,integ5,u},u=Table[uposol1[[n,i]],{i,3},{n,nc}];
integ1=integ[u[[1]]];
integ2=integ[u[[2]]];
integ3=integ[u[[3]]];
integ4=integ[z u[[1]]];
integ5=integ[z u[[2]]];
{integ1,integ2,integ3,integ4,integ5}]];

Contraintes

sigma=
CE[Table[
If[OddQ[n],inst1,inst2];
k3d.(grads[uposol1[[n]]]+defprabs),
{n,nc}]];

divsigma=CE[tr[div[tr[sigma]]]];

Print[CE[sigma[[1,{3,5,6}]]/.z->zc[[1]]]];

Do[Print[CE[(sigma[[n,{3,5,6}]]-sigma[[n+1,{3,5,6}]])/.z->zc[[n+1]]]],
{n,nc-1}];

Print[CE[sigma[[nc,{3,5,6}]]/.z->zc[[nc+1]]]];

Équilibre intérieur

Do[
equ=divsigma[[n]]-secmulti1;
Print[Chop[Expand[equ]]],
{n,nc}];

Contraintes généralisées

sigtr=tr[sigma];
membrane=integ[Transpose[sigtr[[{1,2,4}]]]];
moment=integ[z Transpose[sigtr[[{1,2,4}]]]];
tranchant=integ[Transpose[sigtr[[{5,6}]]]];

intmul1=Chop[Integrate[secmulti1,{z,-0.5,0.5}]];
intmul2=Integrate[Chop[z secmulti1],{z,-1/2,1/2}];

Print[Chop[div2[membrane]-intmul1[[{1,2}]]]];
Print[Chop[div2[moment]-tranchant-intmul2[[{1,2}]]]];
Print[Chop[Expand[divvec2[tranchant]-intmul1[[3]]]]];

Profils des contraintes et des déplacements complémentaires
 dans une région non chargée et de déplacement de membrane nul

Annulation des multiplicateurs de Lagrange

multinul=Chop[Expand[secmulti1]];

ce qui conduit à la matrice d'égalité des polynômes à zéro :

matbrute=
Flatten[
Table[Coefficient[x y z multinul[[i]],x y z^k mon[[j]]s[p]],{i,3},{j,dixi},{k,ddz+1},{p,8 dixi}],2];
Print[dim[matbrute]];

On réduit cette matrice par enlèvement des lignes dépendantes et des colonnes nulles.Détermination de ces dernières.'colut' est la liste des colonnes non nulles. 'mat' est la matrice réduite.

lignenulles[matrice_]:=
Module[{mat,x,lignul},
mat=Chop[matrice];
lignul={};
{d1,d2}=dim[mat];
Do[If[mat[[i]]==Table[0,{d2}],lignul=Append[lignul,i]],
{i,d1}];
lignul];

lignenulles[matbrute]

dim[%]

redu1[system_]:=Module[{colut,d1,d2,la,ligut,mat,n1,n2,system0,system1,x},
n1=lignenulles[system];
n2=lignenulles[tr[system]];
{d1,d2}=dim[system];
colut=Complement[Range[1,d2],n2];
ligut=Complement[Range[1,d1],n1];
system1=Chop[system[[ligut,colut]]];
{d1,d2}=dim[system1];
mat={system1[[1]]};
Off[LinearSolve::"nosol"];
ligut1={1};
Do[
x=system1[[i]];
la=LinearSolve[Transpose[mat],x];
(*on utilise une propriété de 'LinearSolve' : s'il n'a pas trouvé de solution il renvoie un message 'la' qui contient l'expression "LinearSolve".*)
If[StringPosition[ToString[la],"LinearSolve"]!={},
mat=Append[mat,x];ligut1=Append[ligut1,i]],
{i,2,d1}];
(*réactivation des messages d'alerte*)
On[LinearSolve::"nosol"];1
{colut,mat}];

{colut,mat}=redu1[matbrute];Print[dim[mat]];

colut

dim[%]

Il faut maintenant choisir les inconnues c'est-à-dire les s[i] qu'on élimine. A cette fin on commence par compter le nombre d'éléments non nuls de chaque ligne de la matrice 'mat' et on classe les lignes par ordre croissant d'éléments non nuls, ce qui donne la liste 'ord'.
Puis on choisit, ligne par ligne dans l'ordre défini par 'ord', l'inconnue s[j] la plus à droite de numéro j le plus élevé qui soit concernée par cette ligne, et qui n'ait pas été choisie dans une ligne précédente.

inconn[matrice_]:=
Module[{d1,d2,compt,ord,inc,},
{d1,d2}=dim[mat];
compt=Table[0,{d1}];
Do[
compt[[i]]=d2-dim1[tr[Position[mat[[i]],0]][[1]]],
{i,d1}];
ord=Union[compt];
ord=Flatten[Table[tr[Position[compt,ord[[j]]]][[1]],{j,dim1[ord]}]];
inc={};
Do[
Do[
If[Not[Chop[mat[[ord[[i]],j]]]===0] &&Not[MemberQ[inc,j]],inc=Append[inc,j];Break[]],
{j,d2,1,-1}],
{i,d1}];Sort[inc]];

inc=inconn[mat];Print[inc];
inconnues=Table[s[colut[[inc[[j]]]]],{j,dim1[inc]}];Print[inconnues/.substizeta];

Résolution

vari=Table[s[colut[[j]]],{j,dim1[colut]}];
zero=Table[0,{dim1[mat]}];
substi2=Chop[Expand[Simplify[Solve[mat.vari==zero,inconnues][[1]]]]];
Print[substi2/.substizeta]

Print[Chop[Simplify[(mat.vari)/.substi2]]];

Quantités mécaniques après élimination des inconnues s(j) dépendantes des autres

uposol2=Chop[Expand[uposol1/.substi2]];
secmulti2=Chop[Expand[secmulti1/.substi2]];
sigma2=Chop[Expand[sigma/.substi2]];
divsigma2=CE[divsigma/.substi2];
zeta2=Chop[Simplify[zeta/.substi2]/.sub0];
membrane2=membrane/.substi2;
moment2=moment/.substi2;
tranchant2=tranchant/.substi2;

Vérification

Print[Chop[Simplify[div2[membrane2]]]]
Print[Chop[Simplify[div2[moment2]-tranchant2]]];
Print[Chop[Expand[divvec2[tranchant2]]]];
Print[Chop[secmulti2]];

Cas de la contrainte (Subscript[σ, 33])

CE[tr[sigma2][[3,1]]/.z->-0.5]

CE[tr[sigma2][[3,3]]/.z->0.5]

Ayant vérifié les équations d'équilibre et de liaison pour tout couple (x,y)on on se limite maintenant aux valeurs à l'origine après remplacements des déformations généralisées par leur expression en fonction des déplacements généralisés.

substi0={x->0,y->0};

On écrit que le déplacement principal v est nul tandis que w et gradw le sont aussi

subtot=Join[substi0,Table[q1[i]->0,{i,23}]]

substidef=Get["substidef"<>suffixe];

uposol3=CE[(uposol2/.substidef)/.subtot];
secmulti3=CE[(secmulti2/.substidef)/.subtot];
sigma3=CE[(sigma2/.substidef)/.subtot];
divsigma3=CE[(divsigma2/.substidef)/.subtot];
zeta3=CE[(zeta2/.substidef)/.subtot];
membrane3=CE[(membrane2/.substidef)/.subtot];
moment3=CE[(moment2/.substidef)/.subtot];
tranchant3=CE[(tranchant2/.substidef)/.subtot];
gradM3=CE[(grad2[moment2]/.substidef)/.subtot]

Variables q1 qui interviennent dans le gradient du moment

q1inc=Module[{a,b},
b={};
Do[a=Coefficient[gradM3[[j]],q1[k]];
If[a≠0 &&a≠0.,b=Append[b,q1[k]]],
{k,8dixi},{j,6}];
Union[b]]

Expression de ces variables en fonctions des composantes du gradient du moment

solq1=CE[Solve[gradM3=={M111,M112,M221,M222,M121,M122},q1inc][[1]]]

Relations inverses des précédentes

solM=CE[Solve[gradM3=={M111,M112,M221,M222,M121,M122},{M111,M112,M221,M222,M121,M122}][[1]]]

Expression de l'effort tranchant en fonction du gradient du moment

tranchant3

CE[tranchant3/.solq1]

Parfait !

sigpr1=sigpr;sigpr>>sigpr;

Expression des contraintes de cisaillement

Cas général

tr[sigma3][[{5,6}]]

Ce sont les coefficients du degré 2 dans beta qui interviennent

Coefficient[beta,q1[34]]

Coefficient[beta,q1[36]]

Coefficient[beta,q1[45]]

Coefficient[beta,q1[44]]

Coefficient[beta,q1[46]]

Coefficient[beta,q1[35]]

solq1

Expand[tr[sigma3][[{5,6}]]/.solq1]

Comparaison des variantes BD et BN

Les deux choix :

subMBN={M122->-M221/16,M121->-M112/88};
subMDB={M122->M111-M221/8,M121->M222-M112/44};

donnent les contraintes de cisaillement:

cis1=CE[((tr[sigma3][[{5,6}]]/.solq1)/.subMBN)/.solM]
cis2=CE[((tr[sigma3][[{5,6}]]/.solq1)/.subMDB)/.solM]

et les efforts tranchants

t1=integ[tr[cis1]]
t2=integ[tr[cis2]]

qu'on égale à {1,1}

sol1=Solve[t1=={1,1},{q1[34],q1[46]}][[1]]
sol2=Solve[t2=={1,1},{q1[34],q1[46]}][[1]]

Les contraintes de cisaillement deviennent alors :

cisBN=CE[cis1/.sol1]
cisDB=CE[cis2/.sol2]

dim[cisBN]

Fonctions graphiques

Off[Export::"obsalt"];

grcisBN[cis_,couleur_]:=
Block[{cis1,col},
col=Which[couleur==1,RGBColor[1,0,0],couleur==2,RGBColor[0,1,0],couleur==3,RGBColor[0,0,1],couleur==4,RGBColor[1,1,0],couleur==5,RGBColor[1,1,0],couleur==6,RGBColor[0,1,1],True,RGBColor[0,0,0]];
cis1=Which[z<-0.25,cis[[1]],z<0.25,cis[[2]],True,cis[[3]]];
ParametricPlot[{cis1,z},{z,-0.5,0.5},Ticks->{Automatic,{-0.5,-0.25,0.25,0.5}},PlotStyle->{col}(*,DisplayFunction->Identity*)]];

grcisDB[cis_,couleur_]:=
Block[{cis1,col},
col=Which[couleur==1,RGBColor[1,0,0],couleur==2,RGBColor[0,1,0],couleur==3,RGBColor[0,0,1],couleur==4,RGBColor[1,1,0],couleur==5,RGBColor[1,1,0],couleur==6,RGBColor[0,1,1],True,RGBColor[0,0,0]];
cis1=Which[z<-0.25,cis[[1]],z<0.25,cis[[2]],True,cis[[3]]];
ParametricPlot[{cis1,z},{z,-0.5,0.5},Ticks->{Automatic,{-0.5,-0.25,0.25,0.5}},PlotStyle->{col,Dashing[{0.008,0.008}]}(*,DisplayFunction->Identity*)]];

lignes=Graphics[{RGBColor[0,0,0],Dashing[{0.1,0.05}],Line[{{0,0.25},{1.5,0.25}}],Line[{{0,-0.25},{1.5,-0.25}}]}];
leg=Graphics[{RGBColor[1,0,0],Text[StyleForm[Subscript[τ,2],FontFamily->"Times",FontSize->20],{2,0.18}],
RGBColor[0,0,1],Text[StyleForm[Subscript[τ,1],FontFamily->"Times",FontSize->20],{1,0.35}]}];

Show[{grcisBN[cisBN[[1]],3],grcisBN[cisBN[[2]],1],grcisDB[cisDB[[1]],3],grcisDB[cisDB[[2]],1],lignes,leg},Ticks->{Automatic,{-0.5,-0.25,0.25,0.5}},AxesLabel->{StyleForm["τ",FontFamily->"Times",FontSize->20],StyleForm["z",FontFamily->"Times",FontSize->20]},PlotRange->All,AspectRatio->0.7]

Profil des déplacements complémentaires

subcis={q1[36]->0,q1[45]->0,q1[35]->0,q1[44]->0};

ucBN =CE[(tr[tr[uposol3][[{1,2}]]/.sol1])/.subcis];
ucDB=CE[(tr[tr[uposol3][[{1,2}]]/.sol2])/.subcis];

dim[ucBN]

tr[ucBN]

lignesuc=Graphics[{RGBColor[0,0,0],Dashing[{0.1,0.05}],Line[{{0,0.25},{0.36,0.25}}],Line[{{0,-0.25},{-0.36,-0.25}}]}];
leguc=Graphics[{RGBColor[1,0,0],Text[StyleForm[Subscript["u","2"],FontFamily->"Times",FontSize->20],{0.5,0.45}],
RGBColor[0,0,1],Text[StyleForm[Subscript["u","1"],FontFamily->"Times",FontSize->20],{0.5,0.30}],RGBColor[0,0,0]}];

Show[{grcisBN[tr[ucBN][[1]],3],grcisBN[tr[ucBN][[2]],1],grcisDB[tr[ucDB][[1]],3],grcisDB[tr[ucDB][[2]],1],lignesuc,leguc},Ticks->{Automatic,{-0.5,-0.25,0.25,0.5}},AxesLabel->{StyleForm["u",FontFamily->"Times",FontSize->20],StyleForm["z",FontFamily->"Times",FontSize->20]},PlotRange->All,AspectRatio->0.7]

Mise en évidence des degrés concernés par une expression xx

Clear[tabqdif];
tabqdif[exp_,xx_]:=Module[{lis,aa,bb,cc},
lis={};bb=Flatten[xx];cc=Flatten[exp];
Do[aa=Coefficient[bb[[i]],q1[k]];If[Not[aa===0]&&Not[aa===0.]&&Not[MemberQ[lis,q1[k]]],lis=Append[lis,q1[k]]],
{k,3dixi1+1,5dixi1},{i,dim1[bb]}];Print[lis];
Table[Coefficient[cc,lis[[k]]],{k,dim1[lis]}]];

depge=Flatten[Join[vq,wq,beta]];

tabqdif[depge,tr[tr[uposol3][[{1,2}]]]]

Contraintes de membrane et moments de flexion-torsion

On compare les résultats du déplacement complémentaire à ceux de la théorie classique

Contraintes principales

inst0;sigmapr=CE[Table[If[OddQ[n],inst1,inst2];(sigpr/.substidef)/.subtot,{n,nc}]];inst0;

`Différence entre contraintes principales et contraintes calculées

difsig=CE[sigma3-sigmapr];

Mise en évidence des degrés concernés par cette différence

Clear[tabqdif];
tabqdif[exp_,xx_]:=Module[{lis,aa,bb,cc},
lis={};bb=Flatten[xx];cc=Flatten[exp];
Do[aa=Coefficient[bb[[i]],q1[k]];If[Not[aa===0]&&Not[aa===0.]&&Not[MemberQ[lis,q1[k]]],lis=Append[lis,q1[k]]],
{k,3dixi1+1,5dixi1},{i,dim1[bb]}];Print[lis];
Table[Coefficient[cc,lis[[k]]],{k,dim1[lis]}]];

depge=Flatten[Join[vq,wq,beta]];

depge

Différence entre les contraintes de membrane

tabqdif[depge,tr[difsig][[{1,2,4}]]]

Différence entre les contraintes de cisaillement tranverse

tabqdif[depge,tr[difsig][[{5,6}]]]

Différence sur sigma33

tabqdif[depge,tr[difsig][[3]]]

tr[sigmapr][[3]]

Moments des contraintes principales (théorie classique)

mompr=integ[tr[tr[z sigmapr/.subtot][[{1,2,4}]]]]

Moments fournis par la théorie BN

moment3/.subtot

tabqdif[depge,(moment3/.subtot)-mompr]