Déplacement complémentaire dans une plaque stratifiée

degré en (x, y) = deg, degré en z = ddz
données B-D

Introduction

Il s'agit de résoudre le problème de conditions initiales en z = z0 pour l'équation
	div(K (Subscript[grad, s]) u)= -div[K (SuperDagger[B])ζ]+((Subscript[λ, 1])+z(Subscript[λ, 4]))(Subscript[(Overscript[e, ->]), 1])+((Subscript[λ, 2])+z(Subscript[λ, 5]))(Subscript[(Overscript[e, ->]), 2])+(Subscript[λ, 3])(Subscript[(Overscript[e, ->]), 3])

où :
	l'inconnue est le déplacement complémentaire u,
	les données sont ζ et les multiplicateurs de Lagrange (Subscript[λ, j]).
	
	Nota : pour faciliter la lecture on suit quelques règles simples dans le choix des noms de variables :
		le chiffre "0" est réservé aux conditions initiales,
		le chiffre "3" est réservé au vecteur contrainte {(Subscript[σ, 13]),(Subscript[σ, 23]),(Subscript[σ, 33])} qui est une fonction continue de la cote z
		le groupe "pr" aux déformations et aux contraintes principales,
		le groupe "co" aux déformations et aux contraintes complémentaire,
		le groupe "sig" à des contraintes,
		le groupe "def" à des déformations.
	
	On cherche u sous la forme d'un polynômede degré deg en (x,y) et de degré ddz en z,
	On se donne ζ et les multiplicateurs de Lagrange (Subscript[λ, j]) sous la forme de polynômes de degré deg en (x,y).
	Les "conditions initiales" sont, en z = 0, la valeur de u et des trois contraintes (Subscript[σ, 13]),(Subscript[σ, 23]),(Subscript[σ, 33]).
		
	Commentons maintenant la figure suivante. 
	Les inconnues sont les q(i), groupés en une colonne q, qui sont les coefficients des monômes dans le polynôme u cherché.Ils interviennent à la fois dans la représentation de div(K (Subscript[grad, s]) u) et dans celle des conditions initiales où ils se nomment alors q0.
	C'est pourquoi le système linéaire qu'on va construire a pour matrice la superposition des trois sous-matrices : "system"={mat, condu0, condsigco0} qui concernent respectivement div(K (Subscript[grad, s]) u), les conditions initiales en déplacement, les conditions initiales en contraintes.
	Si deg est le degré des polynômes par rapport à (x,y), ddz celui en z, alors il y a dixi= (∑_(d=0)^deg)(d+1) monômes en (x,y) de degré inférieur ou égal à deg, et dixi(ddz+1) monômes en (x,y,z) de degré inférieur ou égal à (deg+1)(ddz+1).
	Il y a donc 3 dixi(ddz+1) coefficients q(j) et donc autant de colonnes pour les matrices mat, condu0 et condsig0. Ces deux dernières comportent chacune 3 dixi lignes. La matrice "mat" n'en doit donc comporter que  
				3dixi(ddz+1)-6dixi = 3 dixi(ddz-1)
lignes indépendantes.C'est ce qu'on vérifiera au passage.


	Les coordonnées cartésiennes sont (x,y,z). Dans la première partie qui traite du problème de Cauchy pour une couche homogène z est une variable relative qui vaut 0 sur la face inférieure. La cote absolue de cette dernière est notée z0, en sorte que z+z0 est la cote absolue courante. Ce changement de variable permet de garder un caractère purement numérique à la matrice "system" qu'on devra inverser, la variable formelle z0 ne restant présente qu'au second membre.

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"

(*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}]];

Suffixe pour les fichiers enregistrés

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

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

Entrez le nombre de couches 'nc'

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}];

Entrez le degré en xy 'deg' et celui en z 'ddz

deg=2;ddz=5;dixi=Sum[1+k,{k,0,deg}];
Print[suffixe];

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}}

Fonctions diverses

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]];

Simplification d'expressions pour comparaisons

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

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]};

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]};

Constitution d'une matrice d'une matrice déduite de 'matbrute' par suppression des lignes nulles ou linéairement dépendantes des précédentes

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];

redu0[system_]:=
Module[{colut,d1,d2,la,ligut,ligut1,ligdep,mat,n1,n2,system0,system1,x},
n1=lignenulles[system];Print[n1];
n2=lignenulles[tr[system]];
{d1,d2}=dim[system];
colut=Complement[Range[1,d2],n2];
ligut=Complement[Range[1,d1],n1];
ligdep={};
system1=Chop[system[[ligut]]];
{d1,d2}=dim[system1];
mat={system1[[1]]};
Off[LinearSolve::"nosol"];
ligut1={1};
Do[
x=system1[[i]];
(*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".
Mais, cette instrution n'étant fiable qu'en numérique il faut d'abord instancier : *)
inst1;
la=LinearSolve[Transpose[mat],x];If[StringPosition[ToString[la],"LinearSolve"]!={},
(*puis désinstancier :*)
inst0;mat=Append[mat,x];ligut1=Append[ligut1,i],
inst0;ligdep=Append[ligdep,i],
inst0],
{i,2,d1}];
ligdep=Table[ligut[[ligdep[[k]]]],{k,dim1[ligdep]}];
ligut=Table[ligut[[ligut1[[k]]]],{k,dim1[ligut1]}];
(*réactivation des messages d'alerte*)
On[LinearSolve::"nosol"];
{ligdep,ligut,colut,mat}];

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}]]]];

Présentation et sélection 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.

subzeta[deg_]:=
Module[{dixi,tabs1},
dixi=Sum[k+1,{k,0,deg}];
tabs1=Flatten[Table[s1[i,k,j-k],{i,8},{j,0,deg},{k,j,0,-1}],2];
Table[s[j]->tabs1[[j]],{j,8dixi}]];
Print[subzeta[deg][[Range[1,2 dixi]]]];

Fonction inverse

subzetinv[deg_]:=
Module[{dixi,tabs1,sub},
dixi=Sum[k+1,{k,0,deg}];
tabs1=Flatten[Table[s1[i,k,j-k],{i,8},{j,0,deg},{k,j,0,-1}],2];
Table[tabs1[[j]]->s[j],{j,dim1[tabs1]}]];

Print[subzetinv[deg][[Range[1,2 dixi]]]];

Sélection, parmi la liste des s définie pour le degré deg, de ceux qui sont associés à un degré inférieur ou égal à deg1

zetdeg[deg_,deg1_]:=
Module[{dixi,dixi1},
dixi=Sum[1+k,{k,0,deg}];
dixi1=Sum[1+k,{k,0,deg1}];
Flatten[Table[s[dixi(j-1)+k],{j,8},{k,dixi1}]]];

Print[zetdeg[2,1]];

Premier membre du système différentiel
div(K (Subscript[grad, s]) u)

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];

Représentation polynomiale du déplacement complémentaire

On représente ce déplacement complémentaire par un polynôme en (x, y, z),
 de degré 'deg' en (x,y) et ddz en z, 
 dont les coefficients q sont des scalaires en nombre 3*dixi*(ddz+1)
L'ordre de numérotation des coefficients résulte de l'aplatissement de l'arborescence suivante :
n° de composante de u (1 à 3)
n° de monôme en (x,y) (1 à dixi)
n° de monôme en z (1 à ddz).
Le déplacement complémentaire est la somme des monômes q[i] monom[i].

polynomes en (x,y,z)

On explicite ci-dessous l'affectation des coefficients q(i)

monom=Flatten[Table[q[(i-1) dixi (ddz+1)+(j-1) (ddz+1)+k] mon[[j]] z^(k-1),{j,dixi},{k,ddz+1},{i,3}],1];
(*et, à titre didactique :*)
Print[monom[[Range[1,6]]]];

Représentation polynomiale du déplacement complémentaire

upo=
Table[Sum[monom[[k,j]],{k,dim1[monom]}],{j,dim2[monom]}];

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

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

Matrice du premier membre "matbrute"

Le premier membre de l'équation différentielle est la divergence "divsigco" qui dépend des 3 dixi (ddz+1)(soit 108 pour deg = 2 et ddz = 5) paramètres q[i].

	On identifie maintenant, dans divsigco, les coefficients des 3 dixi (ddz+1) monomes de la forme z^(k-1) mon[[j]] q[p] où j varie de 1 à dixi, k de 0 à ddz.
Pour faire cette identification il convient de tout multiplier par x y z pour que q[p] soit toujours multiplié par xyz au moins, sinon on oublie les puissances nulles. On va donc identifier en fait les coefficients de x y z^k mon[[j]] q[p] dans x y z divsigco.

	On affiche la durée de ce calcul formel qui prend un peu de temps.

matbrute=Table[0,{p,3 dixi (ddz+1)},{m,3 dixi (ddz+1)}];
tt=Timing[
Block[{interm},
interm=Expand[x y z divsigco];
Do[matbrute[[(i-1) dixi (ddz+1)+(j-1)(ddz+1)+k,p]]=Coefficient[interm[[i]],x y z^k mon[[j]] q[p]],{i,3},{j,dixi},{k,ddz+1},{p,3 dixi (ddz+1)}]]][[1]];
Print[tt];

Print[matbrute[[1]]];

Détermination des lignes indépendantes de matbrute

Calcul formel

Il a le mérite de montrer que les lignes nulles ou dépendantes ne dépendent pas de la raideur.

{tt,{ligdep,ligut,colut,mat}}=Timing[redu0[matbrute]];
Print["dimensions de 'matbrute' : ", dim[matbrute]];
Print["dimensions de 'mat' : ", dim[mat]];
Print[dim1[ligdep]," lignes dépendantes : ",ligdep];
Print[dim2[mat]-dim1[colut]," colonnes nulles : ",Complement[Range[1,dim2[mat]],colut]];
Print[tt];

Calcul numérique

inst1;matbrute1=matbrute;
{ligdep1,ligut1,colut1,mat1}=redu0[matbrute1];
Print[ligdep1];
inst0;

inst2;matbrute2=matbrute;
{ligdep2,ligut2,colut2,mat2}=redu0[matbrute2];
Print[ligdep2];
inst0;

Indispensable dans le cas où l'on a sauté le calcul formel

inst0;mat=matbrute[[ligut1]];

Etude des lignes inutiles de matbrute

numéro de composante de la divergence, 
degrés en x, y, z du monôme d'identification

'expo' est la table des lignes du type (n,expox,expoy,expoz};
Chacune d'elles correspond à une ligne de 'matbrute', c'est-à-dire à une monôme.
'n' est le numéro de la composante de la divergence, expox l'exposant de x etc.

expo=Flatten[Table[Flatten[{i,exposant[deg][[j]],k-1}],{i,3},{j,dixi},{k,ddz+1}],2];

lecturexpo affiche les lignes précédentes, commentées, pour chacune des lignes de la liste 'ex'

lecturexpo[ex_]:=
Block[{lec,d1,x,y,z},
d1=dim1[ex];
lec=Table[{"composante",ex[[i,1]],"  monôme :",x^ex[[i,2]]y^ex[[i,3]]z^ex[[i,4]]},{i,d1}];
Print[lec]];

Lignes nulles de matbrute

Les monômes concernés par matbrute sont de degré inférieur ou égal à deg+ddz-2 puisque on a dérivé deux fois : une fois pour les déformations, une fois pour la divergence. Les lignes correspondantes sont donc nulles.

lignul=Complement[Range[1,dim1[matbrute]],Join[ligut,ligdep]];

lecturexpo[expo[[lignul]]];

Lignes de matbrute non nulles mais combinaisons linéaires des précédentes

inst0;
Print[ligdep];

lecturexpo[expo[[ligdep]]];

Second membre
-div[K (SuperDagger[B])ζ+(Subscript[ε, o])]+((Subscript[λ, 1])+z(Subscript[λ, 4]))(Subscript[(Overscript[e, ->]), 1])+((Subscript[λ, 2])+z(Subscript[λ, 5]))(Subscript[(Overscript[e, ->]), 2])+(Subscript[λ, 3])(Subscript[(Overscript[e, ->]), 3])

Les différents éléments du système à contruire (rappel)

Multiplicateurs de Lagrange, polynômes en x et y, de coefficients r[j]

5 multiplicateurs mis sous forme polynomiale

lag=Table[r[dixi(i-1)+j],{i,5},{j,dixi}].mon;

Contribution des multiplicateurs au second membre

rcol=Table[r[i],{i,5 dixi}];
scol=Table[s[i],{i,8 dixi}];

secmulti={lag[[1]]+(z+z0) lag[[4]],lag[[2]]+(z+z0) lag[[5]],lag[[3]]};
Print[secmulti];

Déformations et contraintes engendrées par la déformation généralisée ζ ;
div[K( (SuperDagger[B])ζ+(Subscript[ε, o]))] avec (Subscript[ε, o]) = 0

Les s[i], pour i de 1 à 8 dixi, sont les coefficients des déformations généralisées zeta mises sous forme polynomiales ; on construit d'abord "zetamat" qui montre les monômes dont la somme est zeta.
membrane :	i=1,2,3 ;
flexion-torsion : i=4,5,6 ;
cisaillement : i=7,8 ;
defpr : déformations principales (SuperDagger[B]) ζ de la théorie ;
sigpr : contraintes tridimensionnelles principales ;
divsigpr : divergence des contraintes tridimensionnelles principales

Off[General::"spell1"];
zetamat=Table[s[(i-1) dixi+j] mon[[j]],{i,8},{j,dixi}];
Print[zetamat];
zeta=Table[Sum[zetamat[[i,j]],{j,dixi}],{i,8}];
defpr={zeta[[1]]+(z+z0) zeta[[4]],zeta[[2]]+(z+z0) zeta[[5]],0,zeta[[3]]+(z+z0) zeta[[6]],zeta[[7]],zeta[[8]]};
sigpr=k3d.defpr;
divsigpr=Expand[div[sigpr]];

Second membre "secmem1" à 18 dixi lignes, associé au premier membre "matnou"

3 lignes polynomiales

inst0;
secsca=Expand[secmulti-divsigpr];

dim[secsca]

une ligne par coefficient

inst0;secmem=Flatten[Table[Coefficient[Expand[x y z secsca[[j]]],x y z^k  mon[[i]]],{j,3},{i,dixi},{k,ddz+1}]];
Print[dim[secmem]];

Dépendance des équations avec second membre

On constate que la structure de secmem est compatible avec celle de 'matbrute'.

secmem[[lignul]]

Second membre réduit à la dimension du premier

secmem1=secmem[[ligut]];
Print[dim1[secmem1]];inst0;

Solution du système avec conditions initiales

4.1 Conditions initiales

On prend les 72 équations précédentes, et on rajoute 6 dixi conditions en z = 0 : 3 dixi pour le déplacement complémentaire et 3 dixi pour les contraintes (Subscript[σ, 13])

72 premières lignes ; organisation du second membre sous la forme
secrr . r + secss . s

secrr=Join[Table[Coefficient[secmem1[[i]],r[j]],{i,dim1[secmem1]},{j,5 dixi}],Table[0,{i,6 dixi},{j,5 dixi}]];
secss=Join[Table[Coefficient[secmem1[[i]],s[j]],{i,dim1[secmem1]},{j,8 dixi}],Table[0,{i,6 dixi},{j,8 dixi}]];

U donné pour z = 0

Premier membre

condu0 = Flatten[Table[Coefficient[Expand[x y upo[[i]]]/.z->0, x y mon[[j]]  q[k]],
		 {i, 3}, {j, dixi}, {k, 3 (ddz+1) dixi}],1];

Second membre 
La colonne [secu0].{qcol0} est le déplacement imposé en z=0.

secu0 = Join[Table[0, {i, dim1[secmem1]}, {j, 3 dixi}],
		IdentityMatrix[3 dixi],
			Table[0, {i, 3 dixi}, {j, 3 dixi}]]; 
qcol0 = Table[q0[i], {i, 3 dixi}];

Vecteur contrainte imposée sur la face inférieure, à fournir par le seul déplacement complémentaire

Premier membre : la colonne [sigco03].{q} est le vecteur des (Subscript[σ, j3]) produites par le déplacement complémentaire, qui est représenté par les q(k)

inst0; 
sigco03 = sigco[[{3,5,6}]]/.z->0;
condsigco0 = Flatten[Table[Coefficient[Expand[x y sigco03[[i]]], x y mon[[j]] q[k]] , 
		{i,3}, {j,  dixi}, {k, 3 (ddz+1) dixi}],1];

La colonne [secsig0].{sig0} est le vecteur des (Subscript[σ, j3]) imposées en z = 0 ;
on en retranche la partie "principale" (Subscript[(K(SuperDagger[B])ζ), z=0]) de la contrainte totale imposée et on incorpore à "secss" le terme correspondant "secsigpr03"

sigcol0 = Table[sig0[i], {i, 3 dixi}];
secsig0 = Join[Table[0, {i, dim1[secmem1]+3 dixi}, {j, 3 dixi}], 
		IdentityMatrix[3 dixi]];
sigpr03=sigpr[[{3,5,6}]]/.z->0;
(* secsigpr03 = -K (SuperDagger[B])ζ *)
secsigpr03 = -Flatten[Table[Coefficient[Expand[x y sigpr03[[i]]], x y mon[[j]] s[k]] , 
		{i,3}, {j,dixi}, {k,8dixi}],1];
secss[[Range[dim1[secss]-3dixi+1,dim1[secss]]]]=secsigpr03;

4.2 Système complet et résolution

inst0;system=Join[mat,condu0,condsigco0];

inst1;
inv=Inverse[system];
solr1=inv.secrr;
sols1=inv.secss;
solq1=inv.secu0;
solsig1=inv.secsig0;
solution1=Chop[solr1.rcol+sols1.scol+solq1.qcol0+solsig1.sigcol0];

inst2;
inv=Inverse[system];
solr2=inv.secrr;
sols2=inv.secss;
solq2=inv.secu0;
solsig2=inv.secsig0;
solution2=Chop[solr2.rcol+sols2.scol+solq2.qcol0+solsig2.sigcol0];

inst0;

Assemblage des couches

Block[{utemp,soltemp,interm,cc,upo,z0,qcol,qcol1,qcol2,sigcol2},
inst0;

upo=
Table[Sum[monom[[k,j]],{k,dim1[monom]}],{j,dim2[monom]}];uposol=Table[0,{i,nc}];
qcol=Table[q[i],{i,3 (deg+1)dixi}];
qcol2=Table[q2[i],{i,3 dixi}];
qcol1=Table[q1[i],{i,3 dixi}];
qcol2=qcol1;
sigcol2=Table[0.,{i,3 dixi}];
Do[
Print["couche ",n];
z0=zc[[n]];
If[OddQ[n],
inst1;{solr,sols,solq,solsig}={solr1,sols1,solq1,solsig1},
inst2;{solr,sols,solq,solsig}={solr2,sols2,solq2,solsig2}];
soltemp=solr.rcol+sols.scol+solq.qcol2+solsig.sigcol2;

utemp=upo/.Table[q[i]->soltemp[[i]],{i,3dixi(ddz+1)}];
utemp=utemp/.sub0;

usup=Chop[Simplify[(utemp/.z->(zc[[n+1]]-z0))]];
qcol2=Flatten[Table[Coefficient[Expand[x y usup[[i]]],x y mon[[j]]],{i,3},{j,dixi}]];
sigsup=Chop[((k3d.(grads[utemp]+defpr))[[{3,5,6}]])/.{z->zc[[n+1]]-z0}];
sigcol2=Flatten[Table[Coefficient[Expand[x y sigsup[[i]]],x y mon[[j]]],{i,3},{j,dixi}]];

uposol[[n]]=Chop[utemp/.z->(z-z0)];,
{n,3}];
inst0];

Vérifications de continuité

Déplacement complémentaire

Do[	Print[Chop[Expand[(uposol[[n+1]]-uposol[[n]])/.z->zc[[n+1]]]]], {n, nc-1}];

Déformations principales en variable z absolue

defprabs=defpr/.{z->z-z0};

Continuité des contraintes (Subscript[σ, i3])

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

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

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

secmulti1=secmulti/.z->z-z0;
Do[
equ=div[sigma[[n]]]-secmulti1;
Print[Chop[Expand[equ]]],
{n,nc}];

Détermination des multiplicateurs de Lagrange
et des déplacements

Liaisons

liaisons = Block[{integ1, integ2, integ3, integ4, integ5, u},
	u = Table[uposol[[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 sur la face supérieure

sigsup=Block[{interm,cc},
If[OddQ[nc],inst1,inst2];
interm=Chop[Expand[(k3d.(grads[uposol[[nc]]]+defprabs))[[{3,5,6}]]]/.z->zc[[nc+1]]];
cc=Table[Coefficient[Expand[x y interm[[i]]],x y mon[[j]]],{i,3},{j,dixi}];
cc=Chop[Flatten[cc]/.sub0];
inst0;cc];

Equations

Mise en vecteur unique "v" des inconnues

inst0;
substi = Join[Table[q1[i]->v[i], {i, 3 dixi}], Table[r[i]->v[i + 3 dixi], {i, 5 dixi}]];

Equations de liaisons, premier et second membre

Block[{interm},
interm = Simplify[Expand[liaisons/.substi]]; 
equali = Flatten[Table[Coefficient[Expand[x y interm[[i]]], x y v[k] mon[[j]]] , {i, 5}, {j, dixi}, {k, 8 dixi}],1];
secli = -Flatten[Table[Coefficient[Expand[x y interm[[i]]], x y s[k] mon[[j]]] , {i, 5}, {j, dixi}, {k, 8 dixi}],1]];

Equations d'équilibre, premier et second membre

interm1 = Expand[Table[sigsup[[i]]/.substi, {i, 3 dixi}]];
equacon =  Table[Coefficient[interm1[[i]], v[k]] , {i, 3 dixi}, {k, 8 dixi}]/.sub0;
seccon =  -Table[Coefficient[interm1[[i]], s[k]] , {i, 3 dixi}, {k, 8 dixi}]/.sub0;

equa=Join[equali,equacon];

secmem2 = Join[secli, seccon];

inv = Inverse[equa];
sol = (inv.secmem2).scol;

(*Put[sol,"solcomp"<>suffixe];*)

Vérification de la solution

substisol=Join[Table[q1[i]->sol[[i]],{i,3 dixi}],Table[r[i]->sol[[i+3 dixi]],{i,5 dixi}]];

uposol1=Chop[Expand[uposol/.substisol]];

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

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

Print[Chop[Expand[sigma1[[1,{3,5,6}]]/.z->zc[[1]]]]];

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

Print[Chop[Expand[sigma1[[nc,{3,5,6}]]/.z->zc[[nc+1]]]]];

Équilibre intérieur

secmulti1=(secmulti/.z->z-z0)/.substisol;
Do[
equ=div[sigma1[[n]]]-secmulti1;
Print[Chop[Expand[equ]]],
{n,nc}];

Enregistrement et vérification de celui-ci :

Put[{uposol1,secmulti1},"ucomp"<>suffixe];

{uposol2,secmulti2}=Get["ucomp"<>suffixe];

uposol2==uposol1

secmulti2==secmulti1