This started as a comment to Demetri's answer, but it got too long so I made it a proper answer.
If you are trying to get non-linear effects for the time variable, you should model it non-linearly. Fixed effects as you propose will help you with that perfectly since each year will have its orthogonal contribution as a parameter. But if you try to predict an out of time sample (a data point in which the year is not among the ones in the training data) that point will not have a corresponding year in the model matrix. Remember, if you have a categorical variable with 3 categories, and say 4 samples/observations, this is how you write the model matrix, where the $A$ category is the reference category:
$$
\begin{bmatrix}
A \\ B \\ C \\ C
\end{bmatrix}_{\text{samples}} =
\begin{bmatrix}
0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 0
\end{bmatrix}_{\text{model matrix}}
$$
In your approach, each category is a year, i.e. a fixed effect. Now, what do you do if a new sample has category $D$? One of two things happen, you either don't have a prediction or you erroneously give it no parameter or the last parameter in.
If your focus is on in-sample prediction, this isn't a problem at all, and I argue that, if you have large enough data, this is the best approach. Even if the standard errors are large (which indicates that this information, year in our case, is noisy and most likely other variables better estimate the variation in the outcome variable) their estimates can be useful in a prediction setting.
However, if out of time sample estimation is a goal or if you can't spare the degrees of freedom (another way of saying you don't have enough data) a continuous modeling approach is the way to go. By continuous we don't mean linear. There are many ways of doing so, the simplest of them is a simple polynomial transformation of the year variable. Another way is through splines which can come in many flavors (natural cubic splines, B-splines, etc...), and can be included in your model via the splines or mgcv packages in R. This class of modeling techniques is also called GAM or generalized additive models in case you'd like to research about it.
I just saw your comment about the years in the prediction data set. Your problem is solved by splines. To use splines (without penalization/ regularization) with the base-r glm function, you'll need the splines package which comes installed with R. This is how to fit the model:
library(splines)
db_tr <- data[data$group=="1",]
db_pr <- data[data$group=="2",]
f <- formula(y ~ x1 + x2 + x2 + bs(year))
logit.training <- glm(f,family=binomial(link="logit"), data= db_tr)
logit.predict <- predict(logit.training, newdata= db_pr, type="response")
bs stands for B-splines and is the R function to create the B-spline transformation of a variable passed to it in the model matrix created by the formula.
However, the preferred way to estimate those types of models is through a gam in the mgcv package. Here you can find more about Gams.
