{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Arches par les équations d'équilibre local (CORRIGE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Données et paramétrisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Données géométriques fixes\n",
    "# Rayon R, épaisseur h, profondeur b, section, section corrigée (eff tranchant), \n",
    "# moments quadratiques, angles origine et extrémité\n",
    "var('R h b S S_1 I J theta_A theta_B alpha'); \n",
    "theta_A = 0; theta_B = alpha; show('Arche: angle origine=', theta_A,'   angle extremite =', theta_B)\n",
    "# Variables géométriques (abscisse curviligne)\n",
    "var('theta eta')\n",
    "assume(theta>0); assume(alpha>0); assume(eta>0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Vecteur position OG (d'abscisse curviligne R*theta)\n",
    "X_s = vector(SR,3,[R*cos(theta),R*sin(theta),0]); show('Vecteur position du point G, X_s =',X_s);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Base cylindrique en G (d'abscisse curviligne R*theta)\n",
    "e_r = vector(SR,3,[cos(theta),sin(theta),0]); e_t = vector(SR,3,[-sin(theta),cos(theta),0]);\n",
    "e_z = vector(SR,3,[0,0,1]);\n",
    "show('e_r =',e_r,'  e_t =', e_t,'  e_z =',e_z); "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Variables sthéniques: efforts à l'extrémité B\n",
    "var('E F_x F_y F_z M_x M_y M_z'); #show(F_x, F_y, F_z, M_x, M_y, M_z)\n",
    "F_B = vector(SR,3,[F_x, F_y, F_z]); show('F_B =', F_B)\n",
    "M_B = vector(SR,3,[M_x, M_y, M_z]); show('M_B =', M_B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Spécialisation problème plan de l'arche avec force verticale en B ou effort vertical réparti\n",
    "var('T P p')\n",
    "p = 0; # effort réparti nul\n",
    "#P = 0; # effort terminal nul\n",
    "F_B = F_B.subs(F_x=T,F_y=-P,F_z=0); show('F_B =', F_B)\n",
    "M_B = M_B.subs(M_x=0,M_y=0,M_z=0); show('M_B =', M_B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Valeurs numériques (pour affichage)\n",
    "# Rayon RR, épaisseur hh, angle en B\n",
    "RR = 1; hh = RR/5; aalpha = pi/2;# aalpha = 7*pi/16;\n",
    "# Surface section, section modifiée, moment quadratique, polaire\n",
    "SS = hh^2; SS_1 = SS; II = hh^4/12; JJ = 2*II; "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Vecteur position initiale\n",
    "X = X_s.subs(R=RR); show('X = ',X)\n",
    "# Ligne neutre initiale\n",
    "Init = parametric_plot((X[0],X[1]), (theta, 0, aalpha),color='black',linestyle='--'); show(Init)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Géométrie 2D reconstruite\n",
    "var('r')\n",
    "# Sections initiale et terminale\n",
    "X_A = X.subs(theta=theta_A) + r*e_r.subs(theta=theta_A); #show('X_A =',X_A)\n",
    "S_A = parametric_plot((X_A[0],X_A[1]), (r, -hh/2, +hh/2),color='black',linestyle='--',thickness=2); \n",
    "X_B = X.subs(theta=aalpha) + r*e_r.subs(theta=aalpha); #show('X_B =',X_B)\n",
    "S_B = parametric_plot((X_B[0],X_B[1]), (r, -hh/2, +hh/2),color='black',linestyle='--',thickness=2);\n",
    "#show(S_A+S_B)\n",
    "# Bord supérieur initial\n",
    "XSup = X + h/2*e_r; # show('XSup = ',XSup)\n",
    "XSup = XSup.subs(R=RR,h=hh); #show('XSup = ',XSup)\n",
    "# Bord inférieur initial\n",
    "XInf = X - h/2*e_r; # show('XInf = ',XInf)\n",
    "XInf = XInf.subs(R=RR,h=hh); #show('XInf = ',XInf)\n",
    "# Visualisation\n",
    "InitSup = parametric_plot((XSup[0],XSup[1]), (theta, 0, aalpha),color='black',linestyle='--'); #show(Init)\n",
    "InitInf = parametric_plot((XInf[0],XInf[1]), (theta, 0, aalpha),color='black',linestyle='--');\n",
    "INIT = Init+InitSup+InitInf+S_A+S_B; show(INIT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Torseur des efforts intérieurs (par intégration des équations d'équilibre)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Equation d'équilibre local en résultante\n",
    "function('R_x')(theta); function('R_y')(theta); function('M_f')(theta);\n",
    "# Résultante des efforts intérieurs R_G\n",
    "R_G = vector(SR,3,[R_x(theta),R_y(theta),0]); #Eq_R = derivative(R_G,theta) == vector(SR,3,[0,0,0]); show(Eq_R)\n",
    "Eq_x = derivative(R_x(theta),theta) == 0; show('Eq equi en x: ', Eq_x)\n",
    "Eq_y = derivative(R_y(theta),theta) - R*p == 0; show('Eq equi en y: ', Eq_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Résolution des équations en résultante\n",
    "Sol_x = desolve(Eq_x,R_x(theta),[alpha,F_B[0]],ivar=theta); R_G[0] = Sol_x;\n",
    "Sol_y = desolve(Eq_y,R_y(theta),[alpha,F_B[1] ],ivar=theta); R_G[1] = Sol_y;\n",
    "show('R_G = ', R_G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Equilibre en moment\n",
    "function('M_x')(theta); function('M_y')(theta); function('M_z')(theta);\n",
    "M_G = vector(SR,3,[0,0,M_z(theta)]);\n",
    "# Variation des moments\n",
    "dM = vector(SR,3,[derivative(M_x(theta),theta),derivative(M_y(theta),theta),derivative(M_z(theta),theta)]); show('dM =', dM)\n",
    "# Moment de la résultante (x R)\n",
    "RtxR = R*e_t.cross_product(R_G); show('R*e_txR_G =', RtxR)\n",
    "# Equation sur le moment fléchissant\n",
    "Eq_z = dM[2] + RtxR[2] == 0; show('Eq equi en z: ', Eq_z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Résolution de l'équation en moment\n",
    "Sol_z = desolve(Eq_z,M_z(theta),[alpha,M_B[2]],ivar=theta); M_G[2] = Sol_z; show('Moments résultants, M_G =', M_G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Efforts normal, tranchant, moment fléchissant en G (d'abscisse curviligne R*theta)\n",
    "N_G = R_G*e_t; show('N_G = ',N_G)\n",
    "T_G = R_G*e_r; show('T_G = ',T_G)\n",
    "Mf_G = M_G*e_z; Mf_G = Mf_G.simplify_full() ; show('Mf_G =',Mf_G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Formules de Bresse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Vecteur position en L (d'abscisse curviligne R*eta)\n",
    "X_l = X_s.subs(theta=eta); show('X_l =',X_l);\n",
    "# Base cylindrique en L (d'abscisse curviligne R*eta)\n",
    "el_r = e_r.subs(theta=eta); el_t = e_t.subs(theta=eta);\n",
    "show('el_r =',el_r,'  el_t =', el_t,'  e_z =',e_z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Efforts généralisés en L (d'abscisse curviligne R*eta)\n",
    "R_L = R_G.subs(theta=eta); show('R_L = ',R_L)\n",
    "M_L = M_G.subs(theta=eta); show('M_L = ',M_L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Efforts normal, tranchant, moment fléchissant en L\n",
    "N_L = N_G.subs(theta=eta); show('N_L = ',N_L)\n",
    "T_L = T_G.subs(theta=eta); show('T_L = ',T_L)\n",
    "Mf_L = Mf_G.subs(theta=eta); show('Mf_L =',Mf_L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Déplacement en G du aux efforts normal et tranchant\n",
    "u_n = R*integrate(N_L/E/S*el_t + T_L/E/S_1*el_r,eta,0,theta); \n",
    "show('u_n =', u_n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Déplacement en G du à la flexion\n",
    "omega_L = 1/E/I*M_L; show('omega_L = ', omega_L); # Rotation relative en L\n",
    "LG = X_s - X_l; show('LG = ',LG)\n",
    "u_f = R*integrate(omega_L.cross_product(LG),eta,0,theta); \n",
    "u_f = u_f.simplify_full(); \n",
    "show('u_f =', u_f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rotation de la section en G\n",
    "Omega = R*integrate(omega_L,eta,0,theta); show('Omega =', Omega)\n",
    "Omega = Omega.simplify_full(); show('Omega =', Omega)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Traitement de la condition de symétrie"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Condition de symétrie en B (avec la seule déformation de flexion)\n",
    "Sym = u_f[0].subs(theta=alpha) == 0; #show('Condition de symétrie: ',Sym)\n",
    "Sym = Sym.simplify_full(); show('Condition de symétrie: ',Sym)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Détermination de T en fonction de P due à la condition de symétrie\n",
    "if P==0: P = alpha*p*R; # actif pour la charge répartie seule (charge intégrée)\n",
    "Sol = solve(Sym,T); show(Sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extraction de T solution\n",
    "T_sol = Sol[0]; show(T_sol)\n",
    "TsP = factor(T_sol.rhs(),P)/P\n",
    "show('T/P =',TsP); "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Impression du rapport T/P pour COPIER-COLLER\n",
    "print(TsP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cas particulier d'une arche demi-circulaire\n",
    "show('T/P pour arche demi circulaire =',TsP.subs(alpha=pi/2),' = ',n(TsP.subs(alpha=pi/2)))   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evolution du rapport T/P en fonction de alpha\n",
    "plot(TsP,alpha,0,pi/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Déplacement du à la flexion AVEC prise en compte de la condition de symétrie\n",
    "u_f = u_f.subs(T=T_sol.rhs()); u_f = u_f.simplify_full() ; show('u_f =', u_f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rotation de la section AVEC prise en compte de la condition de symétrie\n",
    "Omega = Omega.subs(T=T_sol.rhs()); Omega = Omega.simplify_full(); show('Omega =', Omega) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grandeurs numériques restantes pour affichage \n",
    "EE = 1000; GG = E/2; PP = 1/10; pp = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Position de G après déformation\n",
    "x = X + u_f; x = x.subs(R=RR,h=hh,E=EE,I=II,alpha=aalpha,P=PP,p=pp); show('x = ',x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rotation de la section \n",
    "Omegaa = Omega.subs(R=RR,h=hh,E=EE,I=II,alpha=aalpha,P=PP,p=pp); show('Omegaa =', Omegaa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ligne neutre déformée\n",
    "Def = parametric_plot((x[0],x[1]), (theta, 0, aalpha),color='black',linestyle='-'); \n",
    "show(Init+Def)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Géométrie 2D reconstruite\n",
    "# Sections initiale et terminale après déformation\n",
    "var('r')\n",
    "# Section initiale\n",
    "x_A = (x+r*e_r).subs(theta=theta_A); show('x_A =',x_A) # SANS prise en compte de la rotation\n",
    "x_A = (x_A+Omega.cross_product(r*e_r)).subs(theta=theta_A); show('x_A =',x_A);# AVEC prise en compte de la rotation\n",
    "x_A = x_A.subs(R=RR,h=hh,E=EE,I=II,P=PP,p=pp); show('x_A =',x_A); \n",
    "s_A = parametric_plot((x_A[0],x_A[1]), (r, -hh/2, +hh/2),color='black',linestyle='-',thickness=2); \n",
    "# Section terminale\n",
    "x_B = (x+r*e_r).subs(theta=aalpha); show('x_B =',x_B) # SANS prise en compte de la rotation\n",
    "x_B = (x_B+Omegaa.cross_product(r*e_r)).subs(theta=aalpha); show('x_B =',x_B); # AVEC prise en compte de la rotation\n",
    "x_B = x_B.subs(R=RR,h=hh,E=EE,I=II,P=PP,p=pp); show('x_B =',x_B)\n",
    "s_B = parametric_plot((x_B[0],x_B[1]), (r, -hh/2, +hh/2),color='black',linestyle='-',thickness=2); #show(s_A+s_B)\n",
    "# Bord supérieur déformé\n",
    "xSup = x + h/2*e_r ; #show('xSup = ',xSup); SANS prise en compte de la rotation\n",
    "xSup = xSup + Omegaa.cross_product(h/2*e_r); #show('xSup = ',xSup); # AVEC prise en compte de la rotation\n",
    "xSup = xSup.subs(R=RR,h=hh,E=EE,I=II,P=PP,p=pp); show('xSup = ',xSup)\n",
    "# Bord inférieur déformé\n",
    "xInf = x - h/2*e_r ; #show('xInf = ',xInf); # SANS prise en compte de la rotation\n",
    "xInf = xInf - Omegaa.cross_product(h/2*e_r); #show('xInf = ',xInf); # AVEC prise en compte de la rotation\n",
    "xInf = xInf.subs(R=RR,h=hh,E=EE,I=II,P=PP,p=pp); show('xInf = ',xInf)\n",
    "# Visualisation\n",
    "DefSup = parametric_plot((xSup[0],xSup[1]), (theta, 0, aalpha),color='black',linestyle='-'); #show(Init)\n",
    "DefInf = parametric_plot((xInf[0],xInf[1]), (theta, 0, aalpha),color='black',linestyle='-');\n",
    "DEF = Def+DefSup+DefInf+s_A+s_B; show(DEF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualisation\n",
    "show(INIT+DEF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.1",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}