Let's say I have some experimental data (x,y) and I want to find some fitting parameters. But in my mathematical model, the fitting parameter are not in a simple equation that relates x and y. For finding y we should solve some sets of algebraic and differential equations and fitting parameters appears in one of these equations. So currently I am solving these equations using an incremental-iterative numerical method to obtain y. how can I do a curve fit while all the methods in Mathematica want a explicit format of a function.
Edit:
The problem I am trying to solve is the growth of damage in a material under cyclic loading. The equations that should be solved are following:
where bold variables are 3*3 tensors. These equations show a plasticity model coupled with damage. The numerical method used to solve this initial value problem is called elastic predictor-plastic corrector. The curve fitting parameters that I am trying to find are s and S in the last equation. The discretized from of these equations is as follows:
First some material properties are introduced:
Remove["Global`*"]
(* Material Properties *)
EYoung = 10000;
ν = 0.3;
smat = 8.67207899548;
S = 0.0442223584949;
σmeso = 42.4521112168564;
Dc = 0.5;
β = 0.476190476190476;
G = EYoung/(2 (1 + ν));
Cy = 100;
σy = 30;
σf = 10;
h = 0.2;
Two fourth order tensors that show the stiffness and compliance of material are formed:
DElastic =
Table[ EYoung/(
2 (1 + ν))* (KroneckerDelta[i, l]*KroneckerDelta[j, k] +
KroneckerDelta[i, k]*KroneckerDelta[j, l] ) + (
EYoung*ν)/((1 + ν) (1 - 2 ν)) KroneckerDelta[i, j]*
KroneckerDelta[k, l], {i, 1, 3}, {j, 1, 3}, {k, 1, 3}, {l, 1, 3}];
CElastic =
Table[(1 + ν)/(
2 EYoung) (KroneckerDelta[i, l]*KroneckerDelta[j, k] +
KroneckerDelta[i, k]*KroneckerDelta[j, l] ) - ν/
EYoung KroneckerDelta[i, j]*KroneckerDelta[k, l], {i, 1, 3}, {j,
1, 3}, {k, 1, 3}, {l, 1, 3}];
Loading on material is cyclic. One cycle of this loading is enough for curve fitting problem:
LoadingHistory = Table[ 1/2*(Sin[x + 3 π/2] + 1), {x, 0, 2 π, 2 π/80.0}];
Three vectors save the evolution of some variables during each increment of loading. The last one is damage which is what we are looking for:
ϵμp = Table[0, {i, 1, 3}, {j, 1, 3}, {k, 1, Length[LoadingHistory]}];
Xμ = Table[0, {i, 1, 3}, {j, 1, 3}, {k, 1, Length[LoadingHistory]}];
Damage = Table[0, {k, 1, Length[LoadingHistory]}];
Some function are defined to be used:
Devaitoric[tensor_] := Module[{H, tensorD}, H = 1.0/3 Tr[tensor];
tensorD = tensor - 1.0/3 Tr[tensor] IdentityMatrix[3]];
MisesNorm[tensor_] := (
tensorD = Devaitoric[tensor];
Sqrt[3.0/2 Total[tensorD*tensorD, 2]] );
macaulay[x_] := (Abs[x] + x)/2.0;
The main function, solves the set of equations for one step. t1 is the current step where we have the values and t2 is the next step where we are looking to find the values.
mainfun[σmeso_, smat_, S_, t1_, t2_] := (
σ = LoadingHistory[[t2]]*( {
{0.0, 0, 0},
{0, σmeso, 0},
{0, 0, 0}} );
(* meso strain in integration point -tensor inner product*)
ϵ =
Table[Sum[
CElastic[[i, j, k, l]]*σ[[k, l]], {k, 1, 3}, {l, 1,
3}], {i, 1, 3}, {j, 1, 3}];
(*Total Strain at micro level *)
ϵμ = ϵ + β*ϵμp[[All, All,
t1]];
ϵμe = ϵμ - ϵμp[[All, All,
t1]];
σμ =
Table[Sum[
DElastic[[i, j, k, l]]*ϵμe[[k, l]], {k, 1, 3}, {l,
1, 3}], {i, 1, 3}, {j, 1, 3}];
s = σμ - Xμ[[All, All, t1]];
seq = MisesNorm[s];
Δp = 0.0;
f = seq - σf;
If[f > -0.001,
\[ScriptCapitalG] = 3 G (1 - β) + Cy (1 - Damage[[t1]]);
Cs = Table[1, {i, 3}, {j, 3}];
Cp = 1;
While[ Total[Cs, 2] > 1.0*^-13 ∨ Abs[Cp] > 1.0*^-13,
mμ = 3/(2 MisesNorm[s])*Devaitoric[s];
Rs =
s + 2/3 \[ScriptCapitalG] mμ Δp -
Table[Sum[
DElastic[[i, j, k, l]]*ϵ[[k, l]], {k, 1, 3}, {l, 1,
3}], {i, 1, 3}, {j, 1, 3}] +
2 G (1 - β) ϵμp[[All, All, t1]] +
Xμ[[All, All, t1]];
Rp = MisesNorm[s] - σf;
Cp = (Rp - Total[ mμ* Rs, 2])/\[ScriptCapitalG];
Cs = 2/3 ( Total[ mμ* Rs, 2] - Rp) mμ - (Rs *MisesNorm[s] + 2/3 \ [ScriptCapitalG]*Δp*Total[mμ* Rs, 2] mμ)/(MisesNorm[s] + \[ScriptCapitalG] Δp);
s = s + Cs;
Δp = Δp + Cp;
];
(*Updating variables to be used in next increment*)
ϵμp[[All, All, t2]] = ϵμp[[All, All, t1]] + Δp mμ;
Xμ[[All, All, t2]] =2/3 Cy (1 - Damage[[t1]]) ϵμp[[All, All, t2]];
σμ = s + Xμ[[All, All, t2]];
σplus = Select[Eigenvalues[σμ], # >= 0 &];
σminus = Select[Eigenvalues[σμ], # < 0 &];
Yμ = (1 + ν)/(2 EYoung) (Total[ σplus*σplus] + h ( (1 - Damage[[t1]])/(1 - h Damage[[t1]]))^2 Total[ \σminus*σminus]) - ν/(2 EYoung) (macaulay[ Tr[ σμ]]^2 +h ( (1 - Damage[[t1]])/(1 - h Damage[[t1]]))^2 macaulay[ -Tr[ σμ]]^2);
Damage[[t2]] = Damage[[t1]] + (Yμ/S )^smat Δp;
,
Damage[[t2]] = Damage[[t1]];
ϵμp[[All, All, t2]] = ϵμp[[All, All, t1]];
Xμ[[All, All, t2]] = Xμ[[All, All, t1]];
];
Damage
);
Applying this function to the whole loading history will give us the solution over the whole loading:
Do[mainfun[σmeso, smat, S, i, i + 1], {i, 1,
Length[LoadingHistory] - 1}]
Damage evolution and final value of damage:
ListPlot[Damage]
Last@Damage
The experimental data that we are trying to fit this model to are as follows. DExperiment is the damage after one cycle of loading, which we have tried to find above:
σmesovalues = {32.39618739811888300000,
32.88596552115932800000, 33.41106740534313200000,
33.97317756735912300000, 34.57921776339298200000,
35.23385123371712500000, 35.95417637218495300000,
36.74066188224332300000, 37.60823476195416500000,
38.57504805727071300000, 39.66470304378003000000,
40.96920524497549800000, 42.45211121685638000000};
DExperiment = {1.99976*10^-6, 4.08354*10^-6, 7.67889*10^-6,
0.00001303, 0.0000212061, 0.0000310106, 0.0000452827, 0.0000633068,
0.0000858743, 0.000110066, 0.000135709, 0.000169626, 0.000196569};
