1-修正剑桥模型预测sigma3为常数时应力应变关系

#include<stdio.h>
#include<stdlib.h> 
#include<math.h> 
#include<float.h>
#include<time.h>
//常量定义

#define N 50


//主函数
int main()
{
    FILE *fp;
    fp = fopen("剑桥模型sigema3为常数预测.txt", "w");
    //变量定义
    //if ((fp = fopen("剑桥模型sigema3为常数预测.txt", "a")) == NULL)
    //{
    //    printf("cannot open file\n");
    //    exit(0);
    //}


    double lambda = 0.0853, kappa = 0.0188, M = 1.36, e0 = 0.68, poisson = 0.3;
    double p[N] = { 0.0 };//存放平均应力p
    double q[N] = { 0.0 };//存放偏应力q
    q[0] = 0.0;
    p[0] = 196.0;
    double yita[N] = {0.0};//应力比
    yita[0] = q[0] / p[0];
    double Dpp[N] = {0.0};//柔度矩阵元素
    double Dpq[N] = { 0.0 };
    double Dqp[N] = { 0.0 };
    double Dqq[N] = { 0.0 };
    
    double depsilon1[N] = { 0.0 };
    depsilon1[0] = 0.0;
    depsilon1[1] = 0.005;
    depsilon1[2] = 0.010;
    depsilon1[3] = 0.015;
    for (int i = 4; i < N; i++)
    {
        depsilon1[i] = 0.02;

    }
    double depsilon3[N] = { 0.0 };
    double dsigma1[N] = { 0.0 };
    double sigma1[N] = { 196.0 };
    //sigma1[0] = 196.0;
    double dsigma3[N] = { 0.0 };//全部为零
    double sigma3[N] = { 0 };//
    sigma3[0] = 196.0;
    //double dsigma3 = 0.0;
    double dp[N] = { 0.0 };
    double dq[N] = { 0.0 };
    double depsilonv[N] = { 0.0 };
    double epsilonv[N] = { 0.0 };
    double depsilond[N] = { 0.0 };
    double epsilond[N] = { 0.0 };
    double ck = kappa / (1 + e0);
    double cp = (lambda - kappa) / (1 + e0);
    printf("p    q    yita    Dpp    Dpq    Dqp    Dqq    depsl1    depsl3    dsgm1    dsgm3    dp    dq    depslv    epslv    depsld    epsld    sigma1    sigma3 \n");
    fprintf(fp,"p    q    yita    Dpp    Dpq    Dqp    Dqq    depsl1    depsl3    dsgm1    dsgm3    dp    dq    depslv    epslv    depsld    epsld    sigma1    sigma3 \n");
    for (int i = 1; i < N; i++)
    {
        Dpp[i] = ck + cp*(M*M - yita[i - 1] * yita[i - 1]) / (M*M + yita[i - 1] * yita[i - 1]); //             
        Dpq[i] = cp * 2 * yita[i - 1] / (M*M + yita[i - 1] * yita[i - 1]);
        Dqp[i] = Dpq[i];
        Dqq[i] = (2.0 / 9.0) * ck*(1 + poisson) / (1 - 2 * poisson) + cp * 4 * yita[i - 1] * yita[i - 1] / (M*M*M*M - yita[i - 1] * yita[i - 1] * yita[i - 1] * yita[i - 1]);
        dsigma1[i] = 9.0 * p[i - 1] * depsilon1[i] / (Dpp[i] + 3 * Dpq[i] + 3 * Dqp[i] + 9 * Dqq[i]);
        depsilon3[i] = 0.5*depsilon1[i]*(2 * Dpp[i] + 6 * Dpq[i] - 3 * Dqp[i] - 9 * Dqq[i]) / (Dpp[i] + 3 * Dpq[i] + 3 * Dqp[i] + 9 * Dqq[i]);
        dp[i] = 1 / 3.0*(dsigma1[i] + 2 * dsigma3[i]);
        p[i] = p[i - 1] + dp[i];
        dq[i] = dsigma1[i] - dsigma3[i];
        q[i] = q[i - 1] + dq[i];
        yita[i] = q[i] / p[i];
        depsilonv[i] = depsilon1[i] + 2 * depsilon3[i];
        epsilonv[i] = epsilonv[i - 1] + depsilonv[i];
        depsilond[i] = 2 / 3.0*(depsilon1[i] - depsilon3[i]);
        epsilond[i] = epsilond[i - 1] + depsilond[i];
        sigma1[i] = sigma1[i - 1] + dsigma1[i];    
        sigma3[i] = sigma3[i-1] + dsigma3[i];
        printf("%.1lf    %.1lf    %.3lf    %.3lf    %.3lf    %.3lf    %.3lf    %.3lf    %.4lf    %.2lf    %.1lf    %.2lf    %.2lf    %.3lf    %.3lf    %.3lf    %.3lf    %.1lf    %lf\n", p[i], q[i], yita[i], Dpp[i], Dpq[i], Dqp[i], Dqq[i], depsilon1[i], depsilon3[i], dsigma1[i], dsigma3[i], dp[i], dq[i], depsilonv[i], epsilonv[i], depsilond[i], epsilond[i], sigma1[i], sigma3[i]);
        fprintf(fp,"%.1lf    %.1lf    %.3lf    %.3lf    %.3lf    %.3lf    %.3lf    %.3lf    %.4lf    %.2lf    %.1lf    %.2lf    %.2lf    %.3lf    %.3lf    %.3lf    %.3lf    %.1lf    %.0lf    \r\n", p[i], q[i], yita[i], Dpp[i], Dpq[i], Dqp[i], Dqq[i], depsilon1[i], depsilon3[i], dsigma1[i], dsigma3[i], dp[i], dq[i], depsilonv[i], epsilonv[i], depsilond[i], epsilond[i], sigma1[i], sigma3[i]);
        
    }
    fclose(fp);//关闭输出结果文档



}

h