How do I incorporate PDEmodel for uncoupled thermo-elastic analysis?

I’ve been trying to solve the Navier-Lame equation for thermoelasticity in a reinforced concrete beam by incorporating PDEmodel. It’s inconvinient, but I want my Young’s modulus to be a function of temperature. My procedure is solving a thermal transient problem and then taking the obtained temperature field into my c and f coefficients. The geometry is created via multicuboid, then transformed into fegeometry and fed to the thermal transient model. Then the same geometry (before conversion to fegeometry) is assigned to the general pde model (createpde(3)). The mesh is seperate between models.
The solution I get is nonsense: displacements for a 2m beam, fixed on both ends and subjected to temperature gradient of 1500K/m is around 1.2m (looks like torsion). I cannot find any mathematical errors. The attached code shows my definition of both c and f coefficients (it is just an overview). Can I use the interpolateTemperature and evaluateTemperatureGradient effectively inside my coefficients? Perhaps there is another problem with my approach or some PDE Toolbox limitations?
Thanks for any potential help on this.
Best regards
KS
%Body of the structural problem, Rt are the thermal transient reults
model = createpde(3);
model.Geometry = g;
generateMesh(model, ‘Hmax’, 0.01);

ccoeff1 = @(location,state) myfunWithAdditionalArgs1(location,state, Rt, length(tlist));
ccoeff2 = @(location,state) myfunWithAdditionalArgs2(location,state, Rt, length(tlist));
fcoeff1 = @(location,state) myfunWithAdditionalArgs3(location, state, Rt, length(tlist));
fcoeff2 = @(location,state) myfunWithAdditionalArgs4(location,state, Rt, length(tlist));
gcoef = @(location,state) myfunWithAdditionalArgs5(location,state, Rt, length(tlist));

specifyCoefficients(model,’m’,0,’d’,0,’c’, ccoeff1,’a’,0,’f’,[0;0;0], ‘Cell’, 1);
specifyCoefficients(model,m=0,d=0,c=ccoeff2,a=0,f=[0;0;0], Cell=[2, 3]);

applyBoundaryCondition(model, "neumann", g=gcoef, Face=[3, 4, 5, 6]);
applyBoundaryCondition(model,"dirichlet",Face=2,u=[0,0,0]);
applyBoundaryCondition(model,"dirichlet",Face=1,u=[0,0,0]);

result = solvepde(model)

%Defining of the coefficients-just the concrete
ccoeff1 = @(location,state) myfunWithAdditionalArgs1(location,state, Rt, length(tlist));
function cmatrix1 = myfunWithAdditionalArgs1(location, state, T, id)
n1 = 45;
nr = numel(location.x);
cmatrix1 = zeros(n1, nr);
E = 30e9.*(1-0.0001*interpolateTemperature(T, location.x, location.y, location.z, id));
nu = 0.2;
mu = E./(2.*(1+nu));
lambda = E.*nu./((1-2.*nu).*(1+nu));

cmatrix1 (1,:) = 2.*mu+lambda;
cmatrix1 (3,:) = mu;
cmatrix1 (6,:) = mu;
cmatrix1 (8,:) = mu;
cmatrix1 (10,:) = lambda;
cmatrix1 (16,:) = mu;
cmatrix1 (18,:) = 2.*mu+lambda;
cmatrix1 (21,:) = mu;
cmatrix1 (24,:) = mu;
cmatrix1 (28,:) = lambda;
cmatrix1 (36,:) = mu;
cmatrix1 (38,:) = lambda;
cmatrix1 (40,:) = mu;
cmatrix1 (42,:) = mu;
cmatrix1 (45,:) = 2.*mu+lambda;
end

fcoeff1 = @(location,state) myfunWithAdditionalArgs3(location, state, Rt, length(tlist));
function fmatrix1 = myfunWithAdditionalArgs3(location, state, T, id)
n1=3;
nr = numel(location.x);
fmatrix1 = zeros(n1, nr);
alpha = 1.2e-5;
nu = 0.2;
E =30e9.*(1-0.0001*interpolateTemperature(T, location.x, location.y, location.z, id));
[gradTx, gradTy, gradTz] = evaluateTemperatureGradient(T, location.x, location.y, location.z, id);
fmatrix1(1,:) = E.*alpha.*gradTx./(1-2.*nu);
fmatrix1(2,:) = E.*alpha.*gradTy./(1-2.*nu);
fmatrix1(3,:) = E.*alpha.*gradTz./(1-2.*nu);
endI’ve been trying to solve the Navier-Lame equation for thermoelasticity in a reinforced concrete beam by incorporating PDEmodel. It’s inconvinient, but I want my Young’s modulus to be a function of temperature. My procedure is solving a thermal transient problem and then taking the obtained temperature field into my c and f coefficients. The geometry is created via multicuboid, then transformed into fegeometry and fed to the thermal transient model. Then the same geometry (before conversion to fegeometry) is assigned to the general pde model (createpde(3)). The mesh is seperate between models.
The solution I get is nonsense: displacements for a 2m beam, fixed on both ends and subjected to temperature gradient of 1500K/m is around 1.2m (looks like torsion). I cannot find any mathematical errors. The attached code shows my definition of both c and f coefficients (it is just an overview). Can I use the interpolateTemperature and evaluateTemperatureGradient effectively inside my coefficients? Perhaps there is another problem with my approach or some PDE Toolbox limitations?
Thanks for any potential help on this.
Best regards
KS
%Body of the structural problem, Rt are the thermal transient reults
model = createpde(3);
model.Geometry = g;
generateMesh(model, ‘Hmax’, 0.01);

specifyCoefficients(model,’m’,0,’d’,0,’c’, ccoeff1,’a’,0,’f’,[0;0;0], ‘Cell’, 1);
specifyCoefficients(model,m=0,d=0,c=ccoeff2,a=0,f=[0;0;0], Cell=[2, 3]);

applyBoundaryCondition(model, "neumann", g=gcoef, Face=[3, 4, 5, 6]);
applyBoundaryCondition(model,"dirichlet",Face=2,u=[0,0,0]);
applyBoundaryCondition(model,"dirichlet",Face=1,u=[0,0,0]);

result = solvepde(model)

fcoeff1 = @(location,state) myfunWithAdditionalArgs3(location, state, Rt, length(tlist));
function fmatrix1 = myfunWithAdditionalArgs3(location, state, T, id)
n1=3;
nr = numel(location.x);
fmatrix1 = zeros(n1, nr);
alpha = 1.2e-5;
nu = 0.2;
E =30e9.*(1-0.0001*interpolateTemperature(T, location.x, location.y, location.z, id));
[gradTx, gradTy, gradTz] = evaluateTemperatureGradient(T, location.x, location.y, location.z, id);
fmatrix1(1,:) = E.*alpha.*gradTx./(1-2.*nu);
fmatrix1(2,:) = E.*alpha.*gradTy./(1-2.*nu);
fmatrix1(3,:) = E.*alpha.*gradTz./(1-2.*nu);
end I’ve been trying to solve the Navier-Lame equation for thermoelasticity in a reinforced concrete beam by incorporating PDEmodel. It’s inconvinient, but I want my Young’s modulus to be a function of temperature. My procedure is solving a thermal transient problem and then taking the obtained temperature field into my c and f coefficients. The geometry is created via multicuboid, then transformed into fegeometry and fed to the thermal transient model. Then the same geometry (before conversion to fegeometry) is assigned to the general pde model (createpde(3)). The mesh is seperate between models.
The solution I get is nonsense: displacements for a 2m beam, fixed on both ends and subjected to temperature gradient of 1500K/m is around 1.2m (looks like torsion). I cannot find any mathematical errors. The attached code shows my definition of both c and f coefficients (it is just an overview). Can I use the interpolateTemperature and evaluateTemperatureGradient effectively inside my coefficients? Perhaps there is another problem with my approach or some PDE Toolbox limitations?
Thanks for any potential help on this.
Best regards
KS
%Body of the structural problem, Rt are the thermal transient reults
model = createpde(3);
model.Geometry = g;
generateMesh(model, ‘Hmax’, 0.01);

specifyCoefficients(model,’m’,0,’d’,0,’c’, ccoeff1,’a’,0,’f’,[0;0;0], ‘Cell’, 1);
specifyCoefficients(model,m=0,d=0,c=ccoeff2,a=0,f=[0;0;0], Cell=[2, 3]);

applyBoundaryCondition(model, "neumann", g=gcoef, Face=[3, 4, 5, 6]);
applyBoundaryCondition(model,"dirichlet",Face=2,u=[0,0,0]);
applyBoundaryCondition(model,"dirichlet",Face=1,u=[0,0,0]);

result = solvepde(model)

Cart

Cart