SDT-base Contents   Functions      PDF Index |
Purpose
Computation of time and non linear responses.
Syntax
def=fe_time(model) def=fe_time(TimeOpt,model) [def,model,opt]=fe_time(TimeOpt,model) model=fe_time('TimeOpt...',model) TimeOpt=fe_time('TimeOpt...')
Description
fe_time groups static non-linear and transient solvers to compute the response of a FE model given initial conditions, boundary conditions, load case and time parameters. Note that you may access to the fe_time commands graphically with the simulate tab of the feplot GUI. See tutorial (section 4.10) on how to compute a response.
Three types of time integration algorithm are possible: the Newmark schemes, the Theta-method, and the time Discontinuous Galerkin method. Implicit and explicit methods are implemented for the Newmark scheme, depending on the Newmark coefficients β and γ, and non linear problems are supported.
The parameters of a simulation are stored in a time option data structure TimeOpt given as input argument or in a model.Stack entry info,TimeOpt. Initial conditions are stored as a curve,q0 entry.
The solvers selected by the string TimeOpt.Method are
Here is a simple example to illustrate the common use of this function.
model=fe_time('demo bar'); % build the model % Integrated use of TimeOpt model=fe_time('TimeOptSet Newmark .25 .5 0 1e-4 100',model); def=fe_time(model); % compute the response % Define time options and use structure directly opt=fe_time('TimeOpt Newmark .25 .5 0 1e-4 100'); def=fe_time(opt,model); % compute the response % Store as model.Stack entry {'info','TimeOpt',opt} model=stack_set(model,'info','TimeOpt',opt); def=fe_time(model); % compute the response
The TimeOpt data structure has fields to control the solver
getpref('OpenFEM','THRESHOLD',1e-6);
to control the output
This section details the applicable input and the output options.
Initial conditions may be provided in a model.Stack entry of type info named q0 or in an input argument q0. q0 is a data structure containing def and DOF fields as in a FEM result data structure (section 4.10). If any, the second column gives the initial velocity. If q0 is empty, zero initial conditions are taken. In this example, a first simulation is used to determine the initial conditions of the final simulation.
model=fe_time('demo bar'); TimeOpt=fe_time('TimeOpt Newmark .25 .5 0 1e-4 100'); TimeOpt.NeedUVA=[1 1 0]; % first computation to determine initital conditions def=fe_time(TimeOpt,model); % no input force model=fe_case(model,'remove','Point load 1'); % Setting initial conditions q0=struct('def',[def.def(:,end) def.v(:,end)],'DOF',def.DOF); model=stack_set(model,'curve','q0',q0); def=fe_time(TimeOpt,model);
An alternative call is possible using input arguments
def=fe_time(TimeOpt,model,Case,q0)
In this case, it is the input argument q0 which is used instead of an eventual stack entry.
You may define the time dependence of a load using curves as illustrated in section 7.9.
You may specify the time steps by giving the 'TimeVector'
TimeOpt=struct('Method','Newmark','Opt',[.25 .5 ],... 'TimeVector',linspace(0,100e-4,101));
This is useful if you want to use non constant time steps. There is no current implementation for self adaptive time steps.
To illustrate the output options, we use the example of a 2D propagation. Note that this example also features a time dependent DofLoad excitation (see fe_case) defined by a curve, (see fe_curve), here named Point load 1.
model=fe_time('demo 2d'); TimeOpt=fe_time('TimeOpt Newmark .25 .5 0 1e-4 50');
You may specify specific output by selecting DOF indices as below
i1=fe_case(model,'GettDof'); i2=feutil('findnode y==0',model)+.02; TimeOpt.OutInd=fe_c(i1,i2,'ind'); model=stack_set(model,'info','TimeOpt',TimeOpt); def=fe_time(model); % Don't animate this (only bottom line)
You may select specific output time step using TimeOpt.OutputFcn as a vector
TimeOpt.OutputFcn=[11e-4 12e-4]; % opt.OutputFcn='tout=t(1:100:opt.Opt(5));'; % other example TimeOpt=feutil('rmfield',TimeOpt','OutInd'); model=stack_set(model,'info','TimeOpt',TimeOpt); def=fe_time(model); % only two time steps saved
or as a string to evaluate. In this case it is useful to know the names of a few local variables in the fe_time function.
In this example the default output function (for TimeOpt.NeedUVA=[1 1 1]) is used but specified for illustration
TimeOpt.OutputFcn=['out.def(:,j1+1)=u;' ... 'out.v(:,j1+1)=v;out.a(:,j1+1)=a;']; model=stack_set(model,'info','TimeOpt',TimeOpt); def=fe_time(model); % full deformation saved
This example illustrates how to display the result (see feplot) and make a movie
cf=feplot(model,def); fecom('ColorDataEvalA'); fecom(cf,'SetProp sel(1).fsProp','FaceAlpha',1,'EdgeAlpha',0.1); cf.ua.clim=[0 2e-6];fecom(';view2;AnimTime;ch20;scd1e-2;'); st=fullfile(sdtdef('tempdir'),'test.gif'); fecom(['animMovie ' st]);fprintf('\nGenerated movie %s\n',st);
Note that you must choose the Anim:Time option in the feplot GUI.
You may want to select outputs using observations matrix
model=fe_time('demo bar'); Case=fe_case('gett',model); i1=feutil('findnode x>30',model); TimeOpt=fe_time('TimeOpt Newmark .25 .5 0 1e-4 100'); TimeOpt.c_u=fe_c(Case.DOF,i1+.01); % observation matrix TimeOpt.lab_u=fe_c(Case.DOF,i1+.01,'dofs'); % labels def=fe_time(TimeOpt,model);
If you want to specialize the output time and function you can specify the SaveTimes as a time vector indicating at which time the SaveFcn string will be evaluated. A typical TimeOpt would contain
TimeOpt.SaveTimes=[0:Ts:TotalTime];
TimeOpt.SaveFcn='My_function(''Output'',u,v,a,opt,out,j1,t);';
The field FinalCleanupFcn of the TimeOpt can be used to specify what is done just after the time integration.
fe_simul provides a generic clean up function which can be called using
opt.FinalCleanupFcn='fe_simul(''fe_timeCleanup'',model)';
If the output has been directly saved or from iiplot it is possible to load the results with the same display options than for the fe_timeCleanup using fe_simul('fe_timeLoad',fname)';
Some command options can be used:
-restitFeplot adds a .TR field to the output to allow deformation on the fly restitution in feplot. These two options cannot be specified simultaneously.
For the Newmark scheme, TimeOpt has the form
TimeOpt=struct('Method','Newmark','Opt',Opt)
where TimeOpt.Opt is defined by
[beta gamma t0 deltaT Nstep]
beta and gamma are the standard Newmark parameters [46] ([0 0.5] for explicit and default at [.25 .5] for implicit), t0 the initial time, deltaT the fixed time step, Nstep the number of steps.
The default residual is r = (ft(j1,:)*fc'-v'*c-u'*k)'; (notice the sign change when compared to NLNewmark).
This is a simple 1D example plotting the propagation of the velocity field using a Newmark implicit algorithm. Rayleigh damping is declared using the info,Rayleigh case entry.
model=fe_time('demo bar'); data=struct('DOF',2.01,'def',1e6,... 'curve',fe_curve('test ricker dt=1e-3 A=1')); model = fe_case(model,'DOFLoad','Point load 1',data); TimeOpt=struct('Method','Newmark','Opt',[.25 .5 3e-4 1e-4 100],... 'NeedUVA',[1 1 0]); def=fe_time(TimeOpt,model); % plotting velocity (propagation of the signal) def_v=def;def_v.def=def_v.v; def_v.DOF=def.DOF+.01; feplot(model,def_v); if sp_util('issdt'); fecom(';view2;animtime;ch30;scd3'); else; fecom(';view2;scaledef3'); end
The time discontinuous Galerkin is a very accurate time solver [59] [60] but it is much more time consuming than the Newmark schemes. No damping and no non linearities are supported for Discontinuous Galerkin method.
The options are [unused unused t0 deltaT Nstep Nf], deltaT is the fixed time step, Nstep the number of steps and Nf the optional number of time step of the input force.
This is the same 1D example but using the Discontinuous Galerkin method:
model=fe_time('demo bar'); TimeOpt=fe_time('TimeOpt DG Inf Inf 0. 1e-4 100'); TimeOpt.NeedUVA=[1 1 0]; def=fe_time(TimeOpt,model); def_v=def;def_v.def=def_v.v; def_v.DOF=def.DOF+.01; feplot(model,def_v); if sp_util('issdt'); fecom(';view2;animtime;ch30;scd3'); ... else; fecom(';view2;scaledef3'); end
For the non linear Newmark scheme, TimeOpt has the same form as for the linear scheme (method Newmark). Additional fields can be specified in the TimeOpt data structure
Jacobian | string to be evaluated to generate a factored
jacobian matrix in matrix or ofact object ki. The default
jacobian matrix is 'ki=ofact(model.K{3}+2/dt*model.K{2} +4/(dt*dt)*model.K{1});' |
Residual | Defines the residual used for the Newton iterations of each type step. It is typically a call to an external function. The default residual is 'r = model.K{1}*a+model.K{2}*v+model.K{3}*u-fc;' where fc is the current external load, obtained using (ft(j1,:)*fc')' at each time step. |
IterInit | evaluated when entering the correction iterations. This can be used to initialize tolerances, change mode in a co-simulation scheme, etc. |
IterEnd | evaluated when exiting the correction iterations. This can be used to save specific data, ... |
IterFcn | Correction iteration algorithm function, available are iterNewton (default when omitted) or iterNewton_Sec. Details of the implementation are given in the staticNewton below. |
MaxIterSec | for iterNewton_Sec applications (see staticNewton). |
ResSec | for iterNewton_Sec applications (see staticNewton). |
Following example is a simple beam, clamped at one end, connected by a linear spring at other end and also by a non linear cubic spring. The NL cubic spring is modeled by a load added in the residual expression.
% Get simple test case for NL simulation in sdtweb demosdt('BeamEndSpring') model=demosdt('BeamEndSpring'); % simple example building opt=stack_get(model,'info','TimeOpt','GetData'); disp(opt.Residual) opt.Follow=1; % activate simple monitoring of the % number of Newton iterations at each time step def=fe_time(opt,model);
For non linear static problems, the Newton solver iterNewton is used. TimeOpt has a similar form as with the NLNewmark method but no parameter Opt is used.
An increment control algorithm iterNewton_Sec can be used when convergence is difficult or slow (as it happens for systems showing high stiffness variations). The Newton increment Δ q is then the first step of a line search algorithm to optimize the corrective displacement increment ρ Δ q, ρ∈ R in the iteration. This optimum is found using the secant iteration method. Only a few optimization iterations are needed since this does not control the mechanical equilibrium but only the relevance of the Newton increment. Each secant iteration requires two residual computations, which can be costly, but more efficient when a large number of standard iterations (matrix inversion) is required to obtain convergence.
Fields can be specified in the TimeOpt data structure
Jacobian | defaults to 'ki=ofact(model.K{3});' |
Residual | defaults to 'r = model.K{3}*u-fc;' |
IterInit | and IterEnd are supported see fe_time TimeOpt |
IterEnd | |
MaxIterSec | Maximum secant iterations for the iterNewton_Sec iteration algorithm. The default is 3 when omitted. |
ResSec | Residual evaluation for the secant iterations of the iterNewton_Sec iteration algorithm. When omitted, fe_timetries to interpret the Residual field. The function must fill in the secant residual evaluation r1 which two columns will contain the residual for solution rho(1)*dq and rho(2)*dq. The default ResSec field will be then 'r1(:,1) = model.K{3}*(u-rho(1)*dq)-fc; r1(:,2) = model.K{3}*(u-rho(2)*dq)-fc;'. |
Below is a demonstration non-linear large transform statics.
% Sample mesh, see script with sdtweb demosdt('LargeTransform') model=demosdt('largeTransform'); % % Now perform the Newton loop model=stack_set(model,'info','TimeOpt', ... struct('Opt',[],'Method','StaticNewton',... 'Jacobian','ki=basic_jacobian(model,ki,0.,0.,opt.Opt);',... 'NoT',1, ... % Don't eliminate constraints in model.K 'AssembleCall','assemble -fetimeNoT -cfield1', ... 'IterInit','opt=fe_simul(''IterInitNLStatic'',model,Case,opt);')); model=fe_case(model,'setcurve','PointLoad', ... fe_curve('testramp NStep=20 Yf=1e-6')); % 20 steps gradual load def=fe_time(model); cf=feplot(model,def); fecom(';ch20;scc1;colordataEvalZ'); % View shape ci=iiplot(def);iicom('ch',{'DOF',288.03}) % View response
You may want to use numerical damping in a time integration scheme, the first possibility is to tune the Newmark parameters using a coefficient α such that β = (1+α)2/4 and γ=1/2+α. This is known to implement too much damping at low frequencies and is very depending on the time step [46].
A better way to implement numerical damping is to use the HHT-α method which applies the Newmark time integration scheme to a modified residual balancing the forces with the previous time step.
For the HHT-α scheme, TimeOpt has the form
TimeOpt=struct('Method','nlnewmark','Opt',Opt,... 'HHTalpha',alpha)
where TimeOpt.Opt is defined by
[beta gamma t0 deltaT Nstep]
beta and gamma are the standard Newmark parameters [46] with numerical damping, t0 the initial time, deltaT the fixed time step, Nstep the number of steps.
The automatic TimeOpt generation call takes the form [alpha unused t0 deltaT Nstep] and will compute the corresponding β, γ parameters.
This is a simple 1D example plotting the propagation of the velocity field using the HHT-α implicit algorithm:
model=fe_time('demo bar'); TimeOpt=fe_time('TimeOpt hht .05 Inf 3e-4 1e-4 100'); TimeOpt.NeedUVA=[1 1 0]; def=fe_time(TimeOpt,model);
The call
TimeOpt=fe_time('TimeOpt hht .05 Inf 3e-4 1e-4 100');
is strictly equivalent to
TimeOpt=struct('Method','nlnewmark',... 'Opt',[.275625 .55 3e-4 1e-4 100],... 'HHTalpha',.05);
The θ-method is a velocity based solver, whose formulation is given for example in [61, 62]. It considers the acceleration as a distribution, thus relaxing discontinuity problems in non-smooth dynamics. Only a linear implementation is provided in fe_time. The user is nevertheless free to implement a non-linear iteration, through his own IterFcn.
This method takes only one integration parameter for its scheme, θ set by default at 0.5. Any values between 0.5 and 1 can be used, but numerical damping occurs for θ>0.5.
The TimeOpt.Opt takes the form [theta unused t0 deltaT Nstep].
This is a simple 1D example plotting the propagation of the velocity field using the θ-Method:
model=fe_time('demo bar'); TimeOpt=fe_time('TimeOpt theta .5 0 3e-4 100'); def=fe_time(TimeOpt,model);
This method can be used to integrate first order problem of the form . One can use it to solve transient heat diffusion equation (see p_heat).
Integration scheme is of the form
θ can be define in opt.Opt(1).
Explicit Euler (θ=0) is not implemented at this time.
Best accuracy is obtained with θ=1/2 (Crank-Nicolson).
See also