Prog c

Fermé
helmi - 5 déc. 2006 à 03:15
Char Snipeur
Messages postés
9688
Date d'inscription
vendredi 23 avril 2004
Statut
Contributeur
Dernière intervention
2 octobre 2020
- 5 déc. 2006 à 13:05
bonjour:
voila audessous le programme qui prend le dynamique d'un reseau d'atome au cours de temps(nouvelles positions,forces...).mais moi je voulais prendre l'exemple d'une poutre = resau d'atomes et encastrée=ecouches où les atomes sont fixes ainsi comment introduire une force appliquee sur l' extrimité.
merci
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define nat 121
#define dt 0.005

// CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
// DECLARATIONS
//

FILE *pf1,*pf2,*pf3,pf4;
double rx[nat],ry[nat],vx[nat],vy[nat];
double ax[nat],ay[nat],axa[nat],aya[nat];//accel avant et apres
double fx[nat],fy[nat];
double dx,dy,r2,ep,ec,alpha,fij,t,rxi,ryi,etot;
int i,j,npas,ifixtemp;
double v2;

// CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
// CONSTANTES ET PARAMETRES
//

double c1=2.0*dt*dt/3.0,c2=-dt*dt/6.0,c3=dt/3.0;
double c4=5.0*dt/6.0,c5=-dt/6.0;

double ap=4.0,bp=4.0,af=24.0,bf=48.0;


int ifixtemp=1;
double temp=0.1;
int maxpas=10000;

main()
{

// CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
// LECTURE FICHIER CONFIGURATION INITIALE
// ET ECRITURE VISU.IN
//
pf1=fopen("CONF.IN","r");
pf2=fopen("VISU.IN","w");
for (i=0;i<nat;i++)
{
fscanf(pf1,"%lf %lf \n",&rx[i],&ry[i]);
fscanf(pf1,"%lf %lf \n",&vx[i],&vy[i]);
fscanf(pf1,"%lf %lf \n",&ax[i],&ay[i]);
fscanf(pf1,"%lf %lf \n",&axa[i],&aya[i]);
fprintf(pf2,"%lf %lf %d \n",rx[i],ry[i],i);
}
close(pf1);
close(pf2);



// CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
// BOUCLE PRINCIPALE
//
pf1=fopen("ENE","w");
pf2=fopen("TRAJ","w");
pf3=fopen("VIT","w");
for (npas=1;npas<=maxpas;npas++)
{


//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// calcul des nlles positions
//
for (i=0;i<nat;i++){
rx[i] = rx[i] + dt*vx[i] + c1*ax[i] + c2*axa[i];
ry[i] = ry[i] + dt*vy[i] + c1*ay[i] + c2*aya[i];
fx[i]=0.0;
fy[i]=0.0;
}


//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// calcul des nouvelles forces
//
ep=0.0;
for (i=0;i<nat-1;i++){
for (j=i+1;j<nat;j++){
dx=rx[i]-rx[j];
dy=ry[i]-ry[j];
r2=dx*dx+dy*dy;
fij=bf/(r2*r2*r2*r2*r2*r2*r2) - af/(r2*r2*r2*r2);
fx[i]+=fij*dx;
fy[i]+=fij*dy;
fx[j]+=-fij*dx;
fy[j]+=-fij*dy;
ep+= bp/(r2*r2*r2*r2*r2*r2) -ap/(r2*r2*r2);
}
}
ep=ep/nat;


//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// calcul des nouvelles vitesses
// + corection pout fixer T
//

ec=0.0;
for (i=0;i<nat;i++){
vx[i]=vx[i]+c3*fx[i]+c4*ax[i]+c5*axa[i];
vy[i]=vy[i]+c3*fy[i]+c4*ay[i]+c5*aya[i];
ec+=vx[i]*vx[i]+vy[i]*vy[i];
}
ec=ec*0.5/nat;
if (ifixtemp==1) {
alpha=sqrt(temp/ec);
for (i=0;i<nat;i++){
vx[i]=alpha*vx[i];
vy[i]=alpha*vy[i];

}
}

//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// actualisation des tableaux
//
for (i=0;i<nat;i++){
axa[i]=ax[i];
aya[i]=ay[i];
ax[i]=fx[i];
ay[i]=fy[i];
}



//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// sorties
//
t=(float)npas*dt;
if (fmod(npas,1000)==0) {
printf("%d / %d \n",npas,maxpas);
}

if (fmod(npas,10)==0) {
fprintf(pf1,"%d %lf %lf \n",npas,ep,ep+ec);
}

if (fmod(npas,500)==0) {
for (i=0;i<nat;i++) { fprintf(pf2,"%lf %lf \n",rx[i],ry[i]);
}
}
if (fmod(npas,100)==0) {
for (i=0;i<nat;i++) { fprintf(pf3," %lf \n",sqrt(vx[i]*vx[i]+vy[i]*vy[i]));
}
}


// FIN BOUCLE PRINCIPALE
}
close(pf1);
close(pf2);
close(pf3);

// CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
// ECRITURE FICHIER CONFIGURATION FINALE
// ET ECRITURE VISU.OUT
//
//
pf1=fopen("CONF.OUT","w");
pf2=fopen("VISU.OUT","w");
for (i=0;i<nat;i++)
{
fprintf(pf1,"%lf %lf \n",rx[i],ry[i]);
fprintf(pf1,"%lf %lf \n",vx[i],vy[i]);
fprintf(pf1,"%lf %lf \n",ax[i],ay[i]);
fprintf(pf1,"%lf %lf \n",axa[i],aya[i]);
fprintf(pf2,"%lg %lg %d \n",rx[i],ry[i],i);
}
close(pf1);
close(pf2);


// CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
// FIN DU MAIN
}

2 réponses

mamiemando
Messages postés
31292
Date d'inscription
jeudi 12 mai 2005
Statut
Modérateur
Dernière intervention
8 août 2022
7 390
5 déc. 2006 à 09:34
Heu tu sais moi je ne suis pas physicienne donc si tu ne me dis pas précisemment ce que tu veux faire je ne vois pas trop comment t'aider.

Bonne chance
0
Char Snipeur
Messages postés
9688
Date d'inscription
vendredi 23 avril 2004
Statut
Contributeur
Dernière intervention
2 octobre 2020
1 329
5 déc. 2006 à 13:05
Salut.
Pour l'encastrage, c'est simple, il suffit de ne pas modifié une partie des atomes, comme ça ils resterons fixe, ce qui simulera ton encastrage.
Par contre, pour la force, c'est plus difficile. Mais je pense que rajouter un terme dans
//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// calcul des nouvelles forces
//
devrai faire l'affaire.
en fait, ton modèle c'est un chapelet d'atomes, donc un truc linéaire?
à ce moment là, il faut encastré 2 atomes, sinon ta tige va betement pivoté.
Bonne chance
0